Mechanical–chemical coupled modeling of bone regeneration within a biodegradable polymer scaffold loaded with VEGF

Abstract

Repairing critical-size bone defects with engineered scaffolds remains a challenge in orthopedic practice. Insufficient vascularization is a major reason causing the failure of bone regeneration within scaffolds. Loading exogenous vascular endothelial growth factor (VEGF) in biodegradable polymer scaffolds and controlling its release rate can promote vascularization in scaffolds and accelerate bone regeneration during bone repair. In this study, we developed a 3D mechanical–chemical model of bone regeneration, which combines multiple mechanical–chemical factors including mechanical stimulation, scaffold degradation, VEGF release and transportation, vascularization and oxygen delivery. This model simulated the coupled dynamic mechanical–chemical environments during bone regeneration and scaffold degradation and predicted bone growth under different mechanical–chemical conditions. Moreover, the predictive power of the model was preliminarily validated by experimental data in literature. Based on the validated model, the effect of exogenous VEGF doses on bone regeneration and the optimal doses under different mechanical stimulations was investigated. The simulation results suggested that there was an optimal range of VEGF doses, which promoted the efficiency of bone regeneration, and an appropriate mechanical stimulation improved the effect of VEGF on bone regeneration. The present work may provide a useful platform for future design of bone scaffolds to regenerate functional bones.

Appendices

Appendix 1: Numerical implementation of sprouting angiogenesis

We used Euler finite difference approximations to discretize Eq. (14), and the coefficients P₀–P₆ in Eq. (15) have the forms,

$$\begin{aligned} P_{0} = & 1 - \frac{{6D_{{{\text{EC}}}} T_{{{\text{vessel}}}} }}{{l_{{{\text{el}}}}^{2} }} \\ & \quad + \frac{{\theta_{0} T_{{{\text{vessel}}}} }}{{4k_{1} l_{{{\text{el}}}}^{2} \left( {1 + \frac{{C_{{{\text{V}}\,l,m,n}}^{q} }}{{k_{1} }}} \right)^{2} }}[(C_{{{\text{V}}\,l + 1,m,n}}^{q} - C_{{{\text{V}}\,l - 1,m,n}}^{q} )^{2} + (C_{{{\text{V}}\,l,m + 1,n}}^{q} - C_{{{\text{V}}\,l,m - 1,n}}^{q} )^{2} + (C_{{{\text{V}}\,l,m,n + 1}}^{q} - C_{{{\text{V}}\,l,m,n - 1}}^{q} )^{2} ] \\ & \quad - \frac{{\theta_{0} k_{1} T_{{{\text{vessel}}}} }}{{l_{{{\text{el}}}}^{2} (k_{1} + C_{{{\text{V}}\,l,m,n}}^{q} )}}(C_{{{\text{V}}\,l + 1,m,n}}^{q} + C_{{{\text{V}}\,l - 1,m,n}}^{q} + C_{{{\text{V}}\,l,m + 1,n}}^{q} + C_{{{\text{V}}\,l,m - 1,n}}^{q} + C_{{{\text{V}}\,l,m,n + 1}}^{q} + C_{{{\text{V}}\,l,m,n - 1}}^{q} - 6C_{{{\text{V}}\,l,m,n}}^{q} ) \\ P_{1} = & \frac{{D_{{{\text{EC}}}} T_{{{\text{vessel}}}} }}{{l_{{{\text{el}}}}^{2} }} - \frac{{\theta_{0} k_{1} T_{{{\text{vessel}}}} }}{{4l_{el}^{2} (k_{1} + C_{{{\text{V}}\,l,m,n}}^{q} )}}(C_{{{\text{V}}\,l + 1,m,n}}^{q} - C_{{{\text{V}}\,l - 1,m,n}}^{q} ) \\ P_{2} = & \frac{{D_{{{\text{EC}}}} T_{{{\text{vessel}}}} }}{{l_{{{\text{el}}}}^{2} }} + \frac{{\theta_{0} k_{1} T_{{{\text{vessel}}}} }}{{4l_{{{\text{el}}}}^{2} (k_{1} + C_{{{\text{V}}\,l,m,n}}^{q} )}}(C_{{{\text{V}}\,l + 1,m,n}}^{q} - C_{{{\text{V}}\,l - 1,m,n}}^{q} ) \\ P_{3} = & \frac{{D_{{{\text{EC}}}} T_{{{\text{vessel}}}} }}{{l_{{{\text{el}}}}^{2} }} - \frac{{\theta_{0} k_{1} T_{{{\text{vessel}}}} }}{{4l_{{{\text{el}}}}^{2} (k_{1} + C_{{{\text{V}}\,l,m,n}}^{q} )}}(C_{{{\text{V}}\,l,m + 1,n}}^{q} - C_{{{\text{V}}\,l,m - 1,n}}^{q} ) \\ P_{4} = & \frac{{D_{{{\text{EC}}}} T_{{{\text{vessel}}}} }}{{l_{{{\text{el}}}}^{2} }} + \frac{{\theta_{0} k_{1} T_{{{\text{vessel}}}} }}{{4l_{{{\text{el}}}}^{2} (k_{1} + C_{{{\text{V}}\,l,m,n}}^{q} )}}(C_{{{\text{V}}\,l,m + 1,n}}^{q} - C_{{{\text{V}}\,l,m - 1,n}}^{q} ) \\ P_{5} = & \frac{{D_{{{\text{EC}}}} T_{{{\text{vessel}}}} }}{{l_{{{\text{el}}}}^{2} }} - \frac{{\theta_{0} k_{1} T_{{{\text{vessel}}}} }}{{4l_{{{\text{el}}}}^{2} (k_{1} + C_{{{\text{V}}\,l,m,n}}^{q} )}}(C_{{{\text{V}}\,l,m,n + 1}}^{q} - C_{{{\text{V}}\,l,m,n - 1}}^{q} ) \\ P_{6} = & \frac{{D_{{{\text{EC}}}} T_{{{\text{vessel}}}} }}{{l_{{{\text{el}}}}^{2} }} + \frac{{\theta_{0} k_{1} T_{{{\text{vessel}}}} }}{{4l_{{{\text{el}}}}^{2} (k_{1} + C_{{{\text{V}}\,l,m,n}}^{q} )}}(C_{{{\text{V}}\,l,m,n + 1}}^{q} - C_{{{\text{V}}\,l,m,n - 1}}^{q} ). \\ \end{aligned}$$

The subscripts l, m, n specify the location on the grid and the superscript q specifies the time steps. The schematic diagram of the 3D7P angiogenesis model is shown in Fig.

Fig. 12

Schematic diagram of the 3D7P angiogenesis model

12. Acceding to the average elongation rate for individual capillaries observed with 5 μm/h (Utzinger et al. 2015), it takes about 5 h for a new blood vessel to extend from one grid point to its neighboring point; we therefore took the time increment $T_{{{\text{vessel}}}}$ = 0.25 days.

The rules for branching and anastomosis (Anderson and Chaplain 1998) are as follows. The generation of a new sprout (branching) was assumed to occur only from existing sprout tips, and the newly formed sprouts cannot branch immediately. Sufficient local space is also requisite for the formation of a new sprout. Given that the above conditions are satisfied, it was assumed that each sprout tip has a probability of generating a new sprout and this probability is dependent on the local VEGF concentration. Anastomosis was assumed to occur when one sprout tip encounters another sprout or blood vessels. Only one of the original sprouts continues to grow as a result of the tip-to-tip fusions.

Appendix 2: Numerical implementation of VEGF and oxygen transportation

The partial differential equations Eqs. (13) and (16) are discretized by using explicit finite difference methods. The space length was aligned with the length of finite elements l_el = 25 µm. The time increment of VEGF $T_{{\text{V}}}$ was chosen with 3.5 s to satisfy the stability condition $\left| {\frac{{T_{{\text{V}}} }}{{l_{{{\text{el}}}}^{2} }}D_{{\text{V}}}^{{{\text{eff}}}} } \right| \le \frac{1}{6}$. In terms of oxygen, the oxygen concentration was assumed to be in a quasi-steady state during each time increment of the angiogenesis, since the transportation of oxygen is much faster than the characteristic time for cell proliferation and migration, and the number of elements with relative bone density $\stackrel{-}{{\rho }_{\mathrm{b}}}\left(t\right)$>0 changes very little during a time increment of angiogenesis.

Considering the tiny concentration of VEGF in the interstitial tissue at the beginning of its release, we took zero as the initial condition for Eq. (13). As the measurements of bone marrow aspirates from healthy donors had mean oxygen tension values of 54.9 mmHg (Harrison et al. 2002), we therefore took 54.9 mmHg (0.073 nmol mm⁻³) as the initial condition for Eq. (16).

A homogenous Neumann boundary condition was imposed on the left, right, front and back boundaries of the computational domain by assuming zero flux of oxygen and VEGF along these boundaries, since the initial scaffold geometry and the mechanical load were both symmetric with respect to the plane x = 0.5 mm and plane y = 0.5 mm. A Dirichlet boundary condition of oxygen concentration was imposed on the top and bottom boundaries by assuming that the oxygen concentration is always the initial value (0.073 nmol mm⁻³) on these boundaries. Since the initial concentration of VEGF is extremely low in the surrounding environment, VEGF flux outflow will significantly change the concentration in the surrounding environment. Therefore, a dynamic Dirichlet boundary condition of VEGF concentration was imposed on the top and bottom boundaries, as shown in Fig.

Fig. 13

The boundary concentration of VEGF on the top and bottom boundaries

13. C_V^s is the boundary concentration of VEGF imposed on the top or bottom boundaries, which indicates the boundary concentration on the side of surrounding environment. C_Vⁱ is the boundary concentration on the side of implant. The concentration of VEGF in the surrounding environment was assumed to be linear with the distance from the implant, and the further away from the implant, the lower the concentration. Thus, C_V^s can be given by $C_{{\text{V}}}^{{\text{s}}} = \frac{{L_{0} - l_{{{\text{el}}}} }}{{L_{0} }}C_{{\text{V}}}^{i}$. L₀ represents the distance that VEGF diffuses from the implant after scaffold implantation, estimated by $L_{0} = (D_{{\text{V}}}^{{{\text{bulk}}}} t)^{\frac{1}{2}}$, where t refers to the time after implantation and $D_{{\text{V}}}^{{{\text{bulk}}}}$ represents the VEGF diffusion coefficient in the surrounding environment.

Appendix 3: Parameter values of $C_{{{\text{O}}_{2} ,{\text{upper}}}}$, $C_{{{\text{O}}_{2} ,{\text{lower}}}}$, $u_{{{\text{O}}_{2} }}$ and $q_{{{\text{O}}_{2} }}$

In vitro experiments found that the osteoblast mineralization was reduced by 60–70% when the oxygen tension was decreased from normoxia to 2% (about 15 mmHg, 0.02 nmol mm⁻³) (Nicolaije et al. 2012), and the formation of mineralized bone nodule was almost abolished when the oxygen tension was reduced further to 0.2% (about 1.5 mmHg, 0.002 nmol mm⁻³) (Utting et al. 2006). Therefore, we took $C_{{{\text{O}}_{2} ,{\text{lower}}}}$ = 1.5 mmHg and ξ (15 mmHg) = 0.35 $u_{\max }$. $C_{{{\text{O}}_{2} ,{\text{upper}}}}$ was obtained with the linear interpolation method to be 40 mmHg (about 5.3%, 0.053 nmol mm⁻³).

The measurements of oxygen consumption rate by osteoblasts cultured in static T-flasks and encapsulated mediums were 5.56 × 10^–6 and 1.25 × 10^–7 µmol min⁻¹cell⁻¹, respectively (Wang et al. 2011). Hence, we obtained 3.6 × 10^–3 µmol mm⁻³ s⁻¹ element⁻¹ as the estimate of $u_{{{\text{O}}_{2} }}$, which was validated in Sect. 3.2.

We assumed that the blood vessels passed through the ISF elements from the center, as shown in Fig.

Fig. 14

Schematic diagram of an ISF element through which a blood vessel passes (the blue zone is the ISF and the red zone is the blood vessel)

14. Therefore, the second part of Eq. (16), the free oxygen flux across the blood vessel wall, satisfies an equation of the form (Cai et al. 2016):

$$q_{{{\text{O}}_{2} }} (C_{{{\text{O}}_{2} }}^{{{\text{in}}}} - C_{{{\text{O}}_{2} }}^{{{\text{ex}}}} ) = \frac{{J_{{{\text{O}}_{2} ,{\text{penetration}}}} A_{{{\text{ves}}}} }}{{V_{{{\text{el}}}} - V_{{{\text{ves}}}} }}$$

where V_el is the volume of the element, V_ves is the volume occupied by the blood vessel and A_ves is the area of the blood vessel wall. V_el, V_ves and A_ves are computed by V_el = l_el³, $V_{{{\text{ves}}}} = \frac{{\pi d_{{{\text{ves}}}}^{2} l_{{{\text{el}}}} }}{4}$ and A_ves = π d_ves l_el, respectively. l_el is the side length of the element, and d_ves is the diameter of the blood vessel. $J_{{{\text{O}}_{2} ,{\text{penetration}}}}$ is the free oxygen flux across the blood vessel wall, obtained by (Fang et al. 2008)

$$J_{{{\text{O}}_{2} ,{\text{penetration}}}} = - K_{{\text{w}}} \frac{{(C_{{{\text{O}}_{2} }}^{{{\text{ex}}}} - C_{{{\text{O}}_{2} }}^{{{\text{in}}}} )}}{\alpha w}$$

where α is the Bunsen solubility coefficient, w is the blood vessel wall thickness and K_w is the blood vessel wall permeability (Table

Table 5 The parameter values

5). Therefore, the coefficient $q_{{{\text{O}}_{2} }}$ is obtained by

$$q_{{{\text{O}}_{2} }} = \frac{{K_{{\text{w}}} A_{{{\text{ves}}}} }}{{\alpha w(V_{{{\text{el}}}} - V_{{{\text{ves}}}} )}} = \frac{{K_{{\text{w}}} \pi }}{{\alpha w\left( {\frac{{l_{{{\text{el}}}}^{2} }}{{d_{{{\text{ves}}}} }} - \frac{{\pi d_{{{\text{ves}}}} }}{4}} \right)}}.$$