|Home | About | Journals | Submit | Contact Us | Français|
Osteogenetic microenvironment is a complex constitution in which extracellular matrix (ECM) molecules, stem cells and growth factors each interact to direct the coordinate regulation of bone tissue development. Importantly, angiogenesis improvement and revascularization are critical for osteogenesis during bone tissue regeneration processes. In this study, we developed a three-dimensional (3D) multi-scale system model to study cell response to growth factors released from a 3D biodegradable porous calcium phosphate (CaP) scaffold. Our model reconstructed the 3D bone regeneration system and examined the effects of pore size and porosity on bone formation and angiogenesis. The results suggested that scaffold porosity played a more dominant role in affecting bone formation and angiogenesis compared with pore size, while the pore size could be controlled to tailor the growth factor release rate and release fraction. Furthermore, a combination of gradient VEGF with BMP2 and Wnt released from the multi-layer scaffold promoted angiogenesis and bone formation more readily than single growth factors. These results demonstrated that the developed model can be potentially applied to predict vascularized bone regeneration with specific scaffold and growth factors.
Tissue-engineered bone regeneration [1–3] is a complex biological and medical event which includes the scaffolds themselves, cellular and signaling cues, and vascularization, as well as interactions between them . Many in vitro and/or in vivo experimental studies have been conducted to determine strategies promoting osteogenesis and angiogenesis [4, 5]. But due to the complexity of the process of bone regeneration within scaffolds, such coordinated processes involved in bone regeneration have often been examined piecemeal rather than systematically.
Calcium phosphate (CaP) scaffolds (e.g. beta-tricalcium phosphate , hydroxyapatite  and their composites ) are ideal materials for bone repair due to their biocompatibility, adjustable degradation rates, and excellent bioactivity [9, 10]. When used as scaffolds for bone repair, biodegradable CaP scaffolds often contain human primary cells (e.g. mesenchymal stem cells [MSC], osteocytes, and endothelial cells) and growth factors or cytokines to repair bone tissue. Growth factors (e.g. bone morphogenetic protein 2 (BMP2), transforming growth factor β (TGFβ) and Wnt ligands) affect cellular migration and proliferation and osteogenic differentiation of MSC during bone repair. These growth factors can regulate the expression of Runt-related transcription factor 2 (Runx2) and Osterix (Osx) through intracellular proteins or transcription factors including β-Catenin, Smad1/5 and Smad2/3 . Ideally, key growth factors can be programmatically encapsulated and embedded into CaP scaffolds. These cytokines are then released into the microenvironment of the bone graft after being implanted and stimulate the expression of genes responsible for the osteoblastic differentiation from MSC to pre-osteoblasts and then to active osteoblasts  through a variety of signaling pathways [11–14].
Besides osteoblastic differentiation promoted by growth factors released from the porous CaP scaffold, sustained bone formation needs adequate new blood vessel growth to provide nutrients to cells in the interior of the CaP scaffold. Recently, angiogenesis has been a focus of efforts to improve clinical success of bone grafts by increasing osteoblastic cell survival. Among the many growth factors in the bone microenvironment, VEGF and FGF play a critical role in initiating and sustaining vascular growth during bone healing .
Addressing the remaining challenges in the field of bone regeneration requires combining multiple strategies, such as scaffold fabrication, controlled drug release, and vascularization. When combined with associated experimental studies, mathematical and computational modeling potentially provides a systematic rational approach to study bone regeneration. A number of mathematical models of bone regeneration have been developed recently . Geris and co-workers [17–19] proposed a continuum-type model by employing a set of partial differential equations to describe the spatio-temporal evolution of the densities of cells and the concentrations of growth factors, but these models did not include exogenous growth factor release from the biodegradable scaffolds. J.A. Sanz-Herrera and co-workers [20–22] constructed a micro-macro numerical modeling of bone regeneration using a finite element method by integrating two levels: the tissue level or macroscopic scale, and the scaffold pore level or microscopic scale. Checa and co-workers [23–25] developed a mechano-biological model using a lattice approach to simulate cell activity and investigated the effects of cell seeding and mechanical loading on vascularization and tissue formation .
However, none of the above models simulated exogenous growth factor release from the scaffold and cell response to such growth factors. To overcome these limitations, we propose a computational modeling approach that integrates growth factor release from the scaffold and cell response into a vascularized bone regeneration system.
Our previous study  described a multi-scale systems biological model that linked the intracellular and intercellular signaling pathways with dynamic cellular properties to study the combination of effects and optimal doses of cytokines on bone remodeling based on the intracellular signaling pathway. Therein, we mainly focused on responses of cells and bone mass to cytokines (including BMP2, Wnt ligands, and TGFβ) after they were released from the scaffold, but the explicit kinetics of dynamic growth factor release was not studied. Some drug release models [27–30] provide insights into the mass transport and chemical processes involved in growth factor delivery systems, although none of them have considered responses of cell activity and angiogenesis to such drug concentration changes. The present study was designed to develop a computational systems model of bone regeneration within a 3D porous biodegradable CaP scaffold by simulating osteogenesis in response to growth factor release, based on molecular mechanisms and incorporating angiogenesis and nutrient transportation. We first reconstructed the 3D bone regeneration system, and then investigated the effects of pore size and porosity on growth factor release, angiogenesis, and bone formation. Finally, we examined whether the combination of BMP2, Wnt, and VEGF promoted angiogenesis and bone formation better than single growth factors at the same doses.
We developed a 3D multi-scale system model for bone regeneration within porous biodegradable CaP scaffolds. A continuum reaction-diffusion sub-model for exogenous growth factor release and an agent-based sub-model for cell response were coupled in our model. We applied a set of reaction-diffusion equations to model the process of growth factor release from the porous biodegradable CaP scaffold as well as nutrient transportation by the vasculature, and employed an agent-based model  to simulate activities of each osteoblast and endothelial cell. Our model contained four biological/physical scales from micro level to macro level: molecular, cellular, scaffold, and bone tissue scales (Fig. 1). Each scale is described below.
Runx2 and Osx are two crucial transcription factors in osteoblast differentiation and bone regeneration and their expression can be regulated by released growth factors, including BMP2, TGFβ and Wnt, through activation of intracellular proteins or transcription factors such as Smad1/5 (S1), Smad2/3 (S2) and β-Catenin . The molecular regulatory mechanisms involved in the intracellular signaling pathway were modeled using a system of ordinary differential equations as Equations (1.1 – 1.5), please refer to our previous study  for more details.
It was assumed that the binding and/or dissociation reactions in the intracellular signaling pathway are much faster than both cellular phenotype switches and cellular population changes [32, 33], so the quasi-steady state of the intracellular signaling pathway was derived (Supplementary Text S1) in this study to link the short time scale (Minutes / Hours) of the intracellular signaling pathway to the long time scale (Days) of the cellular activity.
MSC and pre-osteoblasts migrate along the gradient of the normalized concentration of growth factors, including BMP2 and TGFβ, and nutrient such as oxygen. The probability () of MSC and pre-osteoblasts migrating along the i-th direction was described in Equation (2). Active osteoblasts were assumed to migrate very little .
where G and N are the concentrations of growth factors and oxygen, respectively, and li is the directional vector along the i-th direction. There are 6 directions for a cell to migrate in our three dimensional lattice.
Activated Runx2 and Osx play different roles in different stages of osteoblastic lineage. Both Runx2 and Osx can promote the differentiation of MSC into pre-osteoblasts, whereas Runx2 can inhibit the differentiation of pre-osteoblasts into active osteoblasts [11, 34]. The probabilities of MSC differentiating into pre-osteoblasts () and pre-osteoblasts differentiating into active osteoblasts () were related to the expression levels of activated Runx2 and Osx. We employed Hill functions to model these situations as Equation (3) and Equation (4).
where and are baseline probabilities of differentiation of MSC into preosteoblasts and differentiation of pre-osteoblasts into active osteoblasts, respectively.
MSC, pre-osteoblasts, and active osteoblasts can proliferate into their daughter cells with different probabilities [23, 35]. Several growth factors, such as IGF1, can regulate the proliferation of osteoblastic cells. In our simulation, the proliferative probabilities (ppro) were set as constants (please refer to Table S2)  without considering the explicit effects of such growth factors.
Hyperbaric oxygen attenuates cell’s apoptosis  and hypoxia makes apoptosis more likely . The probability of apoptosis in osteoblasts was related to the concentration of oxygen (N) as described in Equation (5). Furthermore, if oxygen concentrations are below a threshold ThOxygen, then osteoblasts will die. The equation is
where papop is the baseline probability of the apoptosis of MSC, pre-osteoblasts, or active osteoblasts, Naverage is average oxygen concentration telling hyperbaric oxygen and hypoxia apart, and is a positive constant.
In bone regeneration, growth factors were encapsulated into nanospheres and loaded into porous bio-degradable CaP scaffolds layer by layer. After being implanted into defected bone, calcium phosphate can degrade through hydration reaction and network breakage [38–40].
The diffusion of extracellular liquid such as water and the disintegration of the calcium phosphate network via hydrolysis are described in Equation (6) and Equation (7) using the standard second order rate expression .
where C and M are water concentration and calcium phosphate molecular weight, respectively, DC is water diffusivity, and kC and kM are the degradation rate constants of water and calcium phosphate, respectively. The level-set method  was employed to update the evolving interface of the surface-eroding CaP scaffold. More precisely, at each step in our simulation, if the calcium phosphate molecular weight decreased below a threshold level, then the corresponding bonds of the porous CaP scaffold would degrade.
Meanwhile, BMP2, Wnt ligands, and VEGF were released from the degrading CaP scaffold, and continuously diffused within the scaffold pores. The paracrine or autocrine of these cytokines, such as VEGF, were ignored as we assumed that the concentrations of these cytokines secreted by individual cells were quite low compared to that released from the CaP scaffold. The processes were modeled in the following reaction-diffusion equation (8).
where Gi is the concentration of each growth factor, DGi is the diffusivity of each growth factor (which varies in the regions of calcium phosphate matrix and scaffold pores), Gi,max is the maximum concentration of each growth factor initially loaded into the scaffold, rGi is the release constant, uGMi is the depletion rate of each cytokine by osteoblasts, and dGi is the degradation rate of each cytokine. The dependence of release rate on the concentrations of water and the remaining growth factors in the CaP scaffold was modeled by using a Hill function in the second term according to the Michaelis-Menten Law , where KC is the Michaelis constant.
The time-dependent characteristic function χscaffold(t, x) is equal to 1 in the calcium phosphate matrix; otherwise it is equal to 0 in the pores of the scaffolds. χosteo(t, x) is equal to 1 if a osteoblastic cell is present at x and is equal to 0 elsewhere. χscaffold and χosteo are updated at each simulation step according to the developing profiles of the porous CaP scaffold and osteoblast distribution.
Considering the tiny concentrations of growth factors within scaffold pores at the beginning of their release, the zero initial condition was applied for equation (8). The non-dimensional initial concentration of water was set as 1 in the scaffold pores and 0 in calcium phosphate matrix, while the non-dimensional initial molecular weight of calcium phosphate was set as 1 only in the domain of CaP scaffold. Homogeneous Neumann boundary conditions were imposed for all the above equations by assuming zero flux along the boundary of the considered domain. We solved these equations numerically with the finite difference method .
We assumed that the motion of an individual endothelial cell located at the tip of a capillary sprout governed the motion of the whole sprout, and chemotaxis in response to VEGF gradients guided the motion of the endothelial cells at the capillary sprout tip .
We defined the probability of migration of endothelial cells as
where V is the concentration of VEGF, lk is the directional vector along the k-th direction, α is the chemotactic coefficient, and kV is a positive constant controlling the weight of VEGF concentrations in chemotaxctic sensitivity.
For a given tip endothelial cell at location (i, j, k), the un-normalized migration probabilities P1, P2 … and P6 can be calculated from equation (9). The un-normalized probability P7 (for a tip cell to remain stationary) is the average of P1, P2 … and P6. After normalization, the above probabilities give the likelihood of a given tip endothelial cell to move up, down, right, or left; to the front or the back; or to remain in its current position.
We then normalized the above numbers: , i = 1, 2,…,7 and defined interval I1 = (0, 1], , i = 2, …, 7. For every sprout tip cell, we checked whether the age of vessel exceeded 18 hours and whether there were enough free sites in its nearest neighborhood. If the above conditions were satisfied, two random numbers r1 and r2 between 0 and 1 were generated. If r1 I2 and r2 I3, then we moved two endothelial cells, one down and another to the front of the sprout tip endothelial cell. If the above branching conditions were not satisfied, we generated another random number r between 0 and 1. If r I5, we moved the tip endothelial cell to the right of the sprout tip endothelial cell. If two sprouts encountered each other, a new sprout continued to grow.
Oxygen can be transported by the neo-vasculature to osteoblasts within scaffold pores, which is described in Equation (10).
where N is oxygen concentration, DN is oxygen diffusivity, qN is the vessel permeability for oxygen, Nblood is blood oxygen concentration, and UN is a cell’s uptake rate of oxygen. χves and χoxteo are the characteristic functions as described above. The initial concentration of oxygen was assumed as the average concentration of oxygen within the scaffold. Homogeneous Neumann boundary condition was also set for the above equation by assuming zero flux along the boundary of the considered domain.
The growth rate of bone mass is assumed to be proportional to the number of active osteoblasts, so we used the cumulative number of active osteoblasts as the relative measure of the bone mass at the designated time.
The time-course of bone formation and vascularization within the scaffold was modeled as an iterative process. Here, we summarize our computing algorithm at each step across multi-scales as follows. At the scaffold scale, we solved Equations (6–8) to obtain the spatial distributions of diffusing water concentration, calcium phosphate molecular weight and concentrations of released growth factors. At molecular scale, we used the calculated local BMP2 and Wnt concentrations at each lattice point as the input of the Equation (1) for intracellular signaling pathway for each osteoblast. At the cellular scale, the migration of MSC and pre-osteoblasts was simulated according to Equation (2); meanwhile, the differentiation of MSC and preosteoblasts was regulated by Equation (3), (4), and the proliferation and apoptosis of osteoblastic cells were simulated accordingly. At the bone tissue scale, the spatial concentration distributions of VEGF (Equation 8) will guide the tip endothelial cells' migration and sprout branching (Equation 9). In turn, the remodelled vasculature at tissue scale influences the spatial concentration distributions of oxygen (Equation 10) within CaP scaffold pores.
We performed our simulation on a three-dimensional cubic lattice. The lattice size was set to 50 × 50 × 50, representing a 1.5 mm × 1.5 mm × 1.5 mm portion of the porous biodegradable CaP scaffold. Each lattice point represented a possible space for cells to occupy. The lattice spacing was 30 µm, the approximate diameter of an osteoblast. The endothelial cells are usually smaller, thus we assumed that the neo-vasculature were not particularly closely packed. Each time step of iterative-process simulated the real time length of 1 day. Initially the scaffold was assumed to be loaded with growth factors (e.g. BMP2, Wnt ligands, and VEGF), and the pores of the scaffold were assumed to be filled with granulation tissue and MSCs were seeded inside the construct. The pore sizes (pore diameters) of different scaffolds ranged from 180 µm to 720 µm by intervals of 60 µm for different simulations. The pores in each type of scaffolds were arranged uniformly with equal distance between inter-pores. In addition, two parent blood vessels were initialized on the up and down layers of the lattice.
The program was performed in MATLAB R2007b (MathWorks) on a DELL computer (Intel (R) Core (TM) i5-3470 CPU @ 3.20GHz / 16.0 GB RAM). The average running time of a typical simulation was 23 minutes. Table S1 and S2 list the parameters used in the model. Most were found in the literature  [31, 46, 47], and some parameters involved in the growth factor release module, including degradation rates of water (kC) and molecular weight of scaffold material CaP (kC), and growth factor release constant (rGi) (please refer to Table S2), were determined by fitting them to our experimental data described below.
In the experiment, purified BMP-2 was dissolved with sterile PBS containing 1% bovine serum albumin (BSA) to make a BMP-2 solution, in which the 1% BSA was used to stabilize BMP-2 and inhibit nonspecific adsorption of BMP-2. Each scaffold was immersed in 150 µl BMP-2 solution containing 50 µg of BMP-2 overnight at 4°C. BMP-2 loaded CaP scaffolds were placed into silanized 24-well culture plates. 2 ml of tissue culture medium was added into each well and incubated in a humidified 95% air, 5% CO2 atmosphere at 37°C. 1.2 ml of the tissue culture medium was collected for the BMP elution assay, and then 1.2 ml fresh culture medium was added to the plates to maintain the constant volume (2 ml) of medium at 1, 3, 5, 7 and 14 days. The release was assayed at designated time points. Concentrations of BMP-2 in the removed solution were analyzed using a human BMP-2 ELISA development kit (PeproTEch Inc., Rock Hill, NJ). The degree of BMP-2 release was considered as the BMP-2 concentration of removed solution. Measurements were performed in triplicate for each time point.
First, we reconstructed the 3D bone regeneration system, and then investigated the effects of pore size and porosity of the scaffold on growth factor release, angiogenesis, and bone formation. Finally, we examined the effects of combining growth factors.
First, we reconstructed the system of vascularized bone regeneration within the 3D porous CaP scaffold, which involved the coupled processes of evolution of scaffold degradation, exogenous growth factor release, angiogenesis, MSC differentiation, and cell growth within scaffold pores over time. Fig. 2 shows a typical simulation of our 3D modeling. On day 7, the blood vessels were rare and scattered on the outside surface of the porous CaP scaffold; on days 14 and 21, the newly formed blood vessels grew into the pores located at the periphery of the scaffold. Simultaneously, some blood vessels branched into clusters of tree-branching capillaries. On day 28, a branched vasculature was observed within the peripheral pores of the scaffold. Few blood vessels were predicted to grow into the pores at the center of the scaffold.
The uniformly distributed pores on the scaffolds gradually enlarged due to biodegradation of the scaffold. At early stages, MSC originally recruited from the host filled in the 3D scaffold pores. MSC gradually differentiated into pre-osteoblasts and active osteoblasts, and then active osteoblasts outnumbered MSC and pre-osteoblasts within the scaffold pores.
Fig. 3 shows a 2D cross-sectional view of the spatio-temporal evolution of BMP2, VEGF, and oxygen concentrations; CaP molecular weight; and scaffold/cell profiles near the top layer of the scaffold on days 1, 7, 14, 21, and 28. BMP2 and VEGF were released from the calcium phosphate matrix and continuously diffused in the scaffold pores, and concentrations of these growth factors in the pores of the scaffold increased over time. In our simulation, the maximum concentration of VEGF initially loaded into the biodegradable CaP scaffold was higher in the center of the scaffold and lower in the periphery, so a gradient of VEGF concentration within the pores was produced as it was released from the scaffold. Our simulation also demonstrated that the scaffold with gradient VEGF loading can promote angiogenesis much more than the scaffold with uniform VEGF loading. The oxygen concentration changed along with cells' uptake and the transportation of oxygen in the neo-vasculature. The molecular weight of calcium phosphate decreased from day 1 to day 28 due to hydrolysis.
Fig. 4 shows how numbers of various cells changed over time in the model. At first, MSC, pre-osteoblasts, and active osteoblasts all increased; once they each achieved a peak level, they decreased due to differentiation and/or apoptosis. Numbers of different cell types peaked at different times, consistent with the development of the osteoblast phenotype through osteoblast lineage [48, 49]. The numbers of MSC peaked at around day 3, pre-osteoblasts achieved their peak on day 7, and active osteoblasts peaked on day 14. The number of endothelial cells increased throughout the simulation.
A video showing the dynamic evolution of 3D vascularized bone regeneration within a porous CaP scaffold is included in the Supplementary Materials (Movie S1).
Fig. 5 shows cumulative released BMP2 from scaffolds with a small pore size (480 µm) or a large pore size (720 µm), respectively. The 480 µm pore scaffold has a faster release rate and higher release fraction of BMP2 compared to the 720 µm pore scaffold. This simulation was validated by our experimental data. The mean squared errors of our prediction were 0.0216 for the smaller pore size scaffold and 0.0295 for larger pore size scaffold, respectively. The good agreement between the simulation prediction and experimental data provided an important validation of our model.
Fig. 6 demonstrates the release profiles of BMP2 from scaffolds with pore sizes ranging from 180 µm to 720 µm. Initially, all scaffolds were loaded with the same amount of BMP2. The scaffolds with smaller pore sizes showed a faster release rate and greater amount of BMP2 released, compared to the larger pore sizes. These results suggest that growth factor release rate and release fraction can be tailored, by controlling the pore size of scaffold, to achieve the specific goals of growth factor release or drug delivery in bone tissue engineering.
Scaffold porosity in this study was calculated as the ratio of void pore volume versus total scaffold volume, and the pore size was defined as the diameter of the scaffold pores. As two critical design parameters that can be controlled in the scaffold design, the pore size and porosity were programmed as two independent factors to investigate their effects on bone regeneration. Our simulations have been conducted with different pore sizes and varied values of porosities (see Fig.7a). Fig. 7b shows the number of endothelial cells with different pore sizes after 4 weeks and 8 weeks. There were more endothelial cells when the pore size was 420 µm or 540 µm. Correspondingly, the numbers of active osteoblasts and bone mass were greatest at week 8, when the pore size was 540 µm (Fig. 7c, d). The optimal pore size for both angiogenesis and bone formation, and greatest porosity, was 540 µm.
We also examined the correlations between the scaffold porosity and the number of active osteoblasts, endothelial cells, and bone mass. These all showed trends similar to those of scaffold porosity, and much less correlation with scaffold pore sizes (Fig. 8a). Overall, the numbers of endothelial cells and active osteoblasts, as well as bone mass, increased with scaffold porosity, showing a stronger dependence of angiogenesis and bone mass formation on scaffold porosity than on pore size (Fig. 8b).
The relationship between higher porosity and enhanced bone formation was consistent with previous reports [50–53] and with experimental studies [53–55] that reported no difference in osteogenesis with different pore sizes.
Finally we studied the effects of different growth factors, in single or combined form, released from the scaffold on angiogenesis and bone formation. The total doses of growth factors initially loaded into multi-layer scaffolds were the same. We examined 4 cases: In case a, 10 doses of Wnt were used; in case b, 10 doses of BMP2 were used; in case c, a combination of 4 doses of Wnt, 2 doses of BMP2 and 4 doses of VEGF was used; and in case d, a combination of 2 doses of Wnt, 4 doses of BMP2 and 4 doses of VEGF was used.
Among all four cases, using only Wnt promoted the lowest increase in active osteoblasts (Fig. 9a). When only BMP2 was used, although osteoblasts increased more quickly initially, they dropped markedly after 2 weeks, resulting in a low level of active osteoblasts after 8 weeks. The combination of Wnt, BMP2, and VEGF promoted the most efficient growth of active osteoblasts. Angiogenesis was promoted significantly when exogenous VEGF was introduced (Fig. 9b). After 8 weeks, bone formation was greatest when all 3 growth factors were combined (Fig. 9c). In addition, bone formation was faster with a higher BMP2 dose. Therefore, these results suggest that a high ratio of BMP2 when growth factors are combined may accelerate the progress and the efficiency of bone formation.
Vascularization, important for nutrient transportation and bone implant survival, remains a major challenge in bone tissue engineering. In our simulation, we found that vascularization could be promoted by initially loading more VEGF into the center of the scaffold than in the periphery. Our simulation also demonstrated that neo-vascularization occurred at the periphery of the scaffold. One of the reasons for this was that the scaffold walls obstructed the ingrowth of new capillaries. Therefore wider spaces between the pores of the scaffold, i.e. higher interconnectivity of the pores, may improve angiogenesis in tissue engineered bone regeneration, which could be useful in guiding scaffold design.
Our model successfully captured key kinetic features of growth factor release from the CaP scaffold. We validated the simulation results with experimental data (Fig. 5). Our simulation demonstrated that scaffolds with smaller pore sizes resulted in a faster release rate and higher release of BMP2 compared to scaffolds with bigger pore sizes. This suggests the possibility of tailoring growth factor release rate and release fraction by controlling the pore size of scaffold, to reach the desired goals of drug release or drug delivery in bone tissue engineering.
We interpret the mechanism involved in this phenomenon as follows: diffusional mass transfer is one of the most important processes in drug release and delivery; therefore, diffusion of drug molecules greatly affects drug release rate. Larger pores in a scaffold means there is a larger region in which growth factor molecules can diffuse. Thus, more molecules must be released from the scaffold to fill these larger pore spaces to reach the same concentration level as the smaller pores. As a result, the average concentration of released growth factors within the bigger pores is lower than that within the smaller pores.
This work was designed to model exogenesis growth factor release and its effects on cell behavior within porous scaffold and to reveal optimal strategies for vascularized bone regeneration. Therefore we did not explicitly simulate those factors that were also important in bone regeneration. For instance, it has been shown that mechanical stimuli affects cell proliferation and differentiation as well as angiogenesis involved in bone regeneration , which are also important but beyond the scope of this work. Mechanical stimuli will be integrated into our multi-scale model by employing finite element method [23, 24] in a follow-up study.
The assumption that the pores of the scaffold were uniformly distributed in our model will be substituted by a real scaffold structure where the pores are distributed heterogeneously, which can be achieved by using mico-CT in our future work. Some more experiments will be conducted to estimate the parameters in our model and to validate the predicted results. Some rates or probabilities regulating cell activities including differentiation, proliferation and apoptosis can be estimated from ALP and dsDNA level data. Diffusivities of growth factors and liquids (such as water) in the microenvironment of scaffold pores can be assayed using dye-labeled perfusion MRI. And intravital microscopy will be used for validating spatial distribution of cells and the growth and topology of vasculature.
Our previous study  examined the underlying molecular and cellular mechanisms for the cell response to released cytokines, and predicted that the combination of Wnt and BMP2 can achieve best control of bone remodeling and best bone mass regeneration among the combinations of growth factors. We also explored the optimal dose ratios of given cytokine combinations released from the scaffold to most efficiently control the long-term bone remodeling. In the present study, we show that the combination of VEGF, BMP2, and Wnt could promote osteogenesis more than BMP2 or Wnt alone, consistent with results from previous experiments [57–59].
An ongoing experiment from our lab indicates that in mouse bone marrow stromal cells, adding BMP2 on day 1 followed by IGF1 on day 4 resulted in higher alkaline bone phosphatase concentrations and mineralization levels than other tested candidates of temporal combinations of BMP2 and IGF1. Additionally, published work  indicated that sequential delivery of BMP-2 and IGF-1 using a chitosan gel with gelatin microspheres could enhance early osteoblastic differentiation. These results suggest that sequential delivery of multiple cytokines using a multi-layer biodegradable scaffold can potentially improve vascularized bone regeneration. In our ongoing work, we will test other possible cytokine combinations and optimize drug doses, delivery sequence, and timing.
This study developed a 3D multi-scale model of a bone regeneration system that incorporates several coupled processes, including evolution of scaffold degradation, exogenous growth factor release, osteogenic differentiation and proliferation, angiogenesis, and nutrient transportation over time. Scaffold porosity played a more dominant role in promoting bone formation and angiogenesis compared to pore size, while the pore size could be controlled to tailor the release rate and release fraction of the exogenous growth factors. These predictions from our model were confirmed with our experimental data. Furthermore, with the same total doses, the combination of gradient VEGF plus BMP2 and Wnt released from the multi-layer scaffold promoted angiogenesis and bone formation to a much greater degree than BMP2 or Wnt alone. These results suggest that designing pore size and porosity of the scaffold, as well as combining growth factors, are critical for enhancing vascularized bone regeneration.
We would like to acknowledge the members of Center for Bioinformatics and Systems Biology at Wake Forest University School of Medicine, particularly Dr. Jing Su and Dr. Ruoying Chen, for valuable discussions. This work was supported by NIH R01LM010185-03 (Zhou), NIH U01HL111560-01 (Zhou), NIH 1R01DE022676-01 (Zhou), U01 CA166886-01 (Zhou), DOD-W81XWH-11-2-0168-P4 (Zhou), NIH R01AR057837 (Yang), NIH R01DE021468 (Yang) and DOD W81XWH-10-1-0966 (Yang).
Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.