Computational and Mathematical Methods in Medicine

Volume 2012 (2012), Article ID 748302, 14 pages

doi:10.1155/2012/748302

Research Article

Spongiosa Primary Development: A Biochemical Hypothesis by Turing Patterns Formations

Oscar Rodrigo López-Vaca¹ and Diego Alexander Garzón-Alvarado²

¹Grupo de Investigación en Estudios y Aplicaciones de Ingeniería Mecánica (GEAMEC), Universidad Santo Tomás, Bogotá, Colombia

²Grupo 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.

Abstract.

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 [1–4]. 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, [7–9]. 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, 11–13] 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 [15–17].

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 [7–9, 13, 14, 18–21]. 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. [7–9, 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 [36–38]. 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:

where is the concentration of hypertrophic chondrocytes, , represent the concentrations of VEGF and MMP13, respectively. The remaining terms are model parameters. and 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, regulates the nonlinear interaction between the concentration of MMP13-VEGF quantifying the concentration or molecular inhibition of each molecular factor, and and are the diffusion coefficients of VEGF and MMP13, respectively. In the biological interpretation of the above equations the term represents the nonlinear activation of (production of VEGF by the presence of MMP13) and the nonlinear consumption of (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 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, 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 () and MMP13 (). The time integration was performed using the trapezoidal rule.

3.3.1. Weak Formulation

Let ((1a), (1b), and (1c)) be rewritten as

Using weighted residuals, the system of (2) takes the form

where represents the domain of the problem which is limited by the boundary . , , and are weight functions. Using Green's theorem and weakening the system of (3), the residue of the problem is given by

Defining null flow conditions on the contour of the problem, the border-related terms in (4) are canceled, so the residue is expressed as

To discretize the finite element solution we use approaching functions written as a linear combination of orthogonal functions as shown in the following:

where , , and represent the shape functions which depend only on the space used for the formulation, and are the values of and 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 , and are residue vectors for each equation and is the gradient vector of the approximation functions

Using the time discretization equation based on the Crank-Nicolson equations, the equations in (7) are transformed into the following:

where and are the nodal values and 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:

The nodal values and in time can be approximated by the iterative algorithm of Newton-Raphson, as described in (10). In this equation, and represent the difference of the nodal values in two consecutive iterations

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 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/mm³. 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  [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 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 () and 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, 47–51].

(a)

(b)

Figure 4: Temporal evolution of the concentration of 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.

(a)

(b)

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, 55–58]. 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.

Appendices

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:

where the unknown functions can be interpreted as a reagent of concentration, the term 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 and in the plane , for a system of two chemical components, are cut in the point , 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

It is possible to find a Turing instability if the reaction rates and diffusion coefficients allow the fixed point 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:

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

where we have the following terms of reaction:

For the point of equilibrium, we have

The Jacobian (stabilization matrix) is given by

Therefore,

For a steady point , we have and ; therefore, and are given by

Replacing (A.9) into (A.8), we have

Taking and [9] and replacing in (A.10), we have

Therefore, the stability matrix takes the following structure of signs:

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 , 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 and , and leads to . 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 , 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:

These inequalities define a domain in the parameter space, known as the Turing space, where the uniform steady state 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 , respectively, through , we can obtain

where the parameters can be identified for the two models and their relationship:

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 (, and ). To estimate these values, we consider for this work some experimental evidence:

(i)The typical concentration of VEGF in human tissue is  ng/mL [43].(ii)The study domain is a cubic three-dimensional element of side  mm [44].(iii)Concentration of hypertrophic chondrocytes in the proximal femoral epiphysis is  cells/mm³ [8].(iv)The coefficient of diffusivity for VEGF is  mm²/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 , , and , 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 , and  mm.

γ

Wave mode

8,6123
346,3578
0,1
0,9
(2,2,2)

8.5736
700,4675
0,1
0,9
(4,2,2)

Table 2: (a) Set of values for the process of ossification of the growth plate , b = 0.9, γ = 346,3578, and d = 8,6123. (b) Set of values for the process of ossification of the growth plate , b = 0.9, γ = 700,4675, and d = 8,5756.

(a)

Parameter
Value
Units

0,6562
mm³/cell day

6,56 (10¹¹)
mm⁹/cell day g²

6,56 ()
g/cell day

5,91 ()
g/cell day

5,9 ()
mm²/s

11,7
Minutes

(b)

Parameter
Value
Units

1,327
mm³/cell day

1,327 (10¹²)
mm⁹/cell day g²

1,327 ()
g/cell day

1,19 ()
g/cell day

5,9 ()
mm²/s

11,7
Minutes

Acknowledgment

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

References

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.

F. Shapiro, “Developmental bone biology,” in Pediatric Orthopedic Deformities, p. 953, Academic Press, 2002.

D. Shier, Hole’s Human Anatomy & Physiology, McGraw-Hill, 2001.

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.

K. Pelttari, E. Steck, and W. Richter, “The use of mesenchymal stem cells for chondrogenesis,” Injury, vol. 39, supplement 1, pp. S58–S65, 2008.

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.

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.

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.

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.

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.

S. Provot and E. Schipani, “Molecular mechanisms of endochondral bone development,” Biochemical and Biophysical Research Communications, vol. 328, no. 3, pp. 658–665, 2005.

M. B. Goldring, K. Tsuchimochi, and K. Ijiri, “The control of chondrogenesis,” Journal of Cellular Biochemistry, vol. 97, no. 1, pp. 33–44, 2006.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

S. M. Krane and M. Inada, “Matrix metalloproteinases and bone,” Bone, vol. 43, no. 1, pp. 7–18, 2008.

G. Murphy and H. Nagase, “Progress in matrix metalloproteinase research,” Molecular Aspects of Medicine, vol. 29, no. 5, pp. 290–308, 2008.

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.

E. H. Filvaroff, “VEGF and bone,” Journal of Musculoskeletal Neuronal Interactions, vol. 3, no. 4, pp. 304–307, 2003.

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.

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.

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.

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.

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.

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.

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.

A. Gierer and H. Meinhardt, “A theory of biological pattern formation,” Kybernetik, vol. 12, no. 1, pp. 30–39, 1972.

P. K. Maini, “Mathematical models in morphogenesis,” in Mathematics Inspired by Biology, pp. 151–189, Springer, Berlin, Germany, 1999.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

A. M. Turing, “The chemical basis of morphogenesis,” Philosophical Transactions of the Royal Society B, vol. 237, no. 641, pp. 37–72, 1952.

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.