Evolutionary Model for Bone Adaptation Using Cellular Automata

1. UNIVERSITY OF NOTRE DAME V I T A C E D O - D U L S P E S t = 0 t = 6 t = 3 t = 9 t = 15 t = 12 t = 18 t = 24 t = 21 Evolutionary Model for Bone Adaptation Using Cellular Automata University of Notre Dame – National University of Colombia – COLCIENCIAS Andrés Tovar1Neal Patel1 Gabriel A. Letona2 Glen L. Niebur3 Mihir Sen4 John E. Renaud4 Department of Aerospace and Mechanical Engineering University of Notre Dame, Notre Dame, Indiana 46556, USA Phone: (574)631-8616 Fax: (574)631-8341 E-mail: gniebur@nd.edu - jrenaud@nd.edu (A) (B) (C) (D) Figure 1 – Hierarchical model at the continuum and tissue level. A coronal section from the proximal femur (A) is modeled in finite-element model package FEMLAB (B). The model is divided into a lattice of sub-models of 1x1 mm2 each. The stress/strain field operating on the sub-models is determined by finite element analysis in the continuum model. Initially, each sub-model is a solid piece of material composed of 25x25 cellular automata. The tissue-level CA algorithm dictates the trabecular structure within the sub-models and hence their anisotropic properties and apparent density at the continuum level (C). A detailed view of the femoral neck shows the trabecular architecture simulated in the CA sub-models (D). Motivation and Objectives Simulations of bone functional adaptation can improve our understanding of age-related bone loss and be used in design of tissue-engineered devices. Recently developed simulation models incorporate traditional optimization algorithms to account for the anisotropic nature of trabecular bone. Despite some mechanical aspects reflected by the idealized anisotropic structure, it only represents a mathematical abstraction of the biological architecture of trabecular bone. The objective of this investigation is to develop an algorithm that simulates bone functional adaptation at the tissue level and incorporates biological aspects of bone remodeling. This algorithm is incorporated into a hierarchical model in order to determine apparent density distribution at the continuum level. The technique makes use of finite element analysis, cellular automaton (CA) paradigm, control and structural optimization techniques. Materials and Methods The hierarchical algorithm developed in this investigation includes a continuum-level (global) model that is divided into tissue-level sub-models. The stress/strain field evaluated in the global model defines the load condition in the sub-models. Cellular automata comprising a sub-model modify the topology by processes of formation/resorption on a pre-determined exterior surface. During this process, new surfaces are created resulting in formation of trabeculae. The material is characterized by a continuous volume fraction that varies from 0 (bone marrow) to 1 (fully mineralized tissue). The elastic properties are quantified using a power-law relationship between Young’s modulus and volume fraction. Each automaton in the sub-model contains an osteocyte that senses the mechanical stimuli. The size of an automaton is calculated according to the density of 12,000 to 20,000 osteocytes per cubic millimeter. A cubic millimeter of mineralized tissue is modeled as a lattice of 25 x 25 x 25 = 15,625 cellular automata. The communication network among osteocytes is modeled according to cellular automaton principles. Each CA is assumed to receive information from its immediate external environment and also from its neighbors within a certain proximity. In a 2D model an automaton is assumed to establish communication with up to 8 cells (Moore neighborhood), while in 3D the CA communicates with up to 18 neighbors. The mechanical stimulus perceived by the sensors is equivalent to the averaged strain energy density within the neighborhood. The evolutionary remodeling process utilizes a local control rule of functional adaptation that drives the averaged strain energy density to a target value. The control strategy includes proportional, integral and derivative actions. The rate of change in mass is expressed by The remodeling is constrained to the free surfaces of the mineralized tissue. When convergence is achieved, the density information and the anisotropic properties of the sub-models are fed back to the continuum model. Results and Conclusions A hierarchical functional adaptation model has been developed in this investigation. A continuum-level model evaluates the global stress/strain field. Tissue-level sub-models use the global field values to define local load conditions. Cellular automata in the sub-models apply a control rule to drive an averaged strain energy density to a target value. The convergent results are returned to the global model. Preliminary results show the application of this technique to quantify apparent density distribution in a two-dimensional femoral head (Fig. 1) and to simulate the formation of trabecular structure in a three-dimensional model (Fig. 2). Figure 2 – Evolutionary cellular automaton model at the tissue level. The cubic trabecular model of 8 mm3 is comprised of 125,000 cellular automata. Three rotational moments around the Cartesian axes are applied as distributed forces over the six surfaces of the cubic model. The structure is supported on the lower edge. Trabecular architecture results from formation and resorption on the surfaces during the evolutionary adaptation process. This figure shows some of the first 24 iterations. 1 Graduate Research Assistant 2 Research Assistant 3 Assistant Professor 4 Professor 14th conference of the European Society of Biomechanics July 4 – 7, 2004