Numerical Simulation of Bone Remodeling Coupling the Damage Repair Process in Human Proximal Femur

: Microdamage is produced in bone tissue under the long-term effects of physiological loading, as well as age, disease and other factors. Bone remodeling can repair microdamage, otherwise this damage will undermine bone quality and even lead to fractures. In this paper, the damage variable was introduced into the remodeling algorithm.The new remodeling algorithm contains a quadratic term that can simulate reduction in bone density after large numbers of loading cycles. The model was applied in conjunction with the 3D finite element method (FEM) to the remodeling of the proximal femur. The results showed that the initial accumulation of fatigue damage led to an increase in density but when the damage reached a certain level, the bone density decreased rapidly until the femur failed. With the accumulation of damage, bone remodeling was coupled with fatigue damage to maintain the function of bone. When the accumulation of damage reached a certain level, bone remodeling failed to repair the accumulated fatigue damage in time, and continued cyclic loading significantly weakened the loadbearing capacity of the bone. The new mathematical model not only predicts fatigue life, but also helps to further understand the compromise between damage repair and damage accumulation, which is of great significance for the prevention and treatment of clinical bone diseases.


Introduction
Human bone is a living tissue whose structure and shape are constantly being renewed to maintain its mechanical performance and adapt to its changing mechanical environment [1]. Bone is often subjected to repetitive fatigue loading during daily physical activities. When the loading exceeds a certain stress-strain level, microdamage or microcracks will occur in the bone tissue and with sufficient accumulation of damage may cause fracture [2]. However, bone tissue, as a special biological material, can repair itself in order to adapt to the surrounding mechanical and physiological environment. It does this through bone resorption and bone formation at specific locations to adjust its internal structure and external shape [3,4]. This allows bone to maintain or Physiologically, bone remodeling is a process of bone tissue repair that employs basic multicellular units (BMUs) [11][12][13]. BMUs are discrete anatomical structures in which osteoclasts remove damaged tissue and osteoblasts generate new bone tissue. Fatigue damage first activates the BMUs and eliminates excess damage through the bone remodeling process [14]. Burr [15] proposed, based on spatial and temporal associations with new remodeling, that accumulation of damage in bone could initiate the remodeling process to repair itself. In other words, the bone can sense the damaged area and repair it before fracture occurs [16,17].
Damage plays an important role in the bone remodeling process and may involve substantial changes in mechanical function so it is very important to investigate the relationship of bone remodeling to fatigue damage. In this study, a new damage adaptive remodeling model is proposed. The FE numerical simulation was used to obtain the corresponding damage variable by changing the number of loading cycles, and bone remodeling under different degrees of damage was evaluated. The resorption of damage associated with excessive loading cycles also was simulated. The new mathematical model established in this paper can predict the fatigue life of the femur, and also provides a reference for the prevention and treatment of clinical fractures.

Development of Bone Remodeling Algorithm
Most relevant previous studies on bone remodeling were experimental and analyzed the mechanical properties of bone through in vivo experiments. With the rapid development of computer technology in the 1990s, the numerical simulation of bone remodeling using computer simulation technology became more popular. At present, there are mainly two kinds of algorithms for the numerical simulation of bone remodeling: the mechanical model and the physiological model. The mechanical model is generally expressed as the change of apparent density with a mechanical stimulus, whereas the physiological model is generally expressed by the damage state and porosity change [18].
A conceptual mechanical steady-state theory was proposed by Frost. Based on this theory, Carter et al. [19] presented a qualitative description of bone remodeling that could be expressed as where ρ is the apparent density, which represents the internal structural characteristics of bone; S is the mechanical stimulus; B is the remodeling rate coefficient; K is the steady state reference; and ρ cb is the maximum density of cortical bone. Huiskes [20] proposed a remodeling algorithm which considered the concept of a "lazy zone." Weinans [21] also proposed the concept of a "dead zone" and believed that bone remodeling didn't occur when the stimulus was between K(1 ± ω).
As mentioned above, bone damage and repair are coupled, and bone remodeling can also influence the bone's damage behavior. Based on these conditions, models of bone behavior have been developed and applied to the prediction of fracture and implant loosening [1,[22][23][24][25][26]. The bone damage-repair process is shown in Fig. 1. Under the long-term action of external loading, bone will develop fatigue damage, causing changes in bone material properties. Fatigue damage can activate bone remodeling [27], replacing the damaged bone and maintaining its mechanical properties. Based on the process of bone damage-repair, this paper proposes a new bone remodeling model that considers different values of damage accumulation, which is significant for preventing fractures.

Bone Remodeling Algorithm Based on Damage Repair
In this paper, based on the classical bone remodeling equation, the damage variable D proposed by Hambli [28,29] was introduced into the simulation, and a new control equation for fatigue damage was proposed to simulate bone damage behavior. In the formula, N is the number of loading cycles, N f is the number of cycles to failure, and α and β represent the material fatigue parameters of bone. We assigned α = 0.01 and β = 0.001. We considered the damage variable, caused by the number of loading cycles, to be an isotropic parameter. This formula shows the evolution of damage with the number of loading cycles. The scope of the damage variable D is from 0 to 1. D = 0 corresponds to the undamaged state; when the number of loading cycles is equal to the number of failure cycles, D = 1, and there will be local failure. In other words, the damage variable D is positively correlated with the number of loading cycles N. In the simulation, the damage variable can be changed by changing the number of loading cycles. From this, the new control equations are modified on the basis of the traditional control equation and expressed as: In the formulae the remodeling control quantity ρ is the apparent density of bone; S is the mechanical stimulus; A, B are bone material parameters; K is the steady-state reference value; and ω represents the dead zone range. When S is between (1 ± ω)K, bone remodeling will not occur. Based on [21], the specific values of parameters in Eq. (3) are: Using this equation, the rate of change in density with damage accumulation can be studied. The equation also can be used to calculate the change of bone density values, and the new bone structure can be obtained by constantly updating the density value through iterative calculation. The density change rate curve expressed by this equation is a quadratic curve with a downward opening that varies with the damage. As the number of loading cycles increases, which increases damage, bone density and the density change rate first increase to a maximum value (bone gain), then begin to decrease until they become negative where overload resorption may occur [30]. In the equation, by changing the bone material parameters A and B, the shape and position of the quadratic curve can be adjusted to the corresponding mechanical behavior of different types of the bone.
In the study of bone remodeling, equivalent stress, equivalent strain and strain energy density all can be selected as the mechanical stimulus. According to Weinans, bone is regarded as a continuous material [21], and the strain energy density S = U/ρ is selected as the mechanical stimulus to simulate bone remodeling based on the fatigue damage model. The apparent strain energy density is used as the strain energy density in the calculation, regardless of the stress history. Therefore, the strain energy density is expressed by the current stress and strain [30]. In the simulation here, uniaxial loading is considered, then strain energy can be expressed as In previous studies, scholars defined the properties of bone material related to their apparent density [31]. If we consider bone tissue as an isotropic material, the relationship between the elastic modulus and apparent density of bone material can be defined as: E = E(ρ) = Cρ r , where C and r are the constants, and their values generally are obtained experimentally [32,33]. This equation has been accepted and is used widely. In this paper, we considered the role of the damage variable [34], then the elastic modulus E = C(1 − D)ρ γ , and the constants C = 3790 MPa · (g/cm 3 ) −3 and γ = 3, so the mechanical stimulus can be expressed as a function of stress: Next, substituting Eq. (4) into Eq. (3): We used the improved Euler's method to calculate the differential equation, as shown below:

Application to Femoral Remodeling Using the Finite Element Method
To realize the practical value of the algorithm, the new bone remodeling mathematical model was applied to the FE software Abaqus. By converting the remodeling control equation into the Python program and using iSight for software integration, the control equation was combined with the FE calculation process in Abaqus to improve the efficiency of the cyclic calculation.

FE Simulation of 2D Proximal Femur
Due to the large number of elements in the 3D FE model of the proximal femur and complex bone remodeling algorithm, a simpler 2D FE model was used to simulate bone remodeling and to verify the feasibility of the algorithm. Then, the 3D proximal femur model was simulated with bone remodeling considering fatigue damage.
The CT scan data were imported into Mimics software to establish the 2D shape of the proximal femur in the coronal plane, which was imported into Abaqus for mesh division to obtain a 2D FE model. The 2D FE model assumed that the femur is an isotropic material with an initial elastic modulus of 14 GPa, a Poisson's ratio of 0.3, and an initial density of 1.4 g/cm 3 . Three conditions are shown in Tab. 1 [35] and illustrated in Fig. 2. Based on the simulation of bone remodeling under these three typical conditions, the von Mises stress distribution of the femur was obtained over the loading period. The maximum stresses corresponding to different cycles were extracted, and iterative calculations were carried out by using the algorithm proposed in this paper. The damage variable, elastic modulus and density corresponding to different cycles were calculated successively, and their variation trends were analyzed.
As shown in Figs. 3-5, the von Mises stress corresponded to the increase of loading cycles under the three typical loading conditions. By analyzing the stress nephogram under these three conditions, it was found that the maximum stress was concentrated primarily near the femoral shaft. The stress on the periphery of the femoral shaft was greater, and bone remodeling began in this area initially. Through the repair effect of bone remodeling, the corresponding density of the periphery of the femoral shaft increased and the area with high femoral shaft density tended to thicken. The intracortical area of the femoral shaft was less stimulated by stress, so the repair effect of bone remodeling was relatively weak. The stress distribution can only be seen simply through the cloud map, and the observation of the overall trend change is not very obvious. As the amount of damage was positively correlated with the number of cycles, the maximum von Mises stress of different cycles was extracted and brought into the bone remodeling algorithm for iterative calculation. The variation of elastic modulus and density with the damage was obtained and analyzed. Fig. 6 shows the variation of the elastic modulus under the three typical conditions, and the overall trend is downward. In single-leg stance or adduction, as damage increases, the elastic modulus decreases slowly at first and then rapidly. In abduction, the elastic modulus decreases uniformly. In single-leg stance or adduction, the corresponding damage is different when the trend of elastic modulus changes, indicating that the femur presents a different relationship between modulus and damage under different conditions of loading angle and magnitude. So although the overall trend is roughly the same, the corresponding damage is different when there is a turning point. With the accumulation of damage, bone remodeling is insufficient to repair it, so the decline of elastic modulus is more obvious. As shown in Fig. 7, in single-leg stance, when the damage variable was less than 0.8, although the density change rate fluctuated up and down the density change was greater than zero. Similarly, as shown in Figs. 8 and 9, in abduction when the damage was less than 0.7, and in adduction when the damage was less than 0.8, the overall density change rate was greater than zero. A density change rate greater than zero indicates that the fatigue damage caused by the cycles can be repaired by bone remodeling, and bone deposition would occur on the femoral shaft at this time. When the cycles and damage continue to increase, the density change rate begins to decrease until it is less than zero, indicating that bone resorption would occur at this time. Different accumulations of damage correspond to different remodeling behaviors. We believe there is a critical damage quantity D 0 under different loading conditions. When D < D 0 , bone remodeling is coupled with damage and can repair it. When D > D 0 , bone remodeling can no longer fully repair the damage, and so damage continues to accumulate, eventually leading to fatigue failure.   The variation of density under three typical conditions was further analyzed, as shown in Figs. 10-12. When D < D 0 , the bone density increases, consistent with the trend of density change rate, because bone deposition occurs at this time. When D > D 0 , the damage continues to increase. In single-leg stance or adduction, the density has a significant downward trend, indicating that resorption due to overload occurs at this time, and bone density decreases. The variation trend of the density value corresponding to abduction does not show an obvious decrease. This is because fatigue damage can be repaired by bone remodeling, so the time when density decreases occurs later after more cycles.

3D FE Simulation of Human Proximal Femur
We used the 3D FE model of the proximal femur to simulate the macroscopic mechanical properties and damage behavior of cortical bone based on the assumption that cortical bone at the continuum level is a purely elastic material. Since cancellous bone is a solid-liquid twophase material composed of trabeculae and bone marrow, it presents a more complex irregular structural network with high porosity. Therefore, we studied the mechanical properties only of cortical bone. Mimics software and Hypermesh were combined to establish the proximal femur 3D FE model quickly and accurately.  In Abaqus, local constraints were applied to the bottom of the femur model and loading was set based on human muscle strength and joint force; typical conditions were shown in Tab. 1. The loading details are shown in Fig. 13 and are based on an adult of 70 Kg bodyweight when stationary or walking slowly. The ODB result file was imported into iSight to extract the maximum von Mises stress, which was substituted into the bone remodeling algorithm. The mechanical stimulus and bone density change rate were obtained in turn, and the updated density was obtained through the Euler iterative algorithm.
After the simulation was completed, the performance parameters of femur remodeling were analyzed. The results of bone remodeling in single-leg stance were analyzed first followed by analyses in hip joint abduction and adduction.

Figure 12:
The variation trend of density before and after critical damage in adduction Figure 13: Schematic of the 3D proximal femur model with loads and boundary conditions. A CT scan of a male volunteer with normal femurs was performed. A 3D geometric model of the proximal femur was established using Mimics and geomagic software. The mesh was divided in hypermesh to establish a 3D FE model. Loads and constraints were applied in Abaqus, and iSight was used in combination to study bone remodeling under three loading conditions In this bone remodeling algorithm, the damage variable gradually increases with the number of cycles. The corresponding damage was obtained by changing the number of cycles in the simulation, and variations in the elastic modulus of the femur were analyzed (Fig. 14).
Under cyclic loading the elastic modulus of bone changes in three different stages. In the first stage, damage is produced, and the modulus decreases rapidly. The damage continues to increase in the second stage even though some bone remodeling occurs to repair the damage. Therefore, damage is resorbed and energy dissipated in the second stage, but the elastic modulus does not change significantly. The damage continues to accumulate in the third stage. Bone remodeling continues to occur, but is insufficient to repair all the damage which continues to accumulate, leading to a rapid decline in modulus.

Figure 14:
Elastic modulus of bone plotted against the damage variable in single-leg stance Next, a comprehensive analysis of the bone density changes and rate of change with fatigue damage was performed (Figs. 15, 16). Bone remodeling plays a role in the process of repairing damage. Fig. 15 shows that the density change rate first increases sharply with an increase of damage, followed by some repair through bone remodeling. When the density change rate reaches the maximum, the bone repair ability is greater than damage accumulation. After that, the density change rate gradually decreases, indicating that bone remodeling and damage gradually tends toward an equilibrium. When the density change rate is zero, the process of bone remodeling ends. As the damage accumulation continues to increase, the density change rate drops sharply, and bone resorption would be predicted to occur. The change of density values is shown in Fig. 16. With the increase of damage, due to the effect of bone remodeling, density gradually increases to make up for the bone loss caused by damage. This is relatively consistent with the trend of density change rate. However, as the damage continues to accumulate, the density decreases. When the density drops to zero, overload resorption occurs and the femur can fail locally due to the accumulation of damage.

Figure 16:
The variation trend of the density in single-leg stance before and after critical damage The variation of the elastic modulus in abduction and adduction also was analyzed (Fig. 17). It was found that the overall change trend of elastic modulus is similar between these conditions and shows a downward trend. In abduction, the overall elastic modulus shows a uniform downward trend. In adduction, the elastic modulus first presents a steady downward trend. After the damage exceeds 0.8, the elastic modulus decreases more quickly. The density change rate (Fig. 18) and the variation in bone density (Fig. 19) under the two conditions were analyzed. As the damage increases, the variation of density and its change rate in abduction and adduction is similar to single-leg stance. The damage variable D is proportional to the number of loading cycles and the cycle ratio N/N f , that is, the damage variable represents the ratio of loading cycles to failure cycles. It was found that different damage situations correspond to different bone remodeling behaviors. There was a different critical damage quality D 0 in different loading conditions. When D < D 0 , bone remodeling was coupled with the damage. Although the damage increased, bone remodeling could repair some of the damage at this time. When D > D 0 , with the accumulation of damage, bone remodeling can no longer maintain sufficient repair, eventually leading to fatigue failure.

Discussion
In conclusion, FE simulations were implemented to predict the evolution of the bone damage healing process through remodeling by using a model that combined both strain-adaptive and damage-adaptive mechanisms. This work demonstrates the advantage of bone remodeling simulations to predict the healing behavior in the proximal femur. However, several important simplifications have been performed in this model. The orientation of bone material in the proximal femur is difficult to determine, and the relationship between the elastic modulus of bone and its apparent density assumed an isotropic material. In addition, the actual bone remodeling process is not only affected by the external mechanical environment, but also by the internal physiological environment. In future work, the anisotropy and the interaction between bone and blood, bone marrow and cells, as well as the role of metabolism in the process of bone remodeling and repair damage, should be considered to establish a more comprehensive damage model. Further, resorption and formation rates in the algorithm should be experimentally measured. In the future, it will be necessary to improve these models by including some of these effects.

Conclusions
In this paper, the damage variable was introduced into the bone remodeling study, and a new bone remodeling algorithm considering fatigue damage was proposed. The new mathematical model can simulate behavior of bone remodeling in different damage states and also can describe resorption due to damage caused by cyclic loading. The model found that different damage amounts correspond to different bone remodeling behaviors in the proximal femur.