Research Article  Open Access
Kaysar Rahman, Mamtimin Geni, Mamatjan Mamut, Nijat Yusup, Muhtar Yusup, "Experimental Observation of the Skeletal Adaptive Repair Mechanism and Bionic Topology Optimization Method", Mathematical Problems in Engineering, vol. 2014, Article ID 593867, 12 pages, 2014. https://doi.org/10.1155/2014/593867
Experimental Observation of the Skeletal Adaptive Repair Mechanism and Bionic Topology Optimization Method
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 selfrepairing process in bone defects in rabbit femur. Experimental result shows that during selfrepairing 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 crosstype reactiondiffusion 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 selfrepairing 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 selfrepairing process of bone defect and real crosstype 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 selfrepairing process of bone defect and real crosstype relations between the normalized bone mass volumes of enamel bone and spongy bone are examined. According to the crosstype feature, a new bone remodeling model based on a nonlinear crosstype reactiondiffusion 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 stresssensitive reactiondiffusion equation. The stresssensitivity 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, selfhealing process in bone defects is observed experimentally. In Section 3, new bone remodeling model based on crosstype reactiondiffusion 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 SelfHealing 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 34 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 selfrepair 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 followup measurement benchmark. Repeat steps and to measure the bone metadata of selfhealing process.
(a)
(b)
(c)
(d)
Punch position of the thighbone and result of the selfrepairing of experimental rabbit bone defect is shown in Figure 2. The opening diameter of rabbit thighbone is approximately 4 mm. Evolutionary selfrepairing 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 selfhealing 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 selfrepair process of bone defect are depicted in Figure 4. From Figure 4 it can be obviously seen that during selfrepairing 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 CrossType ReactionDiffusion Equations
Reactiondiffusion 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 socalled 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 crosstype feature showed in Figure 4, a hypothetical model for bone remodeling based on the crosstype reactiondiffusion 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 crosstype phenomena. After a certain time period both factors form a static distribution.
(a) Tezuka’s model [7]
(b) Present crosstype reactiondiffusion model
Kondo and Asai [11] use the following model based on the Turing reactiondiffusion 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 crosstype reactiondiffusion 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 crosstype 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 crosstype reactiondiffusion 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 zeroflux boundary conditions and spatial domain is discretized with and in the reactiondiffusion 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 activatorinhibitor (Turing reactiondiffusion) system and must be in phase and peaks should be at the same place as peaks. However, in the crosstype reactiondiffusion 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 crosstype reactiondiffusion 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 reactiondiffusion 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 crosstype 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 reactiondiffusion 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 selfrepairing 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 crosstype reactiondiffusion equations (5) are shown in Table 2.
Example 1 (bone selfrepairing process). In order to check the bone selfrepairing process, a twodimensional 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 righthand 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 wellknown 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 meshindependence 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 crosstype reactiondiffusion 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 crosstype phenomena in the process of selfhealing bone volume change between the enamel bone and spongy bone. It indicates that the bone volume changes in the relationship between selfrepairing processes of bone defect are caused by the mechanism of bone remodeling process.
According to the real bone crosstype remodeling characteristics in experiments, modified equation (4) by addition of the coefficients , and stress term , bone remodeling model based on a crosstype reactiondiffusion system is proposed. Then the appropriate coefficients , , are defined by tuning the crosstype equation (5) and its stability.
New bionic topology optimization method by combining crosstype reactiondiffusion 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 selfadjoint 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).
References
 G. H. Meyer, “Die architektur der spongiosa, archief fur den anatomischen und physiologischen,” Wissenschaften im Medicin, vol. 27, no. 4, pp. 1389–1394, 1867. View at: Google Scholar
 J. Wolff, The Law of Bone Remodeling, Springer, Berlin, Germany, 1986.
 H. M. Frost, “Bone “mass” and the “mechanostat”: a proposal,” Anatomical Record, vol. 219, no. 1, pp. 1–9, 1987. View at: Publisher Site  Google Scholar
 H. M. Frost, “The Utah paradigm of skeletal physiology: an overview of its insights for bone, cartilage and collagenous tissue organs,” Journal of Bone and Mineral Metabolism, vol. 18, no. 6, pp. 305–316, 2000. View at: Publisher Site  Google Scholar
 L. M. McNamara and P. J. Prendergast, “Bone remodelling algorithms incorporating both strain and microdamage stimuli,” Journal of Biomechanics, vol. 40, no. 6, pp. 1381–1391, 2007. View at: Publisher Site  Google Scholar
 R. Huiskes, R. Rulmerman, G. H. van Lenthe, and J. D. Janssen, “Effects of mechanical forces on maintenance and adaptation of form in trabecular bone,” Nature, vol. 405, no. 6787, pp. 704–706, 2000. View at: Publisher Site  Google Scholar
 K.I. Tezuka, Y. Wada, and M. Kikuchi, “iBone: a reaction diffusion based shape optimization method,” Key Engineering Materials, vol. 243244, pp. 601–606, 2003. View at: Publisher Site  Google Scholar
 K.I. Tezuka, Y. Wada, A. Takahashi, and M. Kikuchi, “Computersimulated bone architecture in a simple boneremodeling model based on a reactiondiffusion system,” Journal of Bone and Mineral Metabolism, vol. 23, no. 1, pp. 1–7, 2005. View at: Publisher Site  Google Scholar
 S. V. Komarova, R. J. Smith, S. J. Dixon, S. M. Sims, and L. M. Wahl, “Mathematical model predicts a critical role for osteoclast autocrine regulation in the control of bone remodeling,” Bone, vol. 33, no. 2, pp. 206–215, 2003. View at: Publisher Site  Google Scholar
 K.I. Tezuka, M. Yasuda, N. Watanabe et al., “Stimulation of osteoblastic cell differentiation by Notch,” Journal of Bone and Mineral Research, vol. 17, no. 2, pp. 231–239, 2002. View at: Publisher Site  Google Scholar
 S. Kondo and R. Asai, “A reactiondiffusion wave on the skin of the marine angelfish Pomacanthus,” Nature, vol. 376, no. 6543, pp. 765–768, 1995. View at: Publisher Site  Google Scholar
 L. Jia, M. Geni, H. Eli, X. Abdikerem, and M. Kikuchi, “Bone formation based on the turing model under compressed loading condition,” Advanced Materials Research, vol. 33–37, pp. 1011–1016, 2008. View at: Publisher Site  Google Scholar
 M. Geni and M. Kikuchi, “Shape optimization of metal welded bellows seal based on the turing reactiondiffusion model coupled with FEM,” Key Engineering Materials, vol. 385–387, pp. 813–816, 2008. View at: Publisher Site  Google Scholar
 J. Mamatrozi, M. Geni, A. Mahmut, M. Mamut, Z. Rui, and A. Hamudun, “Experimental study on embryo skeletal development to postnatal growth process of rabbit thighbone,” Key Engineering Materials, vol. 462463, pp. 1158–1163, 2011. View at: Publisher Site  Google Scholar
 A. Mahemuti, M. Geni, J. Maitirouzi, M. Mamuti, R. Zhang, and Maimaitituxun, “Topology shape and mechanism of variation of bone layer during growing process of femur,” Journal of Shanghai Jiaotong University: Medical Science, vol. 32, no. 3, pp. 263–268, 2012 (Chinese). View at: Publisher Site  Google Scholar
 M. Mamatjan, G. Mamtimin, Y. Nijat et al., “Analysis changing the rule of sclerotin volume duringthe selfrepairing process of bone defect,” Journal of Biomedical Engineering, vol. 29, no. 4, pp. 682–686, 2012 (Chinese). View at: Google Scholar
 J.M. Rossi and S. WendlingMansuy, “A topology optimization based model of bone adaptation,” Computer Methods in Biomechanics and Biomedical Engineering, vol. 10, no. 6, pp. 419–427, 2007. View at: Publisher Site  Google Scholar
 M. Bagge, “A model of bone adaptation as an optimization process,” Journal of Biomechanics, vol. 33, no. 11, pp. 1349–1357, 2000. View at: Publisher Site  Google Scholar
 M. G. Mullender, “A physiological approach to the simulation of bone remodeling as a selforganizational control process,” Journal of Biomechanics, vol. 27, no. 11, pp. 1389–1394, 1994. View at: Publisher Site  Google Scholar
 P. Fernandes, H. Rodrigues, and C. Jacobs, “A model of bone adaptation using a global optimisation criterion based on the trajectorial theory of wolff,” Computer Methods in Biomechanics and Biomedical Engineering, vol. 2, no. 2, pp. 125–138, 1999. View at: Publisher Site  Google Scholar
 K. Cai, B. S. Chen, and H. W. Zhang, “Topology optimization of continuum structures based on a new bionics method,” International Journal of Computational Methods in Engineering Science and Mechanics, vol. 8, no. 4, pp. 233–242, 2007. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 M. P. Bendsøe and O. Sigmund, Topology Optimization: Theory, Methods, and Applications, Springer, New York, NY, USA, 2003. View at: MathSciNet
 X. Huang and Y. M. Xie, Evolutionary Topology Optimization of Continuum Structures: Methods and Applications, John Wiley & Sons, Chichester, UK, 2010.
 S. Shojaee and M. Mohammadian, “Structural topology optimization using an enhanced level set method,” Scientia Iranica, vol. 19, no. 5, pp. 1157–1167, 2012. View at: Publisher Site  Google Scholar
 A. M. Turing, “The chemical basis of morphogenesis,” Bulletin of Mathematical Biology, vol. 52, no. 12, pp. 153–197, 1990. View at: Publisher Site  Google Scholar
 T. Miura and P. K. Maini, “Periodic pattern formation in reactiondiffusion systems: an introduction for numerical simulation,” Anatomical Science International, vol. 79, no. 3, pp. 112–123, 2004. View at: Publisher Site  Google Scholar
 M. Geni, W. Xufei, and M. Kikuchi, “Study on selfconsistent mesh generating method of hexahedron element based on the local waveform method with damping,” Key Engineering Materials, vol. 306–308, pp. 607–612, 2006. View at: Publisher Site  Google Scholar
 A. Kaveh, B. Hassani, S. Shojaee, and S. M. Tavakkoli, “Structural topology optimization using ant colony methodology,” Engineering Structures, vol. 30, no. 9, pp. 2559–2565, 2008. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2014 Kaysar Rahman et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.