Jtam.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 48, 4, pp. 855-870, Warsaw 2010 MULTISCALE MODELING OF OSSEOUS TISSUES Tadeusz Burczyński Silesian University of Technology, Department of Strength of Materials and Computational Mechanics, Gliwice, Poland; Cracow University of Technology, Institute of Computer Science, Kraków, Poland e-mail: tadeusz.burczynski@polsl.pl Wacław Kuś Anna Brodacka Silesian University of Technology, Department of Strength of Materials and Computational Mechanics, Gliwice, Poland e-mail: waclaw.kus@polsl.pl; anna.brodacka@polsl.pl The paper presents a methodology of the multiscale bone modeling in which the task of identification of material parameters plays the crucial role. A two-scale analysis of the bone is considered and the problem of identification, formulated as an inverse problem, is examined as an important stage of the modelling process. The human femur bone, built form cancellous and cortical bone, is taken as an example of an osseous tissue, and the computational multiscale approach is considered. The methodology presented in the paper allows one to analyze the two-scale modelwith theuseof computationalhomogenization.The representative volume element (RVE) is created for the microstructure of the basis of micro-CT scans. Themacro andmicromodel analyses are performed by using thefinite elementmethod.The identificationof trabeculaematerial parameters on themicro-level is considered as theminimization problem which is solved using evolutionary computing. Keywords:multiscalemodeling of bone, computational homogenization, identification of material parameters 1. Introduction The bone is a kind of connective tissue which is the basic constituent of the skeleton. The osseous tissue consists of cells and the intercellular substance. Thebone is a relatively hardand lightweight compositematerial. These excep- tional mechanical properties of bone aremainly due to its specific hierarchical 856 T. Burczyński et al. microstructure. From the macroscale down to the nanometer one, it is po- ssible to distinguish: cancellous and cortical tissues, Haversian systems and osteons, lamellae, collagen fibres, collagen fibrils, and elementary constituents (collagen, mineral, water, etc.). The modeling of biological tissues has long history. Almost always, the traditionalmodels are focusedon singlebiological scales.But in order tomodel a bone, it is important to consider the hierarchical structure and to study how the properties at each scale are governed by the material organization at the lower levels (Sansalone et al., 2009). The main intention of the modelling of cancellous bone is to decrease the extent of necessary laboratory tests. Applications of multiscale approaches, based e.g. on the finite element method (FEM), enable modeling of the re- presentative volume element (RVE) and calculation of the effective material parameters, including analysis of their change with respect to increasing po- rosity (Ilic et al., 2010). The hierarchical multiscale modeling method for the analysis of the cortical bone consists of two boundary value problems, one for the macroscale and another for the microscale. The coupling between these scales is done using the computational homogenization scheme. In this appro- ach, virtually there is no limitation on the geometry andmaterialmodel of the RVE and constituents, so it can be applied to any kind of nanostructuredma- terial. In addition to thedetermination of the global behavior andproperties of bone, themicrostructural results such as stress distribution in the constituents are also available. By considering various RVE’s regarding themicrostructure of bone in various locations of the macroscopic piece of bone, it is possible to simulate the real behavior of the bone (Ghanbari andNaghdabadi, 2009). The multiscale analysis can be also performed with use of the bridging techniques described in Burczyński et al. (2010) on the levels of nano andmicrometers. The multiscale analysis of remodeling with the use of an artificial micro- structure of the cancellous bone can be found in Kowalczyk (2010). Themul- tiscale approach to the modeling of the cortical bone is considered in Hamed et al. (2010). The paper byAgić et al. (2006) describes analytical calculations of cancellous bonematerial coefficients on the basis of experiments. The goal of the present paper is to present amethodology of themultiscale bonemodeling inwhich the task of identification ofmaterial parameters plays the crucial role. The two scale analysis of the bone is considered and the pro- blem of identification, formulated as an inverse problem, is considered as an important stage of the modeling process. The human femur bone is taken as an example of the osseous tissue and the computational multiscale approach is considered. The femur bone is build form cancellous and cortical bone. The Multiscale modeling of osseous tissues 857 cortical bone is very stiff.The cancellous (trabecular) bone is aporous structu- re built from trabeculae. The head of the femur contains cancellous bonewith high density. The methodology presented in the paper allows one to analyze the two-scale model taking into account the cancellous bone structure. The representative volume element geometry (RVE) is created for microstructure on the basis of micro-CT scans. The computational homogenization method (Terada andKikuchi, 2001;Kouznetsova, 2002) is used to obtain averagedma- terial properties of the microstructure. The macro and micro model analyses are performed by using the finite element method (FEM) (Zienkiewicz et al., 2005). The identification of trabeculaematerial parameters on themicro-level is considered as the minimization problem which is solved using evolutionary computing. 2. The multiscale structure of proximal femur bone Theproximal femurbone is shown inFig.1.The cancellous andcortical tissues are the two main components of the bone. The cancellous tissue is a porous structure with complicated geometry (Fig.1). Fig. 1. The proximal femur bone geometry andmicrostructure The geometry of the cancellous bone changes in different locations of the femur bone. The experimental tests of single trabeculae revealed isotropic material behavior (Tsubota et al., 2003), however the material properties of the cancellous tissue are no longer isotropic. The most common model of the material used for the trabecular bone modeling is transversely isotropic and orthotropic. The cortical bone has also hierarchical structure as shown in Hamed et al. (2010). 858 T. Burczyński et al. The material properties depend on the structure of tissues and change depending on the location in the femur bone. They can be obtained on the basis of density obtained fromCT scans in many papers (Writz et al., 2000). The one-scale analysis takes into account average material properties. The computational homogenization allows one to performan analysis in which the material properties depend on the microstructure of the trabecular bone. 3. Computational homogenization of osseous tissue Themultiscale modeling allows one to take into account the dependences be- tween two or more scales (Fig.2). One of the methods used in the multiscale modeling is computational homogenization (Madej et al., 2008). Aheterogene- ousmaterial is replaced with a homogenous one (Fig.3). The homogenization is useful when the microstructure is periodic. The influence between scales in the computational homogenization is obtained on the basis of a numerical solution to the boundary value problem performed in each scale. Fig. 2. Mulsticale modeling The two-scale analysis of the bone system is shown inFig.4. Theboundary value problem forRVE should be solved for eachGauss integration point. The strain values are transferred to the micromodel during the localization stage. The traction, displacements or periodicboundaryconditions are applied to the Multiscale modeling of osseous tissues 859 Fig. 3. Homogenization of material, (a) heterogeneous structure, (b) structure after homogenization Fig. 4. Two-scale computational homogenization microstructure. The stresses obtained after boundary value problem analysis are used to obtain homogenized average values which are transferred after the homogenized stage to the higher scale. The relationship between stresses and strains for an orthotropic elastic material are expressed as follows σ=Cε or ε=Sσ (3.1) where σ= [σ11,σ22,σ33,σ12,σ13,σ23] ⊤ (3.2) ε= [ε11,ε22,ε33,ε12,ε13,ε23] ⊤ are vectors of stresses and strains. C and Sare the stiffness and compliancematrices, of the orthotropic linear elastic material, respectively. They can be written as 860 T. Burczyński et al. C=    c11 c12 c13 0 0 0 c22 c23 0 0 0 c33 0 0 0 sym. c44 0 0 c55 0 c66    (3.3) S=C−1 =    s11 s12 s13 0 0 0 s22 s23 0 0 0 s33 0 0 0 sym. s44 0 0 s55 0 s66    where s11 = 1 E1 s22 = 1 E2 s33 = 1 E3 s44 = 1 G12 s55 = 1 G23 s66 = 1 G13 s12 = −ν21 E2 s13 = −ν31 E3 s21 = −ν12 E1 s23 = −ν32 E3 s31 = −ν13 E1 s32 = −ν23 E2 (3.4) where Ei is Young’s modulus along the axis i, Gij is the shear modulus in the direction j on the plane whose normal is in the direction i, and νkl is Poisson’s ratio that corresponds to contraction in the direction j when an extension is applied in the direction k. Due to symmetry of the stiffness and compliance matrices, the 9 variables are independent in the fully orthotropic elasticmaterial (Bąk andBurczyński, 2009; Schneider et al., 2009). The trabecular tissue is often modeled with the use of 5 parameters, assuming E1 =E2 and ν23 = ν31. The material coefficients in the case of linear problems can be obtained once for each microstructure. The six analyses should be performed for each microstructure to obtain the 9 independent orthotropic material coefficients. The average strains and stresses for RVE are defined as follows εavg = 1 |ΩRVE| ∫ ΩRVE ε dΩRVE σavg = 1 |ΩRVE| ∫ ΩRVE σ dΩRVE (3.5) where ΩRVE is the area of RVE. Multiscale modeling of osseous tissues 861 The constitutive relation between them has the form σavg =C h εavg (3.6) where Ch is the stiffness tensor of the equivalent homogenous material that fulfils the elastic deformation characteristic for the heterogeneous material. A detailed description of the algorithm of computational homogenization is presented in Fig.5. Fig. 5. Computational homogenization algorithm 4. Identification One of themost important stages of themodeling process is the identification of geometrical and material parameters of the real system. Problems of mul- tiscale modeling need a special kind of identification procedures (Burczyński and Kuś, 2009). 862 T. Burczyński et al. One should evaluate some material parameters of structures in one sca- le having measured information in another scale. In the considered problem, one should identify Young’s modulus E and Poisson’s ratio ν of the single trabeculae having information aboutmeasuredmaterial parameters cij in the macroscale. The material properties in macroscale can be obtained by per- forming an identification task on the basis of e.g. strains or displacements measured for the macromodel. The material properties can be also obtained by using the ultrasonic velocity measurement or mechanical tests for micro- specimens. The problem can be formulated as aminimization task min E,ν F (4.1) where the objective function is formulated as follows F = n∑ i=1 |ai− âi| (4.2) where ai are computed homogenizedRVEmaterial properties and âi areRVE homogenized material properties from themacromodel. For the homogenized orthotropic material n=9 and ai = {c11,c22,c33,c12,c13,c23,c44,c55,c66} The objective function is multimodal in most cases and the optimization should be performed with the use of an algorithm which is resistant to local minima. The wide range of bioinspired algorithms allow one to solve global optimization problems. The minimization problem can be solved using the distributed parallel evolutionary algorithm (Burczyński, 2010). The searched material parameters – Young’s modulus E and Poisson’s ratio ν of the single trabeculae create a chromosome ch= [g1,g2] (4.3) where gi (i = 1,2) are genes: g1 – Young’s modulus E, g2 – Poisson’s ratio ν. Evolutionary algorithms are methods which search through the space of solutions basing on an analogy to the biological evolution of species. They operate on populations of individuals (chromosomes) which can be conside- red as a set of problem solutions. Chromosomes consist of genes which play the role of design variables in a minimization problem. In the proposed evo- lutionary computation, the floating-point representation is applied. It means Multiscale modeling of osseous tissues 863 that genes are real numbers. The fitness function is represented by the objec- tive function F. In a typical sequential evolutionary algorithm, the number of fitness function evaluations during the optimization problem amounts to thousands or more. Because the fitness function evaluation for the considered problem takes a lot of time, a distributed and parallel version of the evolu- tionary algorithm has been proposed. The distributed evolutionary algorithm works similarly to many evolutionary algorithms simultaneously operating on subpopulations, and during themigration phase the chromosomes are exchan- ged between subpopulations. The parallel evolutionary algorithm evaluates fitness function values in a parallel way. More information about such a ge- neralization of evolutionary computing can be found in Burczyński (2010). A flowchart of the distributed parallel evolutionary algorithm is presented in Fig.6. The random starting subpopulations are created on the beginning. The evolutionary operators change gene values similarly as in the biological pheno- mena. The fitness function is calculated for everymodified chromosome in the parallel way. The objective function is computed on the basis of FEManalysis for themicrostructure. The exchange of chromosomes occurs in themigration phase.Thebest chromosomesmigrate between subpopulations.Themigration is performed every a few iterations. The immigrated chromosomes are joined with the subpopulation and take part in the selection. The selection process creates an offspring population of chromosomes. The selection is conducted on the basis of fitness function values. The better fitted chromosomes have the highest probability being in the offspring subpopulation. Then stop optimiza- tion condition is checked, and the algorithm iterates in the case when it is not fulfilled. Thefitness function is computedon thebasis of results obtained fromFEM analysis. TheMSC.Nastran was used in the numerical example. The isotropic material parameters were coded into the MSC.Nastran format and stored in a text file. Solutions to direct problems for RVE were performed next. The homogenized orthotropic material properties were calculated on the basis of results of FEManalyses. Thefitness function valuewas calculated on the basis of experimental and numerical orthotropic properties of the RVE model. 5. Numerical example of multiscale modeling of the bone Twomultiscalemodeling problems of the proximal femur bone are considered: analysis and identification. 864 T. Burczyński et al. Fig. 6. The distributed parallel evolutionary algorithm 5.1. Analysis of the multiscale model of bone The geometry of microstructure is created on the basis ofmicro CT scans. A single trabeculae is considered to be isotropic and its material properties are presented in Table 1. The material properties of the trabecular bone in themacromodel are orthotropic and obtained on the basis of FEManalysis of RVE micromodels. The cortical tissue properties are based on the properties given inWirtz et al. (2000) and are shown in Table 2. Table 1.Material properties of single trabeculae used inRVEmicrostructures analysis Material parameter Young’s modulus [MPa] Poisson’s ratio Value 3300.0 0.33 Thegeometry of theproximal femurboneand loading conditions are shown in Fig.7. The femur bone is loaded by a force of 1300N. The discrete model consists of tetrahedron finite elements and has 98000 DOF (degrees of fre- Multiscale modeling of osseous tissues 865 Table 2.Material properties of cortical tissue Material parameter Young’s modulus [MPa] Poisson’s ratio Value 12000.0 0.3 Fig. 7. Themacrostructure – proximal femur bone, (a) boundary conditions, (b) geometry andmesh edom). The geometry of microstructure of the trabecular bone is presented in Fig.8 and Fig.9. The discrete model consists of tetrahedron finite elements and has 66000 DOF. Thedisplacement boundary conditionswere applied in themicrostructures analysis of RVE. Fig. 8. The microstructure of the trabecular bone from femur head The homogenized material properties of the bone obtained from the RVE microstrucutre are presented in Table 3. The RVE is representative only near 866 T. Burczyński et al. Fig. 9. The mesh of RVEmicrostructure the area where the sample of the trabecular bone is extracted from the real bone. To extend the applicability of the obtained material parameters, they shouldbescaledon thebasis of thebonedensity.Thisprocedurecanbeapplied only for bone areas with the similar geometrical structure. Themodel can be significantly improved whenRVEs are created for characteristic geometries in manyareas of thebone.Toprepare suchRVEs,many samples forCTscanning have to be extracted from the real bone. Table 3.Material properties of bone derived fromRVEmicro-model Material c11 c12 c13 c22 c23 c33 c44 c55 c66 parameter Value [MPa] 1002.0 321.0 239.0 902.0 239.0 604.0 633.0 457.0 457.0 The Henckey-von Mises-Huber stresses distribution for the applied load and for the macromodel of the bone are shown in Fig.10. 5.2. Identification of isotropic material properties of a single trabeculae The isotropic material properties for trabeculae can be found on the basis of known orthotropic material properties of the microstructure model. The orthotropic parameters for the microscale can be obtained by performing a tensile test on a trabeculae specimen or on the basis of ultrasonic velocityme- asurements (Trębacz and Gawda, 2001). The orthotropic material properties in themicroscale can be also acquired byperforming identification on the level of macro model. The identification ofmaterial parameters E and ν of the single trabeculae on the micro level has been performed as minimization of the objective func- Multiscale modeling of osseous tissues 867 Fig. 10. HMH stresses in the macromodel of the bone tion F given by (4.2), through the distributed parallel evolutionary algorithm with the following data: • number of subpopulations – 2 • number of chromosomes in each subpopulation – 15 • number of genes – 2 • evolutionary operators: – uniformmutation (10%) – simple crossover +Gaussian mutation (90%) – ranking selection • migration of the best chromosome – every 2 generations • number of iterations – 35 • number of processors used during test – 4 • computations time – 14h Themicrostructure model shown in Section 5.1 was used for computations. The numerical results of identification are presented in Table 4. The history of the objective function F changing during the optimization for two populations is presented inFig.11. Improvements of objective function values in subpopulations after the migration phase can be observed. 868 T. Burczyński et al. Table 4. Actual and foundmaterial parameters of trabecular bone in micro- scale Material parameters Actual Found Error [%] E [MPa] 3300.0 3305.5 0.16 ν 0.330 0.329 0.30 Fig. 11. History of the objective function for two subpopulations Avery good agreement can be seen between the actual and foundmaterial parameters. The identification problem in the multiscale modeling belongs to the newly emergingmethodology which is very useful in the determination of some material parameters in the microscale having information about some measurements from the macroscale. 6. Conclusions The methodology of multiscale modeling of the bone tissue was presented. Two problems of the multiscale modeling were considered for the proximal femur bone: analysis and identification. The studied problemwas solved using FEM and computational homogenization for the fully orthotropic elastic ma- terial with 9 independent material parameters. The model of the trabecular bone microstructure was prepared on the base of micro-CT scans and used duringmultiscale computation. The hierarchical structure of the cortical bone could be taken into account in the similarway. The identification problemwas considered which enabled determination of material parameters of the isotro- pic single trabeculae on the microscale having information about orthotropic parameters in the macroscale. This problem was solved as the minimization task using evolutionary computing. Multiscale modeling of osseous tissues 869 The presentedmethodology can be used in analysis and identification pro- blems of different kinds of bone tissues. Themultiscale approach implemented to the analysis of osseous tissue is very challenging and opens new possibili- ties to the modeling of processes which take place in biological tissues with hierarchical structures. Acknowledgements The research work was financed by the Polish science budget resources in the years 2008-2010 in the frame of a research project. References 1. Agić A., Nikolić V., Mijović B., 2006, The cancellous bone multiscale morphology-elasticity relationship,Collegium Antropologicum, 30, 2, 409-414 2. BąkR.,BurczyńskiT., 2009,Computational Strength ofMaterials (inPolish Wytrzymałość materiałów z elementami ujęcia komputerowego), WNT, War- szawa 3. Burczyński T., 2010, Evolutionary and immune computations in optimal de- sign and inverse problems, Chapter 2 in:Advances of Soft Computing in Engi- neering, Z.Waszczyszyn (Edit.), Springer, 57-132 4. BurczyńskiT.,KuśW., 2009,Microstructureoptimizationand identification in multi-scale modellig, Chapter in: New Computational Challenges in Mate- rials, Structures and Fluids, J. Eberhadstener et al. (Edit.), Springer, 169-181 5. Burczyński T., Mrozek A., Górski R., Kuś W., 2010, Molecular statics coupled with the subregion boundary element method in multiscale analysis, Int. Journal for Multiscale Computational Engineering, 8, 3, 319-331 6. Ghanbari J., NaghdabadiR., 2009,Nonlinear hierarchicalmultiscalemode- ling of cortical bone considering its nanoscale microstructure, Journal of Bio- mechanics, 42, 1560-1565 7. HamedE., LeeY., Jasiuk I., 2010,Multiscalemodeling of elastic properties of cortical bone,Acta Mechanica (online) 8. Ilic S., Hackl K., Gilbert R., 2010, Application of the multiscale FEM to the modeling of cancellous bone,Biomech. Model Mechanobiol., 9, 87-102 9. Kouznetsova V.G., 2002,Computational homogenization for the multi-scale analysis of multi-phase materials, Ph.D. Thesis, TU Eindhoven. 10. KowalczykP., 2010, Simulation of orthotropicmicrostructure remodelling of cancellous bone, Journal of Biomechanics, 43, 563-569 870 T. Burczyński et al. 11. Madej Ł., Mrozek A., Kuś W., Burczyński T., Pietrzyk M., 2008, Concurrent and upscaling methods in multi scale modelling – case studies, Computer Methods in Material Science, 8, 1, AGH, Kraków 12. Sansalone V., Lemaire T., Naili S., 2009, Variational homogenization for modeling fibrillar structures in bone,Mechanics ResearchCommunications,36, 265-273 13. Schneider R., Faust G., Hindenlang U., Helwig P., 2009, Inhomogene- ous, orthotropicmaterialmodel for the cortical structure of longbonesmodeled on the basis of clinical CT or density data,Comput. Methods Appl. Mech. En- grg., 1298, 2167-2174 14. Terada K., Kikuchi N., 2001, A class of general algorithms for multi-scale analyses forheterogeneousmedia,ComputerMethods inAppliedMechanics and Engineering, 190, 5427-5464 15. TrębaczH.,GawdaH., 2001,The estimation of structural anisotropyof tra- becular and cortical bone tissues based on ultrasonic velocity and attenuation, Acta of Bioengineering and Biomechanics, 3, 2, 41-48 16. Tsubota K., Adachi T., Nishiumi S., Tomita Y., 2003, Elastic properties of single trabeculae measured by micro-three-point bending test, Proc. of the International Conference on Advanced Technology in Experimental Mechanics 17. Wirtz D.C., Schiffers N., Pandorf T., Radermacher K., Weichert D., Forst R., 2000, Critical evaluation of known bone material properties to realize anisotropic FE-simulation of the proximal femur, Journal of Biomecha- nics, 33, 1325-1330 18. Zienkiewicz O.C., Taylor R.L., Zhu J.Z., 2005,The Finite Element Me- thod: Its Basis and Fundamentals, 6th Edition, Butterworth-Heinemann Modelowanie wieloskalowe tkanki kostnej Streszczenie W artykule przedstawionometodologię wieloskalowegomodelowania thanki kost- nej, w której zagadnienie identyfikacji parametrówmateriałowych odgrywa kluczową rolę. Rozpatrzono analizę dwuskalową kości, a problem identyfikacji sformułowano jako zagadnienie odwrotne, będąceważnymetapemprocesumodelowania. Jako przy- kład tkanki kostnej rozważono kość udową zbudowaną z kości gąbczastej i korowej. Manuscript received April 21, 2010; accepted for print June 25, 2010