Nonlinear Mathematical Modelling of Bone Damage and Remodelling Behaviour in Human Femur

In this paper, the remodelling model [29, 30] was used to simulate the bone remodelling and damage growth/repair. In this model, the apparent density was used to evaluate the mechanical properties of bone tissues and was considered as a single scalar. The changing rate of bone density is (1) ρ˙=kr˙Svρ^ \dot \rho = k\dot r{S_v}\hat \rho where ρ is the apparent density, which represents the internal structural characteristics of bones. The change in bone density ρ˙ \dot \rho is assumed to be related to the maximum density ρ^ \hat \rho and S_v denotes the internal surface per unit volume. It can be expressed as (2) Sv=0.02876P5-0.10104P4+0.13396P3-0.09304P2+0.03226P {S_v} = 0.02876{P^5} - 0.10104{P^4} + 0.13396{P^3} - 0.09304{P^2} + 0.03226P where P is the porosity of the bone. The relation between the apparent density and porosity is (3) ρ=(1-P)×ρ0 \rho = \left( {1 - P} \right) \times {\rho _0} where ρ₀ is the tissue density, which is a constant; k is the ratio of the available remodelling surface and the total bone internal surface; and r˙ \dot r is the change rate of bone remodelling after considering the dead zone.

Then the damage can be included in the model. First, we define ω as a damage value to evaluate the damage of bones. The damage rate is ω˙(ω˙=1Nf) \dot \omega ( {\dot \omega = {1 \over {{N_f}}}} ) [24] according to Miner's rule, where N_f is the number of cycles to reach failure for the bones at a given stress level. Then it can be calculated using logN_f = H logσⁱ + JT + Kρⁱ + M, where σ is the stress (Mpa), T is the temperature (°C) and ρⁱ is the density of bones (gcm⁻³), respectively [31]. Here K = 2.364, H = −7.789, M = 15.47 and J = −0.0206 are all constants. Then, the damage can be calculated as (3) ω=∫0tω˙dt,ω≤0.98 \omega = \int\limits_0^t \dot \omega dt,\omega \le 0.98 where ω is the local bone to be repaired and t is the remodelling duration.

The damage threshold ω_crit = 2.6 × 10¹⁰ is obtained when the damage reaches 3,500μɛ during one cycle.

If damage ω is less than ω_crit, adaptive remodelling is activated. Strain energy density at location i can be calculated as Ui=εiEiεi2ρi {U^i} = {{{\varepsilon ^i}{E^i}{\varepsilon ^i}} \over {2{\rho ^i}}} , where Eⁱ is the elastic modulus (MPa) at location i (assumed isotropic), ɛⁱ is the strain and the unit of Uⁱ is MPa(gcm⁻³)⁻¹. At a given location i, the stimulus for the element is Sⁱ, defined as Sⁱ = Uⁱ −U_ref.

The changing rate of density is (4) dρ(xj,t)dt=C1Sstrainj {{d\rho \left( {{x^j},t} \right)} \over {dt}} = {C_1}S_{strain}^j where ɛ^ref is dependent on the strain of the element j and ɛ^j is determined based on the following criteria: If the strain of the element ɛ ≤ 1,000μɛ, ɛ^ref = 1,000μɛ and C₁ = 3.87; if the strain1,000μɛ < ɛ ≤ 2,000μɛ, ɛ^ref = 1,000μɛ and C₁ = 0; If the strain 2,000μɛ < ɛ ≤ 3,500μɛ, ɛ^ref = 1,000μɛ and C₁ = 198.

If damage accumulation ω exceeds the critical value ω_crit, the damage process happens.

Thus (5) dρa(xj,t)dt=B∑j=1Nfj(a,xj)Sj(x) {{d{\rho _a}\left( {{x_j},t} \right)} \over {dt}} = B\sum\limits_{j = 1}^N {f^j}\left( {a,{x^j}} \right){S^j}\left( x \right) where ρ_a is the bone density of bone remodelling, parameter B and ρ_max are the proportion coefficients and the maximal bone density. B = 0.05 (gcm⁻³)² MPa⁻¹day⁻¹ and ρ_max = 1.80gcm⁻³. fⁱ (x,xⁱ) function, denoting the decay in signal intensity relative to distance d (x,xⁱ) and decay parameter D. Then E and v can be expressed as (6) E={1,007×ρ2,if 0.01g/cm−3≤ρ<0.25g/cm−3255×ρ,if 0.25g/cm−3≤ρ<0.40g/cm−32,972×ρ2−933×ρ,if 0.40g/cm−3≤ρ<1.20g/cm−31,763×ρ3.25,if 1.20g/cm−3≤ρ≤1.80g/cm−3 E = \left\{ {\matrix{{1,007 \times {\rho ^2},} \hfill & {{\rm{if}}\,0.01{\rm{g/c}}{{\rm{m}}^{ - 3}} \le \rho < 0.25{\rm{g/c}}{{\rm{m}}^{ - 3}}} \hfill \cr {255 \times \rho ,} \hfill & {{\rm{if}}\,0.25{\rm{g/c}}{{\rm{m}}^{ - 3}} \le \rho < 0.40{\rm{g/c}}{{\rm{m}}^{ - 3}}} \hfill \cr {2,972 \times {\rho ^2} - 933 \times \rho ,} \hfill & {{\rm{if}}\,0.40{\rm{g/c}}{{\rm{m}}^{ - 3}} \le \rho < 1.20{\rm{g/c}}{{\rm{m}}^{ - 3}}} \hfill \cr {1,763 \times {\rho ^{3.25}},} \hfill & {{\rm{if}}\,1.20{\rm{g/c}}{{\rm{m}}^{ - 3}} \le \rho \le 1.80{\rm{g/c}}{{\rm{m}}^{ - 3}}} \hfill \cr } } \right. (7) υ={0.2,if 0.01g/cm-3≤ρ<1.20g/cm-30.32,if 1.20g/cm-3≤ρ<1.80g/cm-3 \upsilon = \left\{ {\matrix{{0.2,} \hfill & {{\rm{if}}\,0.01{\rm{g/c}}{{\rm{m}}^{ - 3}} \le \rho < 1.20{\rm{g/c}}{{\rm{m}}^{ - 3}}} \hfill \cr {0.32,} \hfill & {{\rm{if}}\,1.20{\rm{g/c}}{{\rm{m}}^{ - 3}} \le \rho < 1.80{\rm{g/c}}{{\rm{m}}^{ - 3}}} \hfill \cr } } \right.

A human femur bone model was used to predict the change in bone density and damage under different loadings by finite element method (FEM) (Figure 2) to show how the micro-damage-based mechanism works in the clinical practice.

Fig. 2

As is shown in Figure 2, a four-node quadrilateral plane strain element is used to generate mesh. It is assumed that bones are compressible isotropic materials. Therefore, Young's modulus, permeability and Poisson's ratio are preliminarily defined, all by mineral density. These parameters are all changing with time.

Bone tissue is affected by repeated loads during daily activities, and the peak value of cyclic loads is dominant in the remodelling process [32]. Therefore, the static force and fixed constraint are set as the boundary conditions of the model.

Based on the mechanical law, a remodelling and damage model of human proximal femur is established, which takes into account the strain energy and stress of each element under usual stimulation. In the initial model, bone remodelling was considered and the density distribution was uniform. The biomechanical stimulation of each unit was calculated and the bone mineral density (BMD) of each unit was updated. After 1,000 iterative steps, we can clearly see that two high-density cortical layers appear around the model and a low-density area corresponds to the medullary cavity. The density distribution of the femoral head is very complicated. Two high-density areas and two low-density ones can be seen in and around the femoral head. Then the damage mechanism is introduced into the simulation model. The damage region is taken as the initial, and continuous simulation is carried out using the controlled boundary conditions. The initial damage is a constant higher than the damage threshold to activate the damage repair mechanism.

In addition, we divided the bone remodelling process into three stages: absorption period (20 days), recovery (10 days) and formation (90 days) [33]. During each cycle, the FEM was used to calculate and obtain the stresses and strains in the femoral model.

The model of bone remodelling is implemented using the finite element software ANSYS, and each element is regarded as a separate model. Bone remodelling algorithms can be summarized as follows. First, assume the initial properties of the material and apply the boundary condition. The stresses and strains can be obtained according to the formula. Then the mechanical signal of each sensor can be calculated in bones. It will be detected by the tissues and these signals are sent to initiate the bone remodeling. According to the magnitude of these signals, the change in cell density is simulated and the new density can then be calculated. For the element whose density has changed, a new elastic modulus is determined in the next step.

In this case, all the material parameters comply with the above model and the loading condition is as shown in Figure 3 and Table 1.

Fig. 3

Figure 4 gives the change in mechanical property of the human femur during the remodelling and healing process. This phenomenon can be explained as follows. When the resorption is in progress for the first 20 days, the sensor cells can detect the damage and send out corresponding stimulation signals. Beyond the critical value, the damage remodelling will be dominant. In addition, after the absorption activity occurs, it often diffuses in the nearby area. The bone density in the model can be reduced to the minimum. Then a reversal period follows from day 20 to 50. It is assumed that the absorption and formation rates are the same. Finally, the formation stage starts from the 50th to the 100th day. With the deposition of new bone, the mechanical properties and microstructure of bone gradually returned to the original state. As is shown in Figure 3, our simulation result is very close to the real bone. The distribution of bone mass in the bone is similar to that of real ones.

Fig. 4

In order to describe the damage rate represented by elastic modulus, damage D_j is introduced into the model, which is expressed as follows: (8) Dj=1−Ej/E0 {D_j} = 1 - {{{E_j}} \mathord{\left/{\vphantom {{{E_j}} {{E_0}}}} \right.} {{E_0}}} where E_j is the elastic modulus of an element in step j and is the modulus of the element in the undamaged model. The bone damage evolution can be seen in Figure 5. It is shown that when remodelling reduces the elastic modulus of bones, the damage signal gradually diffuses to other areas near this position. We can draw the conclusion that a series of bone remodelling and repair behaviours are carried out not only in the damaged area but also in the whole femoral model. This influence relationship is mainly determined by the spatial influence function, so that when the signal is strong enough, the signal of damage element can be transmitted to the whole model.

Fig. 5

In this case, we will change the loadings to investigate the influence of those loadings on the bone remodelling process as shown in Tables 2 and 3 and Figure 6.

Fig. 6

The results are shown in Figure 7. As is shown, the implant was tightly in contact with the bone when it was inserted into the tissue. But due to the stress shield effect, the loadings on the inner surface of the bone are not enough to support the normal condition of the bone. Then the inner bone began to absorb and led to a risk of loosening and failure of the structure. This is because the elastic modulus of the pin is much bigger than that of the bone. Most loadings are carried by the implant. The bone will fall into a condition called disuse. This is the main reason which causes the failure of a bone implant.

Fig. 7

However, this condition can be improved by changing its loading condition. We can see in Figure 7(a) that the absorption mainly occurred on the left side. Then we adjust the forces on the right side. The results show that the absorption resulting from the disuse is less than that in the first case.

Under the condition of ensuring the support strength of a prosthesis, choosing a material with a smaller elastic modulus can also reduce the bone stress shielding effect. Scanell et al. [63] considered that the higher the elastic modulus, the greater the stress shielding and the more the bone absorption. Sumner et al. [64] found that the material stiffness has great influence on stress shielding. In order to study the influence of stiffness on the service life of a prosthesis, three kinds of prostheses with different elastic modulus are selected in this paper. The specific parameters are E = 200, 100 and 20 GPa. In order to study the stress shielding, the load is still set according to the load distribution before adjustment. After 500 steps of simulation, the results are as shown in Figures 4–20. When the stiffness of the femoral prosthesis is in the range of 20a–200a, we can see that the stiffness of the femoral prosthesis increases slowly when the stiffness of the prosthesis is less than 20A. With a small axial load, the rate of bone formation in the femur is greater than that of bone resorption, and the combination of femur and prosthesis is better. When the prosthesis was 100 GPa, the number of bone units was very close to that of 20 GPa. This indicates that the 100 and 20 GPa prosthesis have almost the same effect in reducing stress shielding. This is consistent with the results of the influence of prosthesis stiffness on bone reconstruction analysed by Qu Chuanyong [10]. In order to ensure that the prosthesis has enough stiffness to provide support for patients, we will use a more moderate 100-GPa prosthesis for analysis of the disuse load of different ages.

Fig. 8

Fig. 9