Spongiosa Primary Development: A Biochemical Hypothesis by Turing Patterns Formations

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:𝜕𝑆_VEGF𝜕𝑡=𝐶_CH𝛼₁−𝜇𝑆_VEGF+𝛾₀𝑆²_VEGF𝑆_MMP13+𝐷_VEGF∇²𝑆_VEGF,(1a)𝜕𝑆_MMP13𝜕𝑡=𝐶_CH𝛼₂−𝛾₀𝑆²_VEGF𝑆_MMP13+𝐷_MMP13∇²𝑆_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. 𝛼₁ 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 𝐷_VEGF and 𝐷_MMP13 are the diffusion coefficients of VEGF and MMP13, respectively. In the biological interpretation of the above equations the term 𝛾₀𝑆²_VEGF𝑆_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𝛼₁−𝜇𝑆_VEGF+𝛾₀𝑆²_VEGF𝑆_MMP13−𝐷_VEGF∇²𝑆_VEGF=0,𝜕𝑆_MMP13𝜕𝑡−𝐶_CH𝛼₂−𝛾₀𝑆²_VEGF𝑆_MMP13−𝐷_MMP13∇²𝑆_MMP13=0,𝜕𝑐_Bone𝜕𝑡−𝜂𝑆^𝑛_VEGF𝑆^𝑛_VEGF+𝑆^𝑛_umbral⋅𝑇^𝑟_𝑎𝑇^𝑟_𝑎+𝑡^𝑟=0.(2)

Using weighted residuals, the system of (2) takes the form _Ω𝑤₁𝜕𝑆_VEGF𝜕𝑡−𝐶_CH𝛼₁−𝜇𝑆_VEGF+𝛾₀𝑆²_VEGF𝑆_MMP13−𝐷_VEGF∇²𝑆_VEGF𝑑Ω=0,_Ω𝑤₂𝜕𝑆_MMP13𝜕𝑡−𝐶_CH𝛼₂−𝛾₀𝑆²_VEGF𝑆_MMP13−𝐷_MMP13∇²𝑆_MMP13𝑑Ω=0,_Ω𝑤₃𝜕𝑐_Bone𝜕𝑡−𝜂𝑆^𝑛_VEGF𝑆^𝑛_VEGF+𝑆^𝑛_umbral⋅𝑇^𝑟_𝑎𝑇^𝑟_𝑎+𝑡^𝑟𝑑Ω=0,(3) 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 𝑟_𝑆VEGF=_Ω𝑤₁𝜕𝑆_VEGF𝜕𝑡𝑑Ω−_Ω𝑤₁𝐶_CH𝛼₁𝑑Ω+_Ω𝑤₁𝐶_CH𝜇𝑆_VEGF𝑑Ω−_Ω𝑤₁𝐶_CH𝛾₀𝑆²_VEGF𝑆_MMP13𝑑Ω+_Ω𝐷_VEGF∇𝑤₁∇𝑆_VEGF𝑑Ω−_Γ𝑤₁𝐷_VEGF∇𝑆_VEGF⋅𝑛𝑑Γ=0,𝑟_SMMP13=_Ω𝑤₂𝜕𝑆_MMP13𝜕𝑡𝑑Ω−_Ω𝑤₂𝐶_CH𝛼₂𝑑Ω+_Ω𝑤₂𝐶_CH𝛾₀𝑆²_VEGF𝑆_MMP13𝑑Ω+_Ω𝐷_MMP13∇𝑤₂∇𝑆_MMP13𝑑Ω−_Γ𝑤₂𝐷_MMP13∇𝑆_MMP13⋅𝑛𝑑Γ=0,𝑟_𝐶Bone=_Ω𝑤₃𝑑𝑐_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=_Ω𝑤₁𝜕𝑆_VEGF𝜕𝑡𝑑Ω−_Ω𝑤₁𝐶_CH𝛼₁𝑑Ω+_Ω𝑤₁𝐶_CH𝜇𝑆_VEGF𝑑Ω−_Ω𝑤₁𝐶_CH𝛾₀𝑆²_VEGF𝑆_MMP13𝑑Ω+_Ω𝐷_VEGF∇𝑤₁∇𝑆_VEGF𝑑Ω=0,𝑟_𝑆MMP13=_Ω𝑤₂𝜕𝑆_MMP13𝜕𝑡𝑑Ω−_Ω𝑤₂𝐶_CH𝛼₂𝑑Ω+_Ω𝑤₂𝐶_CH𝛾₀𝑆²_VEGF𝑆_MMP13𝑑Ω+_Ω𝐷_MMP13∇𝑤₂∇𝑆_MMP13𝑑Ω=0,𝑟_𝐶Bone=_Ω𝑤₃𝜕𝑐_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𝛼₁𝑑Ω+_Ω𝐍^𝑇𝐶_CH𝜇𝑆^𝑒_VEGF𝑑Ω−_Ω𝐍^𝑇𝐶_CH𝛾₀𝑆^𝑒_VEGF²𝑆^𝑒_MMP13𝑑Ω+_Ω𝐷_VEGF∇𝐍^𝑇∇𝑆^𝑒_VEGF𝑑Ω,𝑟^𝑒_𝑆MMP13=_Ω𝐍^𝑇𝜕𝑆^𝑒_MMP13𝜕𝑡𝑑Ω−_Ω𝐍^𝑇𝐶_CH𝛼₂𝑑Ω+𝐍^𝑇+_Ω𝑤₂𝐶_CH𝛾₀𝑆^𝑒_VEGF²𝑆^𝑒_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𝛼₁𝑑Ω+∫_Ω𝐍^𝑇𝐶_CH𝜇𝑆^𝑒_VEGF𝑑Ω−∫_Ω𝐍^𝑇𝐶_CH𝛾₀𝑆^𝑒_VEGF²𝑆^𝑒_MMP13𝑑Ω⎞⎟⎟⎠^𝑡+Δ𝑡+(1−𝛼)⎛⎜⎜⎝∫_Ω𝐷_VEGF∇𝐍^𝑇∇𝑆^𝑒_VEGF𝑑Ω−∫_Ω𝐍^𝑇𝐶_CH𝛼₁𝑑Ω+∫_Ω𝐍^𝑇𝐶_CH𝜇𝑆^𝑒_VEGF𝑑Ω−∫_Ω𝐍^𝑇𝐶_CH𝛾₀𝑆^𝑒_VEGF²𝑆^𝑒_MMP13𝑑Ω⎞⎟⎟⎠^𝑡,𝑟^𝑒_𝑆MMP13=_Ω𝐍^𝑇𝑆^𝑒_MMP13−𝑆^𝑒_MMP13Δ𝑡𝑑Ω+𝛼⎛⎜⎜⎝∫_Ω𝐷_MMP13∇𝐍^𝑇∇𝑆^𝑒_MMP13𝑑Ω−∫_Ω𝐍^𝑇𝐶_CH𝛼₂𝑑Ω+𝐍^𝑇+∫_Ω𝑤₂𝐶_CH𝛾₀𝑆^𝑒_VEGF²𝑆^𝑒_MMP13𝑑Ω⎞⎟⎟⎠^𝑡+Δ𝑡+(1−𝛼)⎛⎜⎜⎝∫_Ω𝐷_MMP13∇𝐍^𝑇∇𝑆^𝑒_MMP13𝑑Ω−∫_Ω𝐍^𝑇𝐶_CH𝛼₂𝑑Ω+𝐍^𝑇+∫_Ω𝑤₂𝐶_CH𝛾₀𝑆^𝑒_VEGF²𝑆^𝑒_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𝛾₀2𝑆^𝑒_VEGF𝑆^𝑒_MMP13𝑑Ω⎞⎟⎟⎠,𝜕𝑟^𝑒_𝑆VEGF𝜕𝑆^𝑒_MMP13=−𝛼_Ω𝐍^𝑇𝐍𝐶_CH𝛾₀𝑆^𝑒_VEGF²𝑑Ω,𝜕𝑟^𝑒_𝑆MMP13𝜕𝑆^𝑒_VEGF=𝛼_Ω𝐍^𝑇𝐍𝐶_CH𝛾₀2𝑆^𝑒_VEGF𝑆^𝑒_MMP13𝑑Ω,𝜕𝑟^𝑒_𝑆MMP13𝜕𝑆^𝑒_MMP13=1Δ𝑡_Ω𝐍^𝑇𝐍𝑑Ω+𝛼_Ω𝐷_MMP13∇𝐍^𝑇∇𝐍𝑑Ω+_Ω𝐍^𝑇𝐍𝐶_CH𝛾₀𝑆^𝑒_VEGF²𝑑Ω,𝜕𝑟^𝑒_𝐶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/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 (𝑆_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, 47–51].

(a)

(b)

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.

(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: 𝜕𝑢_𝑖𝜕𝑡=𝐷_𝑖∇²𝑢_𝑖+𝑓_𝑖𝑢₁,𝑢₂,…,𝑢_𝑚,𝑥∈Ω,𝑡≥0,(A.1) 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 𝑓(𝑢,𝑣)=0 and 𝑔(𝑢,𝑣)=0 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 𝐉≡⎛⎜⎜⎝𝑓_𝑢𝑓_𝑣𝑔_𝑢𝑔_𝑣⎞⎟⎟⎠=⎛⎜⎜⎝++−−⎞⎟⎟_{TipoI}⎠𝑜⎛⎜⎜⎝+−+−⎞⎟⎟_{TipoII}⎠.(A.2)

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: 𝑓_𝑢+𝑔_𝑣<0,𝑓_𝑢𝑔_𝑣−𝑓_𝑣𝑔_𝑢>0,𝑑𝑓_𝑢+𝑔_𝑣>0,𝑑𝑓_𝑢+𝑔_𝑣²≥4𝑑𝑓_𝑢𝑔_𝑣−𝑓_𝑣𝑔_𝑢,(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 𝜕𝑢𝜕𝑡=𝐷₁∇²𝑢+𝑓(𝑢,𝑣),𝜕𝑣𝜕𝑡=𝐷₂∇²𝑣+𝑔(𝑢,𝑣),(A.4) where we have the following terms of reaction: 𝑓(𝑢,𝑣)=𝑎−𝑢+𝑢²𝑣,𝑔(𝑢,𝑣)=𝑏−𝑢²𝑣.(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𝑢𝑣,𝜕𝑔𝜕𝑣=−𝑢².(A.8)

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

Replacing (A.9) into (A.8), we have 𝜕𝑓𝜕𝑢=2𝑏𝑎+𝑏−1,𝜕𝑓𝜕𝑣=(𝑎+𝑏)²,𝜕𝑔𝜕𝑢=−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,𝐷_Ihh≠0), 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)=((𝛼₁+𝛼₂)/𝜇,𝛼₂𝜇²/𝛾₀(𝛼₁+𝛼₂)²). 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: 𝐶_CH2𝛾₀𝑆^∗_VEGF𝑆^∗_MMP13−𝛾₀𝑆^∗_VEGF−𝜇<0,𝐶²_CH𝛾₀𝑆^∗2_VEGF𝜇−2𝛾₀𝑆^∗_VEGF𝑆^∗_MMP13+2𝛾₀²𝑆^∗3_VEGF𝑆^∗_MMP13>0,𝐶_CH𝐷_MMP132𝛾₀𝑆^∗_VEGF𝑆^∗_MMP13−𝜇−𝐷_VEGF𝛾₀𝑆^∗2_VEGF>0,𝐶²_CH𝐷_MMP132𝛾₀𝑆^∗_VEGF𝑆^∗_MMP13−𝜇−𝐷_VEGF𝛾₀𝑆^∗2_VEGF²⋯−4𝐷_VEGF𝐷_MMP13𝐶²_CH𝛾₀𝑆^∗2_VEGF×𝜇−2𝛾₀𝑆^∗_VEGF𝑆^∗_MMP13+2𝛾₀²𝑆^∗3_VEGF𝑆^∗_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 𝜕𝑢𝜕𝑡=𝛾𝑎−𝑢+𝑢²𝑣+∇²𝑢,𝜕𝑣𝜕𝑡=𝛾𝑏−𝑢²𝑣+𝑑∇²𝑣,(B.2) where the parameters can be identified for the two models and their relationship: 𝑎=𝛼₁𝜇𝛾₀𝜇,𝑏=𝛼₂𝜇𝛾₀𝜇,𝑑=𝐷_MMP13𝐷_VEGF,𝑇=𝐿²𝐷_VEGF,𝛾=𝐿²𝐷_VEGF𝜇𝐶_CH,S_VEGF=S_MMP13=𝜇𝛾₀,(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,𝜇,𝛾,𝛼₁, and 𝛼₂). 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/mm³ [8].(iv)The coefficient of diffusivity for VEGF is 6,9×10⁻⁵ 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 𝑆_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.

Acknowledgment

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