Abstract
Bone adaptive repair theory considers that the external load is the direct source of bone remodeling; bone achieves its maintenance by remodeling some microscopic damages due to external load during the process. This paper firstly observes CT data from the whole self-repairing process in bone defects in rabbit femur. Experimental result shows that during self-repairing process there exists an interaction relationship between spongy bone and enamel bone volume changes of bone defect, that is when volume of spongy bone increases, enamel bone decreases, and when volume of spongy bone decreases, enamel bone increases. Secondly according to this feature a bone remodeling model based on cross-type reaction-diffusion system influenced by mechanical stress is proposed. Finally, this model coupled with finite element method by using the element adding and removing process is used to simulate the self-repairing process and engineering optimization problems by considering the idea of bionic topology optimization.
1. Introduction
It is generally agreed that osteoblasts and osteoclasts are mainly involved in trabecular bone remodeling process. Bone is continuously remodeled by bone formation cells (osteoblasts) and resorption cells (osteoclasts) to keep balance between the newly formed bone mass and the missing bone. In this way, bone is able to optimize its biological structure to adapt itself to the natural environments. Based on the work by Meyer [1], Wolff [2] proposed a theory that the trabecular architecture of bone matches the stress trajectories. Frost [3, 4] proposed a feedback mechanism—the “mechanostat” that controls bone mass in response to mechanical loading. The influence of the mechanical stress on the behavior of osteoclasts and osteoblasts is not yet fully understood. However, it seems that mechanosensitive osteocytes are capable of transuding changes in mechanical stimuli into biochemical signals which regulate cellular responses [5]. In recent years, a number of models have been developed to simulate and explain this biological regulatory process at the cellular and microstructure level [6–9]. Huiskes et al. introduced “bone cell plays a mechanical sensor role and it can transform the stress information to osteoblast and osteoclast” [6]. Tezuka et al. [10] found that human bone marrow mesenchyme stem cells formed periodic cell condensates and are similar to Turing patterns when they are forced to differentiate into osteoblasts by the action of notch. Then they constructed a bone remodeling model based on Turing system with linear kinetic reaction term and named this model as “iBone” (imitation Bone) [7, 8]. This Turing system is used as well by Kondo and Asai [11] to simulate skin patterns of tropical angelfish. Jia et al. [12] generalized “iBone” model and applied this model to bone formation and resorption process under different boundary conditions; Geni and Kikuchi [13] applied this “iBone” model to shape optimization of Metal Welded Bellows Seal problems. Mamatrozi et al. [14] and Mahemuti et al. [15] observe skeletal morphology, the distributions of bone density, and volume of the thighbone by using the CT data from the experiment. Mamatjan et al. [16] propose self-repairing process of bone defect and real cross-type relations between the normalized bone mass volumes of enamel bone and spongy bone are observed.
There is a close connection between the objective of computational bone adaptation to mechanical loading and structural topology optimization. In this field, there are two hot topics. The first one considers bone adaptation as an optimization process [17–19] and the other one considers the rule of cancellers bone adaptation as a basic tool for the optimal design of continuum structures [20, 21].
Topology optimization is often regarded as layout optimization or generalized shape optimization and works as a powerful tool assisting designers in the selection of suitable initial structural topology from an optimal point of view. Structural topology optimization has made remarkable progress over the past three decades. It has found its way into industry and is used in a variety of engineering fields. Considerable research and various topology optimization methods, such as homogenization [22], solid isotropic material with penalization (SIMP) [22], evolutionary structural optimization (ESO and BESO) [23], and level set (LS) method, have been proposed [24].
In this study, to understand how bone cells can create a structure adapted to the mechanical environment, the bone remodeling process is discussed experimentally by using CT scan of New Zealand white rabbit thighbone. Through the observation of the growing process, the self-repairing process of bone defect and real cross-type relations between the normalized bone mass volumes of enamel bone and spongy bone are examined. According to the cross-type feature, a new bone remodeling model based on a nonlinear cross-type reaction-diffusion equation influenced by mechanical stress is proposed. This model is coupled with finite element method (FEM) to establish new bionic topology optimization model based on the idea of Huiskes and Tezuka’s work. The optimization of bone distribution can be achieved by adding or removing elements which is calculated by stress-sensitive reaction-diffusion equation. The stress-sensitivity distribution is calculated by using FEM. This approach provides a novel means of understanding of the remodeling processes involved in fracture repair and the treatment of bone diseases.
The organization of the paper is as follows. In Section 2, self-healing process in bone defects is observed experimentally. In Section 3, new bone remodeling model based on cross-type reaction-diffusion equations is proposed. In Section 4, a new bionic topology optimization method is introduced. In Section 5, by using the new method the numerical simulation of bone defect repair processes and topology optimization for continuous structural problem are conducted. Conclusion of the presented paper is given in Section 6.
2. Experimental Observation of Self-Healing Process of Bone Defects
In this study, four healthy New Zealand white rabbits (Xinjiang Animal Test Center provided) were used. There is no gender preference in selecting bone. The rabbit body mass is around ~ when ages are around 4 months; the rabbits were anesthetized with 3% pentobarbital sodium (30 mg/kg) by injecting in the ear vein. Then the anesthetized rabbits were fixed and the legs were shaved and disinfected. The femur anterolateral incision was made and the femur was exposed. Under sterile conditions, a small hole is drilled in the femur anterolateral in the location of upper 1/4. The shape of the hole mouth is a circular shape with diameter of 3-4 mm; the skin is sutured at the end of the operation. Regular CT scan observations of bone defect area were held from the next day. In this experiment, MIMICS (Materialists Interactive Medical Image Control System) software was used and the CT image data can be easily visualized with it; bone element filling method is proposed to observe the bone defect self-repair process. The volume of bone is calculated by using the following equation: where is the volume of bone, is the number of pixel element, is the pixel size, and is the slice thickness. The bone volume in the hole is calculated by using following steps shown in Figures 1(a), 1(b), 1(c), and 1(d), respectively. (1)Start the MIMICS software and input the initial opening CT data.(2)Find the initial opening position of the femoral bone and insert it into a relatively small diameter cylinder (Figure 1(a)), then adjust the cylinder to the hole in the middle (shown in Figure 1(b)), finally, to contain the defective parts of the drilled bone holes, the cylinder is slightly enlarged so that cylinder just matches the hole (Figure 1(c)), and fix this cylinder as the following measurement benchmark.(3)Input the subsequently scanned CT data.(4)Inset benchmark cylinder, and then remove the outer periphery of the excess portion of benchmark cylinder, retaining only inside part of the bone (see Figure 1(d)). Record the number of the bone element, pixel size, and slice thickness for calculating the total bone volume in the hole.(5)Keep the initial cylinder size as a follow-up measurement benchmark. Repeat steps and to measure the bone metadata of self-healing process.
(a)
(b)
(c)
(d)
Punch position of the thighbone and result of the self-repairing of experimental rabbit bone defect is shown in Figure 2. The opening diameter of rabbit thighbone is approximately 4 mm. Evolutionary self-repairing histories of the experimental rabbit bone defect versus day are shown in Figures 3(a)–3(f), respectively. From Figures 2 and 3 it can be seen that self-healing rabbit bone defects of the hole became repaired after 69 days and completely repaired after 159 days. Dimensionless process is held in order to better observe the relationship between the volume changes of the respective bone layer. The volume changes between the spongy bone and enamel bone during the growth days in self-repair process of bone defect are depicted in Figure 4. From Figure 4 it can be obviously seen that during self-repairing process there is an interaction relationship that exists between spongy bone and enamel bone volume changes of bone defect, that is, when volume of spongy bone increases and enamel bone decreases otherwise vice versa.
(a) 1 day
(b) 4 days
(c) 12 days
(d) 25 days
(e) 69 days
(f) 159 days
3. Bone Remodeling Model Based on Cross-Type Reaction-Diffusion Equations
Reaction-diffusion model was first proposed by Turing [25] and it has been widely used to simulate biological pattern formation [11, 26]. The basic Turing model consists of two hypothetical molecules, activator and inhibitor, which interacts with each other and diffuses independently. These molecules produce a system with positive and negative feedbacks and generate stable periodical patterns so-called Turing patterns. Tezuka et al. constructed a bone remodeling model based on Turing system with linear kinetic reaction term and named this model as “iBone” [7, 8]. A hypothetical model of Tezuka is shown in Figure 5(a). When local stress is applied, the local concentration of activator increases first, followed by an increase in that of the inhibitor. In this study, to understand how bone cells can create a structure adapted to the mechanical environment, and according to cross-type feature showed in Figure 4, a hypothetical model for bone remodeling based on the cross-type reaction-diffusion equations is established. The present hypothetical model is showed in Figure 5(b); when local stress is applied, the local concentration of activator (solid line) increases and decreases in that of the inhibitor (dot line) simultaneously, and the relations between the activator and the inhibitor establish a cross-type phenomena. After a certain time period both factors form a static distribution.
![(a) Tezuka’s model [7]](https://static.hindawi.com/articles/mpe/volume-2014/593867/figures/593867.fig.005a.jpg)
(a) Tezuka’s model [7]
(b) Present cross-type reaction-diffusion model
Kondo and Asai [11] use the following model based on the Turing reaction-diffusion model to simulate skin patterns of tropical angelfish: where and denote the time and space dependent concentrations of activator and inhibitor, and both substances yield diffusion reaction with each other and decay, the is a Laplacian operator, and the and are the diffusion constants for the activator and inhibitor, respectively, which describes the speed of diffusion, and are decay constants, and are compensation coefficient for the activator and inhibitor, respectively, and the set of the , and are the parameters for the reaction terms for the activator and inhibitor, respectively.
Tezuka et al. [7] modify Kondo’s model by adding stress term and propose a hypothetical model of bone remodeling as shown in (3) to simulate bone architecture. There are some works that used the Tezuka model and give a good agreement in the structural optimization problems [13]:
Miura and Maini [26] use the cross-type reaction-diffusion model (4) to simulate periodic pattern formation and give a good agreement for pattern formation process:
In this study, in order to account for real bone cross-type remodeling characteristics from experiments as shown in Figure 4, (4) is modified by adding the coefficients and stress term , and then bone remodeling model based on cross-type reaction-diffusion equations is proposed as shown in the following equation: where and are decay constants for the activator and inhibitor, respectively, the set of the are similar to (2) and (3) and denote reaction terms compensation parameters for the activator and inhibitor, respectively, and the term is the local Von Misses stress at each element which is calculated by using FEM.
The conditions for the emergence of patterns in this model system are that the inhibitor diffuses more rapidly than the activator; that is, , and inhibitor adapts rapidly to changes of the activator in case that it decays more rapidly than the activator; then it should be . Each parameter is set so that concentration of both activator and inhibitor are evenly distributed when no mechanical loads are applied. Suppose that spatial domain size is with zero-flux boundary conditions and spatial domain is discretized with and in the reaction-diffusion equations (3) and (5). In that case, we have 20 pieces of tissue which contain both activator and inhibitor with random numbers. This is described in Microsoft Excel spreadsheet by using Euler forward difference scheme. Abscissa represents the concentration of activator or inhibitor in 20 pieces of tissue . Ordinate means distribution of activator and inhibitor at a later time of Tezuka model (3) and present model (5) and their 30th iteration step with no mechanical loads are shown in Figures 6(a) and 6(b), respectively. Typical parameters are used for the calculation of Tezuka model in (3) as shown in Table 1, and typical parameters for the model in (5) are shown in Table 2.
(a) Tezuka model (3)
(b) Present model (5)
The solid line with square symbol denotes the distribution of (activator) and the dash line with circle symbol denotes the distribution of (inhibitor). In the original activator-inhibitor (Turing reaction-diffusion) system and must be in phase and peaks should be at the same place as peaks. However, in the cross-type reaction-diffusion equations (5) the peaks are at valleys. This is consistent with experimental result in Figure 4 and this may be very interesting and very important research topic.
4. Bionic Topology Optimization Method
In this paper a new bionic topology optimization method based on cross-type reaction-diffusion equation is presented. The major idea of the present approach is to consider the structure to be optimized as a piece of bone which obeys Wolff’s law, and the process of finding the optimum topology of a structure is equivalent to the bone remodeling process. For the sake of simplicity of the model, we delineate the complex bone remodeling system via the actions of two hypothetical molecules; furthermore, the local bone formation and resorption activity are calculated, by the concentrations of activator and inhibitor at each position. In this hypothetical model, bone formation will occur when the ration of activator to inhibitor is over the threshold, and the bone resorption will occur when ratio is over the threshold in the circumferential region near the stress center. Interactions among bone cells are simulated by conventional finite element analysis (FEA) by using the pixel element based on the reaction-diffusion equation (5). This equation is used to calculate the changes in concentrations of these molecules such as activator and inhibitor. For the numerical simulation of bone remodeling process such as changing of the shape of the numerical models by considering the bone formation and resorption phenomenon, the element removing and adding method is used. The displacement value for element is calculated by considering the combined effects of local concentration of activator and inhibitor and expressed as . The parameter “” denotes the threshold of the activator/inhibitor . It is noted that “” the single global parameter affects the balance between bone formation and bone resorption. Shapes of the models change via addition or removal pixel of elements. If , pixel elements will be removed, and if , new pixel elements will be added. Details are shown in Figure 7(b).
(a) Flowchart of bionic topology optimization
(b) Element adding and removing process under the cross-type phenomena
Figure 7(a) shows the flowchart of bionic optimization process, outlined as follows.(1)Create initial design domain and discretize the design domain using a finite element mesh.(2)Give boundary condition and material property and carry out finite element analysis for the structure.(3)Give parameters of reaction-diffusion equation (5) and compute and ,then decide adding or removing elements, and then generate new topology of model.(4)If topology of model is convergent, generate final topology model; otherwise, repeat steps - until the constraint volume (V) is achieved and the total strain energy of the structure reaches a prescribed limit.
5. Numerical Examples and Results
In this section, two numerical examples are presented to demonstrate the effectiveness of the proposed method. In all examples the objective is to minimize the total strain energy with a different material volume constraint. Numerical simulation of self-repairing process of bone defect model and cantilever beam problem, one of the widely used examples in structural topology optimization, are carried out by using presented bionic topology optimization method. No additional treatments (such as sensitivity filtering) for numerical instability control are applied during the optimization. The models are created and the results are visualized by FAST_wmg [27]. In all examples the parameters used for cross-type reaction-diffusion equations (5) are shown in Table 2.
Example 1 (bone self-repairing process). In order to check the bone self-repairing process, a two-dimensional model (10 mm 10 mm and 1 mm = 10 pixels) with single hole (diameter of the hole 4 mm) similar to experimental rabbit defect model was made and used. The initial punch position of rabbit femur and numerical design domain under biaxial uniform tensile force () are shown in Figures 8(a) and 8(b), respectively. In this example the bone is assumed to be an isotropic material and its property of 5 GPa Young’s modulus and Poisson’s ratio of 0.3 are used. The objective function is to minimize the structural total strain energy while the material usage is limited to 100%. Stress distribution and shape change of the hole after each cycle of calculation is observed and results at different iteration steps are visualized and shown in Figures 9(a)–9(f), respectively. From Figure 9, it can be seen that after 30 steps of calculation the shape of the hole becomes repaired and even stress distribution is formed and results are close to the experimental one. The changes of the strain energy and the volume constraint over the iterations are depicted in Figures 10(a) and 10(b), respectively. The strain energy is optimized from 2.5863 to 1.0000 and a uniform distribution of the strain energy along structural design boundary is achieved after 30 iterations. From this figure, we can see that the strain energy converges in a fast and stable way due to the present bionic optimization method.
(a) Step 1
(b) Step 5
(c) Step 10
(d) Step 20
(e) Step 25
(f) Step 30
(a) Strain energy versus iteration
(b) Volume fraction versus iteration
Example 2 (cantilever beam problem). In order to check if the bionic optimization method is applicable for engineering optimization problems, cantilever beam problem, one of the widely used examples of structural topology optimization, is solved. Figures 11(a) and 11(b) show the design domain and initialized topology with uniformly distributed holes of a cantilever beam. The boundary of the left side is fixed and a vertical concentrated force is loaded in the center of the right-hand side. The size of the design domain is and finite element mesh model is made by the size pixel element. The material volume constraint is set to be 45% of the whole design domain. The material properties as Young’s modulus and Poisson’s ratio are used, respectively. The final topologies by using present method after 50 iterations are shown in Figure 11(c). The layout is also optimized by using ant colony optimization (ACO) [28], BESO (bidirectional evolutionary structural optimization) with [23], and SIMP (solid isotropic material with penalization) [22] methods which are shown in Figures 12(a), 12(b), and 12(c), respectively. With the comparison of these approaches, it can be seen that the final design of this study is similar to the well-known examples of the known literature [22–24, 28]. Topology optimization results obtained with different mesh resolutions (, , ) are also presented in Figures 13(a), 13(b), and 13(c), respectively. It is found that the final topologies in all three cases are very similar. These results indicate that there might be less mesh-independence property of the proposed method. Variation of the strain energy and volume fraction versus iteration steps throughout the optimization process is presented in Figures 14(a) and 14(b), respectively. It also indicates that the strain energy increases initially as the volume fraction decreases and then it is convergent to an almost constant value after the objective volume is achieved. In this case, the objective function converges to the final value of 1.52 and the volume fraction ratio remains 0.45 after 50 iteration steps.
(a) Initials design domain
(b) Initialized topology with uniformly distributed holes
(c) Final topologies
(a) Ant colony method (ACO)
(b) BESO method
(c) SIMP method
| (a) |
| (b) |
| (c) |
(a) Strain energy versus iteration
(b) Volume fraction versus iteration
6. Conclusion
A new bionic topology optimization method based on the cross-type reaction-diffusion equation is proposed in this paper. In this method bone adaptation process is the decision for an element to be eliminated or to be added. On the one hand, the bone remodeling phenomenon is a biological realization of the optimal structure principle, requiring equal value of the energy distribution on the structure surface. On the other hand, the optimization scenario is based on the osteoblasts and osteoclasts activity. Several conclusions are summarized as follows.
Experimental results show that there exist cross-type phenomena in the process of self-healing bone volume change between the enamel bone and spongy bone. It indicates that the bone volume changes in the relationship between self-repairing processes of bone defect are caused by the mechanism of bone remodeling process.
According to the real bone cross-type remodeling characteristics in experiments, modified equation (4) by addition of the coefficients , and stress term , bone remodeling model based on a cross-type reaction-diffusion system is proposed. Then the appropriate coefficients , , are defined by tuning the cross-type equation (5) and its stability.
New bionic topology optimization method by combining cross-type reaction-diffusion equations and describing bone adaptation process is implemented by using finite element analysis. The effectiveness of the proposed bionic topology optimization method is verified by some numerical examples of simple bone defects medical problems and common structure. The result shows that the bionic optimal configurations are close to that of the experimental one and the conventional structural optimization results which are obtained from other mature topology optimization approaches.
However, the present bionic topology optimization method can be applied to the self-adjoint optimization problems with simple volume constraints. Further work for more complex optimization problems is needed, and we will pursue in future.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The authors wish to acknowledge the support of the China Natural Science Foundation (no. 50775193, no. 81160542, and no. 11261057) and the National Key Basic Research and Development Program (973 Program no. 2011CB706600).