Table of Contents Author Guidelines Submit a Manuscript
Computational and Mathematical Methods in Medicine
Volume 2012 (2012), Article ID 748302, 14 pages
Research Article

Spongiosa Primary Development: A Biochemical Hypothesis by Turing Patterns Formations

1Grupo de Investigación en Estudios y Aplicaciones de Ingeniería Mecánica (GEAMEC), Universidad Santo Tomás, Bogotá, Colombia
2Grupo de Modelado y Métodos Numéricos en Ingeniería (GNUM), Departamento de Ingeniería Mecánica y Mecatrónica, Facultad de Ingeniería, Universidad Nacional de Colombia, Bogotá, Colombia

Received 2 April 2012; Revised 24 May 2012; Accepted 5 June 2012

Academic Editor: Jérôme Noailly

Copyright © 2012 Oscar Rodrigo López-Vaca and Diego Alexander Garzón-Alvarado. 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.


We propose a biochemical model describing the formation of primary spongiosa architecture through a bioregulatory model by metalloproteinase 13 (MMP13) and vascular endothelial growth factor (VEGF). It is assumed that MMP13 regulates cartilage degradation and the VEGF allows vascularization and advances in the ossification front through the presence of osteoblasts. The coupling of this set of molecules is represented by reaction-diffusion equations with parameters in the Turing space, creating a stable spatiotemporal pattern that leads to the formation of the trabeculae present in the spongy tissue. Experimental evidence has shown that the MMP13 regulates VEGF formation, and it is assumed that VEGF negatively regulates MMP13 formation. Thus, the patterns obtained by ossification may represent the primary spongiosa formation during endochondral ossification. Moreover, for the numerical solution, we used the finite element method with the Newton-Raphson method to approximate partial differential nonlinear equations. Ossification patterns obtained may represent the primary spongiosa formation during endochondral ossification.

1. Introduction

Most of the long bones of the mammals skeletal system are developed from a process called endochondral growth [14]. This process ends with the gradual production of bone from cartilage tissue during fetal development and postnatal growth. The process of ossification occurs from a hyaline cartilage tissue mold, which has a similar shape to the bone in a mature stage. The cartilage tissue molds are formed through the condensation of mesenchymal cells [5] followed by their differentiation into chondrocytes (cells that produce and maintain cartilage matrix) and secretion of typical components of the extracellular matrix of cartilage [6]. Once the mold of cartilage is formed, it is invaded initially in its center and then at each end by a mixture of cells that give rise to the appearance of primary and secondary centers of ossification, respectively, [79]. The ossification centers invade the cartilage gradually until it is completely replaced by bone tissue, except the articular surfaces. In this way, and eventually the bones reach their skeletal maturity [10]. The processes of endochondral development, growth, and elongation of the bones are made by the continuous addition of cartilage and subsequent replacement by bone tissue.

During the chondrocytes differentiation process, the matrix composition changes dramatically through the production of other components such as collagen type X, the expression of metalloproteinases, and subsequent calcification. At the same time the blood vessels invade the calcified cartilage bringing osteoblasts which build immature bone. Chondrocytes in the growth plate are subjected to the influence of excess extracellular factors including systemic and soluble local factors, as well as, extracellular matrix components. Several studies [9, 1113] provide evidence that the proliferation of chondrocytes in the growth plate is under the control of a local closed loop that depends on the spatial and temporal location; and that mainly involves the molecular signals synthesized by chondrocytes: parathyroid hormone-related peptide (PTHrP), Indian hedgehog (Ihh), transforming growth factor (TGFβ), bone morphogenetic proteins (BMPs), vascular endothelial growth factor (VEGF), matrix metalloproteinase type 9, known as gelatinase-B (MMP9), and the transcription factor RUNX2. They interact together in a feedback loop to regulate the rate at which chondrocytes leave the proliferative zone, differentiate hypertrophic cells and give way to immature bone formation [10, 11, 14]. The inappropriate balance in the expression of these molecules along with the ones that encode the collagens and other growth factors have been subject of studies as possible causes of impaired bone formation by the endochondral ossification mechanism [1517].

The process of endochondral ossification has been studied for several years, and different models have been developed in silico, verified by histological reports, and in vivo experiments, which have tried to explain the process of bone formation through this mechanism [79, 13, 14, 1821]. For example, Courtin et al. [18] in their work performed the comparison between the sequence of morphological events involved in embryonic bone formation and spatiotemporal characteristics of self-organization generated by a reaction-diffusion model related to the metabolism of the periosteal bone mineralization. In that article, 3D structures are obtained (by computer simulation) with a close resemblance to the primary internal architecture of the periosteum of long bones. The hypothesis of Courtin et al. is based on the role of self-organization of mineralization of bone metabolism, which gives rise to a well-organized space architecture. Subsequent research such as those by Garzón-Alvarado et al. [79, 22] have raised different hypotheses about the interaction of mechanical, cellular, and molecular factors that lead to the formation of secondary ossification centers in the epiphyses of long bones, and it also helped the development and bone growth and the primary bone formation. These hypotheses suggest that biological processes and interactions between different factors can be represented using mathematical models, where the chemical feedback among molecular reagent factors through reaction-diffusion mechanisms may explain the stable spatial patterns found in the origin of the secondary ossification centers and in the formation of cartilage canals. As far as the authors know, no research has been conducted on the action of different cellular, mechanical, and molecular factors on the development and production of primary spongiosa architecture during the endochondral ossification process, which is the basis for the production of trabecular bone. Similarly, there are no biochemical models involving reaction-diffusion systems with Turing instabilities and reaction equations based on the Schnakenberg model, that allow enlarge the knowledge and understanding of the development in the primary stage of the trabecular bone.

Therefore, this paper presents a hypothesis on the development of trabecular bone architecture. Starting from the assumption that the interaction of two molecular factors expressed by hypertrophic chondrocytes which through a reaction-diffusion mechanism generate a stable spatial-temporal pattern. This patterns lead to the formation of trabeculae present in the primary spongiosa tissue. This is the first model that attempts to explain the formation of primary spongiosa, which serves as the basis for defining a complete model of trabeculae formation at an early stage of skeletal development.

2. Molecular Mechanisms Present in the Endochondral Ossification

2.1. Molecular Factors

The sequential changes in the behavior of chondrocytes in the growth plate are highly regulated by systemic factors and the production of local factors. Growth hormone (GH) and thyroid hormone are systemic factors involved in endochondral ossification that regulate the behavior of the chondrocytes. GH is a peptide hormone-based protein that stimulates growth, cell reproduction and tissue regeneration. It is an important regulator of longitudinal bone growth [23]. The main effect of GH on chondrocytes, it is to stimulate their proliferation [24]. However, thyroid hormone is considered another systemic regulator of bone growth, and it stimulates production of collagen type II and X and alkaline phosphatase (ALP), which act as markers of bone mineralization.

On the other hand, local factors act as receptors to carry out intracellular signaling and selective activation of transcription factors of chondrocytes, such as insulin-like growth factors (IGF), which act as a local mediator of the effects of GH on cartilage growth. These factors are essential for embryonic skeletal development [25, 26], and in chondrocyte proliferation and/or hypertrophy. Parathyroid hormone-related peptide (PTHrP) is expressed by perichondrial cells at the initial stage of the chondrocytes proliferation. The PTHrP diffuses out of its place of production to act on cells carrying the receptor PTH/PTHrP [27]. The PTHrP keeps chondrocytes in a proliferative state and prevents hypertrophy [28]. The Indian Hedgehog (Ihh) is a local factor produced by the expression of prehypertrophic chondrocytes that stimulates the proliferation of chondrocytes and inhibits its hypertrophy. [14]. The Ihh stimulates osteoblastic differentiation of mesenchymal cells, which is essential for the formation of the periosteum surrounding the zone of hypertrophic chondrocytes. The formation of the periosteum precedes the formation of primary ossification center and it is maintained through its expansion [13, 29]. Bone morphogenetic proteins (BMPs) are another local factor. These are members of the transforming growth factor beta (TGFβ superfamily), which is capable of inducing strongly immature bone formation, cartilage, and connective tissue [30]. Finally, within the local factors, there is the vascular endothelial growth factor (VEGF), which stimulates the process known as angiogenesis and acts as a vasodilator by increasing vascular permeability [31]. VEGF acts on vascular endothelial cells through specific tyrosine kinase membrane receptors, thereby regulating functions such as proliferation, differentiation, and migration of chondroblasts, osteoblasts, and osteoclasts [14]. During chondrocytes hypertrophy in the growth plate, VEGF is released and the extracellular matrix surrounding the hypertrophic cells, which begin the process of calcification. Then the extracellular matrix is invaded by blood vessels which provide nutrients and attract osteoblasts and osteoclasts that help the formation of the trabecular bone [32].

Transcription factors are mostly specific to a particular cell lineage and act as growth regulators of cell differentiation. They are predominantly expressed during skeletal development and its main function is to control cell proliferation or differentiation [33]. The transcription factor Runx2 is also called Cbfa1/Osf2/AML3/Til1/PPB2αA and is an essential protein in the differentiation of chondrocytes and osteblasts as well as the morphogenesis of the skeletal system [34]. The Runx2 controls bone mineralization of growing bones by stimulating osteoblast differentiation, promoting chondrocyte hypertrophy, and contributing to the migration of endothelial cells and vascular invasion. The Runx2 is expressed by chondrocytes in the early stages of hypertrophy and is maintained until terminal hypertrophic differentiation [35].

2.2. Regulation of Cartilage Matrix Degradation during Endochondral Ossification

The increase in cell volume experienced by the chondrocytes submitted to hypertrophy requires the degradation of the matrix that surrounds these cells. Moreover, the invasion of the ossification front requires an extensive (but selective) degradation of cartilage transverse columns surrounding the hypertrophic cells in the final state [10, 28]. There have been several studies to identify the proteolytic enzymes responsible of these events of matrix degradation and the cells responsible for its synthesis [3638]. These studies have emphasized on the enzymes capable to degrade the two major protein components of the cartilage matrix, collagen type II and aggrecan. Within the growth plate is expressed selectively MMP13 by hypertrophic chondrocytes, which degrades collagen fibers and aggrecan [38, 39].

The MMP13 expression by chondrocytes is a prerequirement for the invasion of the growth plate by blood vessels, osteoclasts, and osteogenic cells. Therefore, these cells cannot enter in the empty gaps created by the death of hypertrophic chondrocytes until it degrades the septa of the cavities by the MMP13 [10]. Blood vessels invade the growth plate at the same time that the osteoblasts do it, which are necessary for the establishment of the primary ossification center. Thus, in the absence of MMP13 most of the cartilage matrix is not removed and there is no cell invasion in the bone marrow or bone matrix deposition in the remaining cartilage. The vascular invasion of the growth plate is facilitated by the vascular endothelial growth factor (VEGF), which is expressed by chondrocytes, and regulated with the hypertrophy under the control of Runx2 [10, 11]. On the other hand, the VEGF in the endochondral ossification, increases bone formation and decreases bone resorption [40, 41], indicating that VEGF regulates the production of MMP13.

3. Materials and Methods

3.1. Hypothesis Required for the Development of Primary Spongiosa Using Reaction-Diffusion Systems

The main hypothesis of this paper is based on the existence, within the endochondral ossification process, of the controlled interaction of two signaling molecules that diffuse and react chemically in the cartilage extracellular matrix, to carry out the formation of primary spongiosa from the growth plate. Accordingly, we assume the existence of a reaction-diffusion system where two primary molecules are involved, such as VEGF and MMP13, which can lead to a stable pattern in time and unstable in space, similar to the patterns present in the structure of the trabecular bone during endochondral ossification.

The presence of MMP13, which is released by hypertrophic chondrocytes, allows the degradation of cartilage matrix components (collagen and aggrecan) and leads to vascular invasion and ossification front [10, 36, 37]. This vascular invasion is facilitated by the presence of VEGF expressed by hypertrophic chondrocytes [31, 40, 42]. This means that when MMP13 and VEGF exist in all regions of the epiphyseal cartilage, having a high concentration of VEGF, there will be an adequate control of the invasion of endothelial cells, osteoclasts, chondroclasts, and osteoblasts, which are present in the primary ossification development [40]. Similarly, in those areas where there is a high concentration of MMP13, it will completely degrade the cartilage, giving rise to the trabecular bone architecture. Last statements are supported on Hiltunen et al. [43] studies. In their work it was injected a saline solution containing VEGF in the distal femur of white rabbits. Their results demonstrate that VEGF induces bone formation by increasing osteoblast activity and decreasing the resorption process. The resorption process is produced by both osteoclasts in bone and metalloproteinases (MMPs) in the growth plate. Therefore, it can be supposed that in the development of the architecture of the primary spongiosa, there must exist a regulation of MMP13 by the VEGF (inhibitory mechanism). So, it stops the degradation of the cartilage and begins the invasion front of ossification.

3.2. Model Description

The regulatory process proposed in this model is outlined in Figure 1, and it is based on an activator-substrate system (also called exhaustion model) (see Appendix A). The process indicates that there is a control loop between VEGF (activating factor) and the MMP13 (substrate), where the VEGF is self-activated and inhibits the production of MMP13 stopping the degradation process and giving a way to the mineralization of the remaining cartilage matrix [40]. On the other hand, we assume that the MMP13 is self-inhibited, but enables the production of VEGF, thereby this loop is called positive feedback system. The VEGF helps vascular invasion and bring with it the osteogenic cells that allow the construction of the primary spongiosa. The MMP13 allows the degradation of the cartilage matrix and the subsequent invasion of the cartilage by the ossification front.

Figure 1: Control system of the molecular process. VEGF: vascular endothelial growth factor; MMP13: matrix metalloproteinase 13, the figure shows the ratio of the molecular signals produced by the hypertrophy of chondrocytes. The solid lines mean activation, dotted lines inhibition, continuous curved lines self-activation, and dotted curved lines self-inhibition.

The regulatory mechanism is modeled by reaction-diffusion equations. The reaction term (synthesis of soluble extracellular factors) is considered dependent on the concentration of the reactants and the presence of hypertrophic chondrocytes. According to this, the hypothesis is based on the origin of the patterns presented in the primary spongiosa. It could correspond, from a mathematical point of view, to the patterns that occur in the Turing space when two chemical reactants interact.

The definition of the relationships shown in Figure 1 can be quantified by means of equations which provide local changes of the extracellular soluble factors and the production rate of bone:𝜕𝑆VEGF𝜕𝑡=𝐶CH𝛼1𝜇𝑆VEGF+𝛾0𝑆2VEGF𝑆MMP13+𝐷VEGF2𝑆VEGF,(1a)𝜕𝑆MMP13𝜕𝑡=𝐶CH𝛼2𝛾0𝑆2VEGF𝑆MMP13+𝐷MMP132𝑆MMP13,(1b)𝜕𝑐Bone𝜕𝑡=𝜂𝑆𝑛VEGF𝑆𝑛VEGF+𝑆𝑛umbral𝑇𝑟𝑎𝑇𝑟𝑎+𝑡𝑟,(1c)where 𝐶CH is the concentration of hypertrophic chondrocytes, 𝑆VEGF, 𝑆MMP13 represent the concentrations of VEGF and MMP13, respectively. The remaining terms are model parameters. 𝛼1 and 𝛼2 are terms that quantify the production of each molecular factor by hypertrophic chondrocytes, 𝜇 is a constant that quantifies inhibition of the VEGF production by it excess, 𝛾0 regulates the nonlinear interaction between the concentration of MMP13-VEGF quantifying the concentration or molecular inhibition of each molecular factor, and 𝐷VEGF and 𝐷MMP13 are the diffusion coefficients of VEGF and MMP13, respectively. In the biological interpretation of the above equations the term 𝛾0𝑆2VEGF𝑆MMP13 represents the nonlinear activation of 𝑆VEGF (production of VEGF by the presence of MMP13) and the nonlinear consumption of 𝑆MMP13 (by the presence of VEGF). Equation (1c) represents the activation of bone production rate by the presence of high amounts of VEGF, which is regulated as time goes on. In this equation 𝑐Bone indicates the production of bone per unit of volume due to the concentration and distribution of VEGF within the domain. 𝜂 is a constant that regulates the production of bone over time, 𝑆𝑛Umbral represents the value of the concentration of VEGF with which begins the process of ossification, 𝑇𝑎 is the time required for the process of cartilage calcification, and 𝑡𝑟 represents the time that limits the production of bone.

3.3. Solution of the Reaction-Diffusion Equations System Using the Finite Element Method

To solve the set of ((1a), (1b), and (1c)), we used the finite element method, with tetrahedral elements. Due to the nonlinearity of the terms included in the model, we used the Newton-Raphson method to solve the problem of time evolution of the concentration of VEGF (𝑆VEGF) and MMP13 (𝑆MMP13). The time integration was performed using the trapezoidal rule.

3.3.1. Weak Formulation

Let ((1a), (1b), and (1c)) be rewritten as 𝜕𝑆VEGF𝜕𝑡𝐶CH𝛼1𝜇𝑆VEGF+𝛾0𝑆2VEGF𝑆MMP13𝐷VEGF2𝑆VEGF=0,𝜕𝑆MMP13𝜕𝑡𝐶CH𝛼2𝛾0𝑆2VEGF𝑆MMP13𝐷MMP132𝑆MMP13=0,𝜕𝑐Bone𝜕𝑡𝜂𝑆𝑛VEGF𝑆𝑛VEGF+𝑆𝑛umbral𝑇𝑟𝑎𝑇𝑟𝑎+𝑡𝑟=0.(2)

Using weighted residuals, the system of (2) takes the form Ω𝑤1𝜕𝑆VEGF𝜕𝑡𝐶CH𝛼1𝜇𝑆VEGF+𝛾0𝑆2VEGF𝑆MMP13𝐷VEGF2𝑆VEGF𝑑Ω=0,Ω𝑤2𝜕𝑆MMP13𝜕𝑡𝐶CH𝛼2𝛾0𝑆2VEGF𝑆MMP13𝐷MMP132𝑆MMP13𝑑Ω=0,Ω𝑤3𝜕𝑐Bone𝜕𝑡𝜂𝑆𝑛VEGF𝑆𝑛VEGF+𝑆𝑛umbral𝑇𝑟𝑎𝑇𝑟𝑎+𝑡𝑟𝑑Ω=0,(3) where Ω represents the domain of the problem which is limited by the boundary Γ. 𝑤1, 𝑤2, and 𝑤3 are weight functions. Using Green's theorem and weakening the system of (3), the residue of the problem is given by 𝑟𝑆VEGF=Ω𝑤1𝜕𝑆VEGF𝜕𝑡𝑑ΩΩ𝑤1𝐶CH𝛼1𝑑Ω+Ω𝑤1𝐶CH𝜇𝑆VEGF𝑑ΩΩ𝑤1𝐶CH𝛾0𝑆2VEGF𝑆MMP13𝑑Ω+Ω𝐷VEGF𝑤1𝑆VEGF𝑑ΩΓ𝑤1𝐷VEGF𝑆VEGF𝑛𝑑Γ=0,𝑟SMMP13=Ω𝑤2𝜕𝑆MMP13𝜕𝑡𝑑ΩΩ𝑤2𝐶CH𝛼2𝑑Ω+Ω𝑤2𝐶CH𝛾0𝑆2VEGF𝑆MMP13𝑑Ω+Ω𝐷MMP13𝑤2𝑆MMP13𝑑ΩΓ𝑤2𝐷MMP13𝑆MMP13𝑛𝑑Γ=0,𝑟𝐶Bone=Ω𝑤3𝑑𝑐Bone𝑑𝑡𝜂𝑆𝑛VEGF𝑆𝑛VEGF+𝑆𝑛umbral𝑇𝑟𝑎𝑇𝑟𝑎+𝑡𝑟𝑑Ω=0.(4)

Defining null flow conditions on the contour of the problem, the border-related terms in (4) are canceled, so the residue is expressed as 𝑟𝑆VEGF=Ω𝑤1𝜕𝑆VEGF𝜕𝑡𝑑ΩΩ𝑤1𝐶CH𝛼1𝑑Ω+Ω𝑤1𝐶CH𝜇𝑆VEGF𝑑ΩΩ𝑤1𝐶CH𝛾0𝑆2VEGF𝑆MMP13𝑑Ω+Ω𝐷VEGF𝑤1𝑆VEGF𝑑Ω=0,𝑟𝑆MMP13=Ω𝑤2𝜕𝑆MMP13𝜕𝑡𝑑ΩΩ𝑤2𝐶CH𝛼2𝑑Ω+Ω𝑤2𝐶CH𝛾0𝑆2VEGF𝑆MMP13𝑑Ω+Ω𝐷MMP13𝑤2𝑆MMP13𝑑Ω=0,𝑟𝐶Bone=Ω𝑤3𝜕𝑐Bone𝜕𝑡𝜂𝑆𝑛VEGF𝑆𝑛VEGF+𝑆𝑛umbral𝑇𝑟𝑎𝑇𝑟𝑎+𝑡𝑟𝑑Ω=0.(5)

To discretize the finite element solution we use approaching functions written as a linear combination of orthogonal functions as shown in the following: 𝑆𝑒VEGF=𝐍𝑣𝐒VEGF,𝑆𝑒MMP13=𝐍𝑚𝐒MMP13,𝑐𝑒Bone=𝐍𝑐𝐜Bone,(6) where 𝐍𝑣, 𝐍𝑚, and 𝐍𝑐 represent the shape functions which depend only on the space used for the formulation, 𝐒VEGF and 𝐒MMP13 are the values of 𝑆VEGF and 𝑆MMP13 in the nodal points, and the superscript 𝑒 indicates the discretization of the finite element variable. For the weighting functions we used the Bubnov-Galerkin formulation, indicating that the functions 𝑤 take the same form of approximation functions 𝐍.

Substituting (6) in (5), we obtain the residual vector in its discrete form (7), where 𝑟𝑒VEGF, 𝑟𝑒MMP13 and 𝑟𝑒𝐶Bone are residue vectors for each equation and 𝐍 is the gradient vector of the approximation functions 𝑟𝑒𝑆VEGF=Ω𝐍𝑇𝜕𝑆𝑒VEGF𝜕𝑡𝑑ΩΩ𝐍𝑇𝐶CH𝛼1𝑑Ω+Ω𝐍𝑇𝐶CH𝜇𝑆𝑒VEGF𝑑ΩΩ𝐍𝑇𝐶CH𝛾0𝑆𝑒VEGF2𝑆𝑒MMP13𝑑Ω+Ω𝐷VEGF𝐍𝑇𝑆𝑒VEGF𝑑Ω,𝑟𝑒𝑆MMP13=Ω𝐍𝑇𝜕𝑆𝑒MMP13𝜕𝑡𝑑ΩΩ𝐍𝑇𝐶CH𝛼2𝑑Ω+𝐍𝑇+Ω𝑤2𝐶CH𝛾0𝑆𝑒VEGF2𝑆𝑒MMP13𝑑Ω+Ω𝐷MMP13𝐍𝑇𝑆𝑒MMP13𝑑Ω,𝑟𝑒𝐶Bone=Ω𝐍𝑇𝜕𝑐Bone𝜕𝑡𝜂𝑆𝑛VEGF𝑆𝑛VEGF+𝑆𝑛umbral𝑇𝑟𝑎𝑇𝑟𝑎+𝑡𝑟𝑑Ω.(7)

Using the time discretization equation based on the Crank-Nicolson equations, the equations in (7) are transformed into the following: 𝑟𝑒𝑆VEGF=Ω𝐍𝑇𝑆𝑒VEGF𝑆𝑒VEGFΔ𝑡𝑑Ω+𝛼Ω𝐷VEGF𝐍𝑇𝑆𝑒VEGF𝑑ΩΩ𝐍𝑇𝐶CH𝛼1𝑑Ω+Ω𝐍𝑇𝐶CH𝜇𝑆𝑒VEGF𝑑ΩΩ𝐍𝑇𝐶CH𝛾0𝑆𝑒VEGF2𝑆𝑒MMP13𝑑Ω𝑡+Δ𝑡+(1𝛼)Ω𝐷VEGF𝐍𝑇𝑆𝑒VEGF𝑑ΩΩ𝐍𝑇𝐶CH𝛼1𝑑Ω+Ω𝐍𝑇𝐶CH𝜇𝑆𝑒VEGF𝑑ΩΩ𝐍𝑇𝐶CH𝛾0𝑆𝑒VEGF2𝑆𝑒MMP13𝑑Ω𝑡,𝑟𝑒𝑆MMP13=Ω𝐍𝑇𝑆𝑒MMP13𝑆𝑒MMP13Δ𝑡𝑑Ω+𝛼Ω𝐷MMP13𝐍𝑇𝑆𝑒MMP13𝑑ΩΩ𝐍𝑇𝐶CH𝛼2𝑑Ω+𝐍𝑇+Ω𝑤2𝐶CH𝛾0𝑆𝑒VEGF2𝑆𝑒MMP13𝑑Ω𝑡+Δ𝑡+(1𝛼)Ω𝐷MMP13𝐍𝑇𝑆𝑒MMP13𝑑ΩΩ𝐍𝑇𝐶CH𝛼2𝑑Ω+𝐍𝑇+Ω𝑤2𝐶CH𝛾0𝑆𝑒VEGF2𝑆𝑒MMP13𝑑Ω𝑡,𝑟𝑒𝐶Bone=Ω𝐍𝑇̃𝑐Bone𝑐BoneΔ𝑡𝑑ΩΩ𝐍𝑇𝜂𝑆𝑛VEGF𝑆𝑛VEGF+𝑆𝑛umbral𝑇𝑟𝑎𝑇𝑟𝑎+𝑡𝑟𝑑Ω𝑡+Δ𝑡Ω𝐍𝑇𝜂𝑆𝑛VEGF𝑆𝑛VEGF+𝑆𝑛umbral𝑇𝑟𝑎𝑇𝑟𝑎+𝑡𝑟𝑑Ω𝑡,(8)

where 𝑆𝑒VEGF and 𝑆𝑒MMP13 are the nodal values 𝑆𝑒VEGF and 𝑆𝑒MMP13 evaluated over time 𝑡+Δ𝑡. 𝛼 is a parameter characteristic of the integration method.

Using (8) is possible to determine each of the terms of the tangent stiffness matrix, as shown in the following: 𝜕𝑟𝑒𝑆VEGF𝜕𝑆𝑒VEGF=1Δ𝑡Ω𝐍𝑇𝐍𝑑Ω+𝛼Ω𝐷VEGF𝐍𝑇𝐍𝑑Ω+Ω𝐶CH𝜇𝐍𝑇𝐍𝑑ΩΩ𝐍𝑇𝐍𝐶CH𝛾02𝑆𝑒VEGF𝑆𝑒MMP13𝑑Ω,𝜕𝑟𝑒𝑆VEGF𝜕𝑆𝑒MMP13=𝛼Ω𝐍𝑇𝐍𝐶CH𝛾0𝑆𝑒VEGF2𝑑Ω,𝜕𝑟𝑒𝑆MMP13𝜕𝑆𝑒VEGF=𝛼Ω𝐍𝑇𝐍𝐶CH𝛾02𝑆𝑒VEGF𝑆𝑒MMP13𝑑Ω,𝜕𝑟𝑒𝑆MMP13𝜕𝑆𝑒MMP13=1Δ𝑡Ω𝐍𝑇𝐍𝑑Ω+𝛼Ω𝐷MMP13𝐍𝑇𝐍𝑑Ω+Ω𝐍𝑇𝐍𝐶CH𝛾0𝑆𝑒VEGF2𝑑Ω,𝜕𝑟𝑒𝐶Bone𝜕̃𝑐Bone=1Δ𝑡Ω𝐍𝑇𝐍𝑑Ω.(9)

The nodal values 𝑆𝑒VEGF and 𝑆𝑒MMP13 in time 𝑡+Δ𝑡 can be approximated by the iterative algorithm of Newton-Raphson, as described in (10). In this equation, Δ𝑆𝑒VEGF and Δ𝑆𝑒MMP13 represent the difference of the nodal values in two consecutive iterations 𝜕𝑟𝑒𝑆VEGF𝜕𝑆𝑒VEGF𝜕𝑟𝑒𝑆VEGF𝜕𝑆𝑒MMP13𝜕𝑟𝑒𝑆MMP13𝜕𝑆𝑒VEGF𝜕𝑟𝑒𝑆MMP13𝜕𝑆𝑒MMP13Δ𝑆𝑒VEGFΔ𝑆𝑒MMP13=𝑟𝑒VEGF𝑟𝑒MMP13.(10)

3.4. Numerical Implementation

The set of ((1a), (1b), and (1c)) were implemented and solved numerically using the finite element method with a Newton-Raphson scheme. The two examples given were solved in a laptop of 4096 MB and 800 MHz processor speed. Computer simulation was carried out in an incremental iterative scheme that allows solving, computationally, the evolution of both the concentration of molecular factors (𝑆VEGF,𝑆MMP13) and the production of immature bone. Initially the growth plate is assumed as a structural matrix with an initial concentration of chondrocytes in a hypertrophic stage 65.000 cell/mm3. The initial concentrations of VEGF and MMP13 are distributed randomly in the growth plate, with a disturbance of 10% on the steady-state, the concentration is given by (𝑆VEGF,𝑆MMP13)=(1.0,0.9) [ng/mL] [43] (see Appendix A). The selection of random initial conditions around the steady state is similar to the event of molecular expression of the hypertrophic chondrocytes in an area of ossification. The flow condition for each molecular factor in the boundary is assumed to be null, this is because these conditions are assumed periodic over the domain. The parameter values used are shown in Tables 2(a) and 2(b), the justification of all the parameters used in the illustrated examples is given in Appendix B.

4. Results and Discussion

To verify the potential of the proposed model in predicting the primary spongiosa architecture, two numerical tests were conducted in a three-dimensional cubic element with a length of 0,22 mm. This is an average length of a minimal element of the trabecular bone in postmineralized fetal stage, as presented by Ruimerman et al. [44] in their work of modeling and adaptation of trabecular bone. The parameters of the reaction-diffusion model were selected in order to obtain structures with a periodicity in accordance with the one present in the trabecular bone [28, 44, 45]. In the finite element mesh we employed 17.756 nodes and 16.625 tetrahedral elements. In all the simulations we used incremental steps of Δ𝑡=0.1 therefore this means that the simulation time is measured in seconds and every step Δ𝑡 in the biological process corresponds to 7 seconds.

As a result of chemical interaction between the two molecular factors (reactants) and by the numerical results it was determined spatial patterns stable over time. The concentration of molecular factors in the cartilage and the action of the diffusive processes allows the formation of a pattern that is replicated throughout the domain. The architecture of the primary spongiosa obtained by the proposed RD model depends on the parameters used in ((1a), (1b), and (1c)) (see Tables 2(a) and 2(b)), so you can get structures with wave number (2,2,2) as shown in Figure 2 and with wave number (4,2,2) in Figure 3  (see Appendix B). The wave numbers allow either to define or to condition the frequency and distribution of the number of pores in a specific direction [9, 22, 46]. The appearance of either microstructure of the primary spongiosa, depends on the location of the parameters in the reaction diffusion equation ((1a), (1b), and (1c)) in the space of Turing. The location of certain points in this space, determine spatial patterns as shown in the results of this article.

Figure 2: (a) Architecture of the primary spongiosa obtained by the proposed R-D model with a wave number (2,2,2). (b) Spatial distribution of VEGF concentration at the end of the ossification process. (c) Spatial distribution of MMP13 concentration at the end of the ossification process.
Figure 3: (a) Architecture of the primary spongiosa obtained by the proposed R-D model with a wave number (4,2,2). (b) Spatial distribution of VEGF concentration at the end of the ossification process. (c) Spatial distribution of MMP13 concentration at the end of the ossification process.

The results in Figure 2 show the formation of two half-waves in each of the 𝑥, 𝑦, and 𝑧, directions, while Figure 3 shows the formation of four half-waves in the 𝑥 direction and 2 half-waves in the 𝑦 and 𝑧 directions, respectively. Figures 2(b), 2(c), 3(b), and 3(c) show the results for the organization of VEGF and MMP13 after the stabilization of the reaction-diffusion process, note that in the areas of greatest concentration of VEGF is produced cartilage calcification and in areas of greatest presence of MMP13 degradation occurs (empty space).

From the reaction-diffusion mechanism, it can be determined the change in the concentration of VEGF and MMP13 for each time step, as shown in Figure 4. The concentrations of VEGF (𝑆VEGF) and MMP13 (𝑆MMP13) within the cartilage evolve according to its diffusivity, its interaction, and its expression by the hypertrophic chondrocytes. So, the VEGF and the MMP13 are concentrated in high amounts in specific areas of the growth plate, allowing the formation of regular patterns similar to what happened in different biological models [9, 4751].

Figure 4: Temporal evolution of the concentration of VEGF (𝑆VEGF) by the reaction-diffusion mechanism. (a) Wave mode (2,2,2). (b) Wave mode (4,2,2).

The architecture of the primary spongiosa in a cubic element of length 0.44 mm is shown in Figure 5. In this figure we can observe the regular patterns for two different wave modes (2,2,2) and (4,2,2). Likewise, it details the advance in the front of ossification in different time steps, allowing the invasion of the cartilage by osteogenic cells that produce degradation and calcification. This will promote the formation of primary trabeculae, which subsequently undergo the bone remodeling processes, product of the stress distribution on bone tissue.

Figure 5: Primary spongiosa architecture produced by the interaction of VEGF and MMP13 by means of a reaction-diffusion mechanism. (a) Architecture generated by the use of a wave mode (2,2,2). (b) Architecture generated by the use of a wave mode (4,2,2). The blue areas represent ossified cartilage, and the red areas represent unossified cartilage.

The development of primary spongiosa has been the subject of study in recent years, but there is no clarity about the biological and mechanical factors affecting its formation. It was found that the structure of the trabecular tissue from its beginning (fetal age) to adulthood has different patterns in its formation, as well as, a variation in the trabecular density product of the bone remodeling process and the effect of mechanical loads. The proposed model assumes that the formation of patterns is due to the interaction through a reaction-diffusion system of two molecules (VEGF, MMP13) during endochondral ossification. The results presented show that the patterns self-organize along the domain used, as shown in Figures 2, 3, and 5. These structures represent the architecture of the primary spongiosa considering only biochemical effects. The results obtained in this work can be compared to the structure used by Ruimerman et al. [44], which, show self-organized repetitive patterns that serve as the basis for the maintenance and adaptation of mature trabecular structure.

The production of molecular factors that act as activator-substrate by the differentiation of prehypertrophic chondrocyte, are not necessarily the only factors expressed by these chondrogenic cells, which probably affect considerably the process of ossification, even the chondrocytes are not the only cells that act in this process. However, the proposed model only focuses on the formation of primary spongiosa architecture and not the entire calcification process, in which bone cells also operate as osteoclasts and osteoblasts. In the complete calcification process the model should not only incorporate chemical influences (bioregulatory model) but also involves loads and restrictions at the boundary (mechanical effects) as well as additional biochemical factors that must be taken into account. For example, Ruimerman et al. [44], Jang and Kim [52], Renders et al. [53], and Coelho et al. [54] have proved in their works that mechanical factors play an important role in the development, adaptation, and maintenance of the trabecular bone structure. Moreover, from the viewpoint of the biochemical factors, works such as those of Brouwers et al.[14] have evaluated the potential of three growth factors, PTHrP, Ihh, and VEGF, that interact and regulate the tissue differentiation and development of a long bone.

Much has been learned in recent years about the cellular and molecular mechanisms that guide the different events which allow the production of immature bone through endochondral ossification mechanism [6, 10, 11, 14, 19, 27, 31, 33, 45, 5558]. Nevertheless, there are still concerns about the relationship and interaction of different events to allow ossification and endochondral growth.

In this paper we presented a bioregulatory model based on a set of reaction-diffusion equations to predict the formation of primary spongiosa architecture. The application of the reaction-diffusion model with parameters in the Turing space is an area of constant work and controversy in biology. Garzón-Alvarado et al. [8, 9, 22, 46], Courtin et al. [18] and Cramping and Maini [47], used in their researches reaction-diffusion models to simulate different biological processes, finding in their results, that the use of these systems may explain many complex biological phenomena where pattern formation is a constant variable.

5. Conclusions

In this paper we presented the development of a biochemical model involving reaction-diffusion systems with instabilities in the Turing space. This model attempts to explain the generation of primary spongiosa in the process of endochondral ossification, an event that is not yet fully understood, due to the amount of biological, mechanical and biochemical effects that are present. The model involves the controlled interaction of two important molecular factors, such as VEGF and MMP13, present in the development and bone formation.

The work presented in this paper illustrates and supports the validity of the reaction-diffusion models to describe the processes occurring during a complex event of pattern formation in bone biology. From the results presented, it can be concluded that the chemical feedback between the two reactants molecular factors (activator-substrate) could be an explanation from a set of possible factors that determine the complex spatial patterns found in the origin of the architecture of the primary spongiosa. However, it is clear that these results have been obtained with a mathematical model based on assumptions and simplifications that should be discussed.

The hypothesis presented suggests that the origin of the primary spongiosa is internally controlled by cartilage cells. This is achieved through two biochemical reagents, VEGF and MMP13. These are not the only factors acting in the endochondral ossification, there are many others, among which count Ihh, PTHrP, Runx2, BMP, [10, 11, 14, 31] and likely influence similarly to the trabecular bone formation. Until now there has been a great effort to fully understand the role of each of these substances, how they interact and what processes they regulate. It is possible that VEGF-MMP13 are not the factors that control the entire process of endochondral ossification, but the existence of an activator-substrate mechanism ensures high stability for the development of this biological process. By the other side, the hypothesis raised is based entirely on the cycle of bone formation through the endochondral ossification mechanism, as presented in [6, 10, 28]. Afterwards, in the bioregulatory model (see Figure 1) it is taken the MMP13 as an agent responsible for the degradation of cartilage septa and as a promoter of the activation of VEGF, which fosters vascularization of the cartilage and subsequent calcification. However, we do not discard the possibility of a new bioregulatory model where it considers the counter case presented. This indicates that the activation of MMP13 may be produced by the action of VEGF, suggesting that cartilage degradation initially requires vascular invasion, as presented in the work of Pufe et al. [59].

In the development of the model it is assumed for the initial conditions that the activating factor is released by hypertrophic chondrocytes, as well as the substrate; however, the type of spatial instability obtained is independent of the initial conditions. Nevertheless, this model is very stable and robust with respect to the initial conditions and the range of parameters.

Finally, despite all the limitations and simplifications, the proposed mathematical model is able to reproduce in detail the architecture of the primary spongiosa, allowing variation in porosity and thickness of the trabeculae. The proposed model will serve as the basis for the formation of the secondary spongiosa architecture, from the bone remodeling process, observing the action of bone cells and the different mechanical effects that determine the orientation of the trabeculae.


A. Reaction-Diffusion Equations

Turing [60] proved, theoretically, that a chemical system of reaction diffusion could spontaneously, evolve in a heterogeneous spatial pattern from a uniform initial state in response to infinitesimal perturbations [51, 59]. His model is a system of two special partial differential equations, which are known as the two components of a reaction-diffusion system. In the most general case where several components are present the system has the form: 𝜕𝑢𝑖𝜕𝑡=𝐷𝑖2𝑢𝑖+𝑓𝑖𝑢1,𝑢2,,𝑢𝑚,𝑥Ω,𝑡0,(A.1) where the unknown functions 𝑢𝑖(𝑥,𝑡) can be interpreted as a reagent of concentration, the term 𝐷𝑖2𝑢𝑖 describes the diffusion reagent and 𝑓𝑖 is a smooth function (usually polynomial or rational in 𝑢𝑖) that describes a non-linear chemical interaction between the reagents. Ω refers to the bounded domain and in the system operates some initial boundary conditions. Turing also introduced some key concepts such as an activator and an inhibitor, so it was assumed that cellular states are discrete and can be modified by special chemical reagents.

It is known that for two components of a reaction-diffusion system, the Turing instability, which leads to the formation of patterns [60], it is possible to find it for two types of systems characterized by the signs of the stability matrix (Jacobian) in a positive homogeneous stationary system. The signs of the coefficients of the linearized dynamic system around a fixed point, namely the defined curves by the functions 𝑓(𝑢,𝑣)=0 and 𝑔(𝑢,𝑣)=0 in the plane 𝑢𝑣, for a system of two chemical components, are cut in the point (𝑢0,𝑣0), contain relevant information on the mechanisms of destabilization of the homogeneous solution. It can be shown that when the signs of the stability matrix are given by 𝐉𝑓𝑢𝑓𝑣𝑔𝑢𝑔𝑣=++TipoI𝑜++TipoII.(A.2)

It is possible to find a Turing instability if the reaction rates and diffusion coefficients allow the fixed point (𝑢0,𝑣0) to be stable to small homogeneous perturbations, and be able to become unstable to inhomogeneous perturbations, these conditions can be summarized mathematically by four inequalities: 𝑓𝑢+𝑔𝑣<0,𝑓𝑢𝑔𝑣𝑓𝑣𝑔𝑢>0,𝑑𝑓𝑢+𝑔𝑣>0,𝑑𝑓𝑢+𝑔𝑣24𝑑𝑓𝑢𝑔𝑣𝑓𝑣𝑔𝑢,(A.3) where 𝑑 represents the nondimensional coefficient of the diffusivities 𝑑=𝐷𝑣/𝐷𝑢. These inequalities are based on the linear stability analysis, that is, from the study of the eigenvalues (λ) of the linearized dynamics.

The stability matrix indicates what type of Turing instability can be found, for the case presented in this paper, the matrix corresponds to type 1 which identifies a reaction-diffusion system activator-substrate.

According to (A.1), for a reaction-diffusion system of two components, we have 𝜕𝑢𝜕𝑡=𝐷12𝑢+𝑓(𝑢,𝑣),𝜕𝑣𝜕𝑡=𝐷22𝑣+𝑔(𝑢,𝑣),(A.4) where we have the following terms of reaction: 𝑓(𝑢,𝑣)=𝑎𝑢+𝑢2𝑣,𝑔(𝑢,𝑣)=𝑏𝑢2𝑣.(A.5)

For the point of equilibrium, we have 𝑓(𝑢,𝑣)=0,𝑔(𝑢,𝑣)=0.(A.6)

The Jacobian (stabilization matrix) is given by 𝐉=𝜕𝑓𝜕𝑢𝜕𝑓𝜕𝑣𝜕𝑔𝜕𝑢𝜕𝑔𝜕𝑣.(A.7)

Therefore, 𝜕𝑓𝜕𝑢=2𝑢𝑣1,𝜕𝑓𝜕𝑣=𝑢2,𝜕𝑔𝜕𝑢=2𝑢𝑣,𝜕𝑔𝜕𝑣=𝑢2.(A.8)

For a steady point (𝑢0,𝑣0), we have 𝑓(𝑢,𝑣)=0 and 𝑔(𝑢,𝑣)=0; therefore, 𝑢 and 𝑣 are given by 𝑢=𝑎+𝑏,𝑣=𝑏(𝑎+𝑏)2.(A.9)

Replacing (A.9) into (A.8), we have 𝜕𝑓𝜕𝑢=2𝑏𝑎+𝑏1,𝜕𝑓𝜕𝑣=(𝑎+𝑏)2,𝜕𝑔𝜕𝑢=2𝑏𝑎+𝑏,𝜕𝑔𝜕𝑣=(𝑎+𝑏)2.(A.10)

Taking 𝑎=0.1 and 𝑏=0.9 [9] and replacing in (A.10), we have 𝜕𝑓𝜕𝑢=0.8,𝜕𝑓𝜕𝑣=1,𝜕𝑔𝜕𝑢=1.8,𝜕𝑔𝜕𝑣=1.(A.11)

Therefore, the stability matrix takes the following structure of signs: 𝐉𝑓𝑢𝑓𝑣𝑔𝑢𝑔𝑣=++.(A.12)

Indicating that the substance 𝑢 is self-activated and inhibits the substance 𝑣, while, the substance 𝑣 is self-inhibited and activates the substance 𝑢, forming an activator-substrate system.

B. Estimation of the Values for the Parameters

The set of ((1a), (1b), and (1c)) corresponds to a coupled system, where the equations that correspond to the molecular factors (1a) and (1b) are extended reaction-diffusion equations, similar to a Turing system, which has a controlled diffusion due to instabilities. For (𝐷VEGF,𝐷Ihh0), the controlled diffusion by instabilities appear for some combination of parameters [9, 47, 61], this defines a domain in the parameters space called Turing space. To get a Turing space is necessary a linear stability analysis of the reaction-diffusion system on the homogeneous solution, which is obtained by forcing (𝜕𝑆VEGF/𝜕𝑡)(𝐷VEGF=0)=0 and (𝜕𝑆MMP13/𝜕𝑡)(𝐷MMP13=0)=0, and leads to (𝑆VEGF,𝑆MMP13)=((𝛼1+𝛼2)/𝜇,𝛼2𝜇2/𝛾0(𝛼1+𝛼2)2). The linear analysis allows you to find the spatial patterns of the linearized solution and the range of parameters that ensure the emergence of such specific patterns [61]. Therefore, the solution can be expressed as (𝑆VEGF,𝑆MMP13)=(𝑢+𝑆VEGF,𝑣+𝑆MMP13), where 𝑢 and 𝑣 are small perturbations in each molecular factor, respectively. From (1a) and (1b) the results of the linear analysis allow to write the following inequalities: 𝐶CH2𝛾0𝑆VEGF𝑆MMP13𝛾0𝑆VEGF𝜇<0,𝐶2CH𝛾0𝑆2VEGF𝜇2𝛾0𝑆VEGF𝑆MMP13+2𝛾02𝑆3VEGF𝑆MMP13>0,𝐶CH𝐷MMP132𝛾0𝑆VEGF𝑆MMP13𝜇𝐷VEGF𝛾0𝑆2VEGF>0,𝐶2CH𝐷MMP132𝛾0𝑆VEGF𝑆MMP13𝜇𝐷VEGF𝛾0𝑆2VEGF24𝐷VEGF𝐷MMP13𝐶2CH𝛾0𝑆2VEGF×𝜇2𝛾0𝑆VEGF𝑆MMP13+2𝛾02𝑆3VEGF𝑆MMP13>0.(B.1)

These inequalities define a domain in the parameter space, known as the Turing space, where the uniform steady state (𝑆VEGF,𝑆MMP13) is linearly unstable.

If we express (1a) and (1b) in a nondimensional form (Schnakenberg equation) and as a function of small perturbations of the molecular factors (𝑆VEGF,𝑆MMP13), respectively, through (𝑢,𝑣), we can obtain 𝜕𝑢𝜕𝑡=𝛾𝑎𝑢+𝑢2𝑣+2𝑢,𝜕𝑣𝜕𝑡=𝛾𝑏𝑢2𝑣+𝑑2𝑣,(B.2) where the parameters can be identified for the two models and their relationship: 𝑎=𝛼1𝜇𝛾0𝜇,𝑏=𝛼2𝜇𝛾0𝜇,𝑑=𝐷MMP13𝐷VEGF,𝑇=𝐿2𝐷VEGF,𝛾=𝐿2𝐷VEGF𝜇𝐶CH,SVEGF=SMMP13=𝜇𝛾0,(B.3) where 𝑇 is the characteristic time of the biological process and 𝐿 is the characteristic length of the dimensional model. Therefore, defining (𝛾,𝑑,𝑎, and 𝑏) it is possible to obtain the eigenvalues and eigenvectors of the set of (B.3) and from them, the different spatial patterns corresponding to different wave numbers.

In the case of the proposed dimensional model it is necessary to define the nondimensional parameters (𝐿,𝐷VEGF,𝐷MMP13,𝐶CH,𝜇,𝛾,𝛼1, and 𝛼2). To estimate these values, we consider for this work some experimental evidence:(i)The typical concentration of VEGF in human tissue is 𝑆VEGF=𝑆MMP13=1 ng/mL [43].(ii)The study domain is a cubic three-dimensional element of side 𝐿=0,22 mm [44].(iii)Concentration of hypertrophic chondrocytes in the proximal femoral epiphysis is 𝐶CH=65000 cells/mm3 [8].(iv)The coefficient of diffusivity for VEGF is 6,9×105 mm2/s [14].

To reproduce the patterns present in the primary trabecular bone architecture with the proposed model, it is necessary that all parameters are in the Turing space and that they satisfy the restrictions (A.3). Therefore, taking the values for 𝑆VEGF, 𝐿, and 𝐷VEGF, and the set of values for 𝛾, 𝑑, 𝑎, and 𝑏 [8, 22] (see Table 1) that satisfy the restrictions of Turing, and using the relations of (24), we get the set of values to be used for the solution of ((1a), (1b), and (1c)) (see Tables 2(a) and 2(b)).

Table 1: Vibration modes for the nondimensional model of Schakenberg with 𝑎=0.1,𝑏=0.9,𝑇=1, and 𝐿=0.22 mm.
Table 2: (a) Set of values for the process of ossification of the growth plate 𝑎=0.1, b = 0.9, γ = 346,3578, and d = 8,6123. (b) Set of values for the process of ossification of the growth plate 𝑎=0.1, b = 0.9, γ = 700,4675, and d = 8,5756.


This work was supported by the Santo Tomás University and the program of support to research projects (FODEIN-USTA).


  1. D. R. Carter and G. S. Beaupré, Skeletal Function and form Mechanobioloy of Skeletal Development, Aging and Regeneration, Cambridge University Press, New York, NY, USA, 2001.
  2. F. Shapiro, “Developmental bone biology,” in Pediatric Orthopedic Deformities, p. 953, Academic Press, 2002. View at Google Scholar
  3. D. Shier, Hole’s Human Anatomy & Physiology, McGraw-Hill, 2001.
  4. W. E. Roberts and J. K. Hartsfield, “Bone development and function: genetic and environmental mechanisms,” Seminars in Orthodontics, vol. 10, no. 2, pp. 100–122, 2004. View at Publisher · View at Google Scholar · View at Scopus
  5. K. Pelttari, E. Steck, and W. Richter, “The use of mesenchymal stem cells for chondrogenesis,” Injury, vol. 39, supplement 1, pp. S58–S65, 2008. View at Publisher · View at Google Scholar · View at Scopus
  6. N. Ortega, D. J. Behonick, and Z. Werb, “Matrix remodeling during endochondral ossification,” Trends in Cell Biology, vol. 14, no. 2, pp. 86–93, 2004. View at Publisher · View at Google Scholar · View at Scopus
  7. D. A. Garzón-Alvarado, L. M. P. Cortés, and R. P. C. Sandoval, “A mathematical model of epiphyseal development: hypothesis on the cartilage canals growth,” Computer Methods in Biomechanics and Biomedical Engineering, vol. 13, no. 6, pp. 765–772, 2010. View at Publisher · View at Google Scholar · View at Scopus
  8. D. A. Garzón-Alvarado, J. M. García-Aznar, and M. Doblaré, “A reaction-diffusion model for long bones growth,” Biomechanics and Modeling in Mechanobiology, vol. 8, no. 5, pp. 381–395, 2009. View at Publisher · View at Google Scholar · View at Scopus
  9. D. A. Garzón-Alvarado, J. M. García-Aznar, and M. Doblaré, “Appearance and location of secondary ossification centres may be explained by a reaction-diffusion mechanism,” Computers in Biology and Medicine, vol. 39, no. 6, pp. 554–561, 2009. View at Publisher · View at Google Scholar · View at Scopus
  10. E. J. Mackie, Y. A. Ahmed, L. Tatarczuch, K. S. Chen, and M. Mirams, “Endochondral ossification: how cartilage is converted into bone in the developing skeleton,” International Journal of Biochemistry and Cell Biology, vol. 40, no. 1, pp. 46–62, 2008. View at Publisher · View at Google Scholar · View at Scopus
  11. S. Provot and E. Schipani, “Molecular mechanisms of endochondral bone development,” Biochemical and Biophysical Research Communications, vol. 328, no. 3, pp. 658–665, 2005. View at Publisher · View at Google Scholar · View at Scopus
  12. M. B. Goldring, K. Tsuchimochi, and K. Ijiri, “The control of chondrogenesis,” Journal of Cellular Biochemistry, vol. 97, no. 1, pp. 33–44, 2006. View at Publisher · View at Google Scholar · View at Scopus
  13. D. A. Garzón-Alvarado, J. M. García-Aznar, and M. Doblare, “The early bone epiphysis formation: a numerical simulation,” Journal of Biomechanics, vol. 39, pp. S642–S642, 2006. View at Google Scholar
  14. J. E. M. Brouwers, C. C. van Donkelaar, B. G. Sengers, and R. Huiskes, “Can the growth factors PTHrP, Ihh and VEGF, together regulate the development of a long bone?” Journal of Biomechanics, vol. 39, no. 15, pp. 2774–2782, 2006. View at Publisher · View at Google Scholar · View at Scopus
  15. U. Kornak and S. Mundlos, “Genetic disorders of the skeleton: a developmental approach,” American Journal of Human Genetics, vol. 73, no. 3, pp. 447–474, 2003. View at Publisher · View at Google Scholar · View at Scopus
  16. R. Pogue, E. Sebald, L. King, E. Kronstadt, D. Krakow, and D. H. Cohn, “A transcriptional profile of human fetal cartilage,” Matrix Biology, vol. 23, no. 5, pp. 299–307, 2004. View at Publisher · View at Google Scholar · View at Scopus
  17. C. A. Francomano, I. McIntosh, and D. J. Wilkin, “Bone dysplasias in man: molecular insights,” Current Opinion in Genetics and Development, vol. 6, no. 3, pp. 301–308, 1996. View at Publisher · View at Google Scholar · View at Scopus
  18. B. Courtin, A. M. Perault-Staub, and J. F. Staub, “A reaction-diffusion model for trabecular architecture of embryonic periosteal long bone,” Complexity International, vol. 4, pp. 1–17, 1997. View at Google Scholar · View at Scopus
  19. M. J. F. Blumer, S. Longato, and H. Fritsch, “Structure, formation and role of cartilage canals in the developing bone,” Annals of Anatomy, vol. 190, no. 4, pp. 305–315, 2008. View at Publisher · View at Google Scholar · View at Scopus
  20. L. M. Peinado Cortés, J. C. Vanegas Acosta, and D. A. Garzón Alvarado, “A mechanobiological model of epiphysis structures formation,” Journal of Theoretical Biology, vol. 287, no. 12, pp. 13–25, 2011. View at Publisher · View at Google Scholar · View at Scopus
  21. C. C. Van Donkelaar and R. Huiskes, “The PTHrP-Ihh feedback loop in the embryonic growth plate allows PTHrP to control hypertrophy and Ihh to regulate proliferation,” Biomechanics and Modeling in Mechanobiology, vol. 6, no. 1-2, pp. 55–62, 2007. View at Publisher · View at Google Scholar · View at Scopus
  22. D. A. Garzón-Alvarado, M. A. Velasco, and C. A. Narváez-Tovar, “Self-assembled scaffolds using reactiondiffusion systems: a hypothesis for bone regeneration,” Journal of Mechanics in Medicine and Biology, vol. 11, no. 1, pp. 231–272, 2011. View at Publisher · View at Google Scholar · View at Scopus
  23. O. Nilsson, R. Marino, F. De Luca, M. Phillip, and J. Baron, “Endocrine regulation of the growth plate,” Hormone Research, vol. 64, no. 4, pp. 157–165, 2005. View at Publisher · View at Google Scholar · View at Scopus
  24. E. B. Hunziker, J. Wagner, and J. Zapf, “Differential effects of insulin-like growth factor I and growth hormone on developmental stages of rat growth plate chondrocytes in vivo,” Journal of Clinical Investigation, vol. 93, no. 3, pp. 1078–1086, 1994. View at Google Scholar · View at Scopus
  25. R. Randhawa and P. Cohen, “The role of the insulin-like growth factor system in prenatal growth,” Molecular Genetics and Metabolism, vol. 86, no. 1-2, pp. 84–90, 2005. View at Publisher · View at Google Scholar · View at Scopus
  26. W. Schlegel, D. Halbauer, A. Raimann et al., “IGF expression patterns and regulation in growth plate chondrocytes,” Molecular and Cellular Endocrinology, vol. 327, no. 1-2, pp. 65–71, 2010. View at Publisher · View at Google Scholar · View at Scopus
  27. K. Huch, S. Kleffner, J. Stöve, W. Puhl, K. P. Günther, and R. E. Brenner, “PTHrP, PTHr, and FGFR3 are involved in the process of endochondral ossification in human osteophytes,” Histochemistry and Cell Biology, vol. 119, no. 4, pp. 281–287, 2003. View at Google Scholar · View at Scopus
  28. G. V. Segre and K. Lee, “Endochondral bone formation regulation by parathyroid hormone-related peptide, Indian hedgehog, and parathyroid hormone,” in The Parathyroids: Basic and Clinical Concepts, pp. 245–260, Academic Press, San Diego, Calif, USA, 2001. View at Google Scholar
  29. D. A. Stevens and G. R. Williams, “Hormone regulation of chondrocyte differentiation and endochondral bone formation,” Molecular and Cellular Endocrinology, vol. 151, no. 1-2, pp. 195–204, 1999. View at Publisher · View at Google Scholar · View at Scopus
  30. P. De Biase and R. Capanna, “Bone morphogenetic proteins and growth factors: emerging role in regenerative orthopaedic surgery,” Journal of Orthopaedics and Traumatology, vol. 8, no. 1, pp. 43–48, 2007. View at Publisher · View at Google Scholar · View at Scopus
  31. M. F. Carlevaro, S. Cermelli, R. Cancedda, and F. D. Cancedda, “Vascular endothelial growth factor (VEGF) in cartilage neovascularization and chondrocyte differentiation: auto-paracrine role during endochondral bone formation,” Journal of Cell Science, vol. 113, no. 1, pp. 59–69, 2000. View at Google Scholar · View at Scopus
  32. B. C. J. Van Der Eerden, M. Karperien, and J. M. Wit, “Systemic and local regulation of the growth plate,” Endocrine Reviews, vol. 24, no. 6, pp. 782–801, 2003. View at Publisher · View at Google Scholar · View at Scopus
  33. X. Yang and G. Karsenty, “Transcription factors in bone: developmental and pathological aspects,” Trends in Molecular Medicine, vol. 8, no. 7, pp. 340–345, 2002. View at Publisher · View at Google Scholar · View at Scopus
  34. Z. L. Deng, K. A. Sharff, N. Tang et al., “Regulation of osteogenic differentiation during skeletal development,” Frontiers in Bioscience, vol. 13, no. 6, pp. 2001–2021, 2008. View at Publisher · View at Google Scholar · View at Scopus
  35. V. Lefebvre and P. Smits, “Transcriptional control of chondrocyte fate and differentiation,” Birth Defects Research Part C, vol. 75, no. 3, pp. 200–212, 2005. View at Publisher · View at Google Scholar · View at Scopus
  36. L. Troeberg and H. Nagase, “Proteases involved in cartilage matrix degradation in osteoarthritis,” Biochimica et Biophysica Acta, vol. 1824, no. 1, pp. 133–145, 2011. View at Publisher · View at Google Scholar · View at Scopus
  37. S. M. Krane and M. Inada, “Matrix metalloproteinases and bone,” Bone, vol. 43, no. 1, pp. 7–18, 2008. View at Publisher · View at Google Scholar · View at Scopus
  38. G. Murphy and H. Nagase, “Progress in matrix metalloproteinase research,” Molecular Aspects of Medicine, vol. 29, no. 5, pp. 290–308, 2008. View at Publisher · View at Google Scholar · View at Scopus
  39. T. E. Cawston and A. J. Wilson, “Understanding the role of tissue degrading enzymes and their inhibitors in development and disease,” Best Practice and Research: Clinical Rheumatology, vol. 20, no. 5, pp. 983–1002, 2006. View at Publisher · View at Google Scholar · View at Scopus
  40. E. H. Filvaroff, “VEGF and bone,” Journal of Musculoskeletal Neuronal Interactions, vol. 3, no. 4, pp. 304–307, 2003. View at Google Scholar · View at Scopus
  41. M. Murata, K. Yudoh, and K. Masuko, “The potential role of vascular endothelial growth factor (VEGF) in cartilage. How the angiogenic factor could be involved in the pathogenesis of osteoarthritis?” Osteoarthritis and Cartilage, vol. 16, no. 3, pp. 279–286, 2008. View at Publisher · View at Google Scholar · View at Scopus
  42. C. Maes and G. Carmeliet, “Vascular and nonvascular roles of VEGF in bone development,” in VEGF in Development, pp. 79–90, Springer, Austin, Tex, USA, 2008. View at Google Scholar
  43. M. O. Hiltunen, M. Ruuskanen, J. Huuskonen et al., “Adenovirus-mediated VEGF-A gene transfer induces bone formation in vivo.,” The FASEB Journal, vol. 17, no. 9, pp. 1147–1149, 2003. View at Google Scholar · View at Scopus
  44. R. Ruimerman, P. Hilbers, B. Van Rietbergen, and R. Huiskes, “A theoretical framework for strain-related trabecular bone maintenance and adaptation,” Journal of Biomechanics, vol. 38, no. 4, pp. 931–941, 2005. View at Publisher · View at Google Scholar · View at Scopus
  45. M. J. Olszta, X. Cheng, S. S. Jee et al., “Bone structure and formation: a new perspective,” Materials Science and Engineering R, vol. 58, no. 3–5, pp. 77–116, 2007. View at Publisher · View at Google Scholar · View at Scopus
  46. D. A. Garzán-Alvarado and A. M. Ramírez Martinez, “A biochemical hypothesis on the formation of fingerprints using a turing patterns approach,” Theoretical Biology and Medical Modelling, vol. 8, no. 1, article 24, 2011. View at Publisher · View at Google Scholar · View at Scopus
  47. E. J. Cramping and P. K. Maini, “Reaction-diffusion models for biological pattern formation,” Methods and Applications of Analysis, vol. 8, no. 2, pp. 415–428, 2001. View at Google Scholar
  48. A. Gierer and H. Meinhardt, “A theory of biological pattern formation,” Kybernetik, vol. 12, no. 1, pp. 30–39, 1972. View at Publisher · View at Google Scholar · View at Scopus
  49. P. K. Maini, “Mathematical models in morphogenesis,” in Mathematics Inspired by Biology, pp. 151–189, Springer, Berlin, Germany, 1999. View at Google Scholar
  50. P. K. Maini, K. J. Painter, and H. N. P. Chau, “Spatial pattern formation in chemical and biological systems,” Journal of the Chemical Society, vol. 93, no. 20, pp. 3601–3610, 1997. View at Google Scholar · View at Scopus
  51. J. D. Murray, “Pattern formation in integrative biology a marriage of theory and experiment,” Comptes Rendus de l'Academie des Sciences, vol. 323, no. 1, pp. 5–14, 2000. View at Google Scholar · View at Scopus
  52. I. G. Jang and I. Y. Kim, “Computational simulation of trabecular adaptation progress in human proximal femur during growth,” Journal of Biomechanics, vol. 42, no. 5, pp. 573–580, 2009. View at Publisher · View at Google Scholar · View at Scopus
  53. G. A. P. Renders, L. Mulder, G. E. J. Langenbach, L. J. van Ruijven, and T. M. G. J. van Eijden, “Biomechanical effect of mineral heterogeneity in trabecular bone,” Journal of Biomechanics, vol. 41, no. 13, pp. 2793–2798, 2008. View at Publisher · View at Google Scholar · View at Scopus
  54. P. G. Coelho, P. R. Fernandes, H. C. Rodrigues, J. B. Cardoso, and J. M. Guedes, “Numerical modeling of bone tissue adaptation-A hierarchical approach for bone apparent density and trabecular structure,” Journal of Biomechanics, vol. 42, no. 7, pp. 830–837, 2009. View at Publisher · View at Google Scholar · View at Scopus
  55. R. T. Ballock and R. J. O'Keefe, “Current concepts review: the biology of the growth plate,” Journal of Bone and Joint Surgery, vol. 85, no. 4, pp. 715–726, 2003. View at Google Scholar · View at Scopus
  56. D. R. Carter, M. C. H. Van Der Meulen, and G. S. Beaupré, “Mechanical factors in bone growth and development,” Bone, vol. 18, supplement 1, pp. 5S–10S, 1996. View at Google Scholar · View at Scopus
  57. B. De Crombrugghe, V. Lefebvre, R. R. Behringer, W. Bi, S. Murakami, and W. Huang, “Transcriptional mechanisms of chondrocyte differentiation,” Matrix Biology, vol. 19, no. 5, pp. 389–394, 2000. View at Publisher · View at Google Scholar · View at Scopus
  58. B. De Crombrugghe, V. Lefebvre, and K. Nakashima, “Regulatory mechanisms in the pathways of cartilage and bone formation,” Current Opinion in Cell Biology, vol. 13, no. 6, pp. 721–727, 2001. View at Publisher · View at Google Scholar · View at Scopus
  59. T. Pufe, V. Harde, W. Petersen, M. B. Goldring, B. Tillmann, and R. Mentlein, “Vascular endothelial growth factor (VEGF) induces matrix metalloproteinase expression in immortalized chondrocytes,” Journal of Pathology, vol. 202, no. 3, pp. 367–374, 2004. View at Publisher · View at Google Scholar · View at Scopus
  60. A. M. Turing, “The chemical basis of morphogenesis,” Philosophical Transactions of the Royal Society B, vol. 237, no. 641, pp. 37–72, 1952. View at Publisher · View at Google Scholar
  61. A. Madzvamuse, P. K. Maini, and A. J. Wathen, “A moving grid finite element method for the simulation of pattern generation by turing models on growing domains,” Journal of Scientific Computing, vol. 24, no. 2, pp. 247–262, 2005. View at Publisher · View at Google Scholar · View at Scopus