Changes in scaffold porosity during bone tissue engineering in perfusion bioreactors considerably affect cellular mechanical stimulation for mineralization.

Bone tissue engineering (BTE) experiments in vitro have shown that fluid-induced wall shear stress (WSS) can stimulate cells to produce mineralized extracellular matrix (ECM). The application of WSS on seeded cells can be achieved through bioreactors that perfuse medium through porous scaffolds. In BTE experiments in vitro, commonly a constant flow rate is used. Previous studies have found that tissue growth within the scaffold will result in an increase of the WSS over time. To keep the WSS in a reported optimal range of 10-30 mPa, the applied external flow rate can be decreased over time. To investigate what reduction of the external flow rate during culturing is needed to keep the WSS in the optimal range, we here conducted a computational study, which simulated the formation of ECM, and in which we investigated the effect of constant fluid flow and different fluid flow reduction scenarios on the WSS. It was found that for both constant and reduced fluid flow scenarios, the WSS did not exceed a critical value, which was set to 60 mPa. However, the constant flow velocity resulted in a reduction of the cell/ECM surface being exposed to a WSS in the optimal range from 50% at the start of culture to 18.6% at day 21. Reducing the fluid flow over time could avoid much of this effect, leaving the WSS in the optimal range for 40.9% of the surface at 21 days. Therefore, for achieving more mineralized tissue, the conventional manner of loading the perfusion bioreactors (i.e. constant flow rate/velocity) should be changed to a decreasing flow over time in BTE experiments. This study provides an in silico tool for finding the best fluid flow reduction strategy.

Introduction

Mechanical stimulation in terms of fluid-induced wall shear stress (WSS) on bone cells can regulate extracellular matrix (ECM) mineralization in the presence of osteoinductive media (Giorgi et al., 2016; Wittkowske et al., 2016). Previous experimental studies sought to investigate this mechanobiological response of bone cells (osteoprogenitors/osteoblasts/osteocytes) by applying fluid-induced wall shear stress (WSS) on cells that were seeded on 2D substrates (Delaine-Smith et al., 2012; Mai et al., 2013; Michael Delaine-Smith et al., 2015). It was found that more mineral was deposited when applying a WSS in the range of 51–1200 mPa, compared to static culturing (Delaine-Smith et al., 2012; Mai et al., 2013; Michael Delaine-Smith et al., 2015). In 3D bone tissue engineering (BTE) in vitro, WSS is applied on cells seeded on scaffolds typically by perfusing a medium through the scaffold pores (Vetsch et al., 2015; Yeatts and Fisher, 2011). Previous BTE experiments have found that a WSS in a range of 10–30 mPa (Sikavitsas et al., 2003), or 0.55–24 mPa (Vetsch et al., 2017) can stimulate the cells to deposit mineralized ECM within 3D scaffolds. According to other studies, excessively high WSS in 3D scaffolds (i.e. >60 mPa) can cause cell death (Mccoy et al., 2012; Olivares et al., 2009). Therefore, controlling the WSS within the stimulating range is important for BTE in vitro. However, the WSS on cells within scaffolds depends not only on the loading conditions of the bioreactors (i.e. applied flow rate), but also on the scaffold geometry (i.e. pore size, pore shape and porosity) (Zhao et al., 2016, Zhao et al., 2018), and may be highly non-uniform. To avoid excessive trial-and-error experiments to obtain the best flow rate conditions for ECM mineralization, computational modelling approaches have been used for calculating the WSS on cells in perfusion bioreactors (Vetsch et al., 2017; Zhao et al., 2018). Some researchers idealized the scaffold geometry and used analytical estimations to calculate the local WSS (Blecha et al., 2010; Goldstein et al., 2001). Specifically, the scaffold pores were idealized as a material with cylindrical holes and constant porosity, and analytical equations then were used to calculate the WSS on the cylindrical hole surfaces. More detailed WSS quantifications were conducted by applying Computational Fluid Dynamics (CFD) approaches on scaffolds with more complex pore geometries (Marin and Lacroix, 2015; Olivares et al., 2009; Papantoniou et al., 2014; Zhao et al., 2016, Zhao et al., 2019). The geometries were obtained either from computer aided design (CAD) (Olivares et al., 2009; Zhao et al., 2016) or from micro-computed tomography (μCT) scanning (Marin and Lacroix, 2015; Santamaría et al., 2013; Zhao et al., 2019). The geometry of these computational models represented the empty scaffolds without considering cells/ECM growth in the pores. Recent in silico studies, which incorporated CFD approaches with ECM growth models have found that the WSS on mixed cell/ECM surface will change with the growth of ECM/cells (Guyot et al., 2015, Guyot et al., 2016a; Nava et al., 2013; Zhao et al., 2015). The maximum WSS on cells increased 2-fold with the pore volume filling from 10% to 60% by cell/ECM (Guyot et al., 2015). However, it is still unclear whether this WSS variation will exceed the optimal range for stimulating mineralization or not.

In previous BTE experiments, which used perfusion bioreactors to apply WSS on cells, typically flow rates were applied to the bioreactors by a peristaltic pump (Li et al., 2009; Sikavitsas et al., 2003; Vetsch et al., 2017). A recent computational study has optimized the external flow for maximizing the mineralization in BTE in vitro, and it was predicted that the optimal fluid peak velocity was in the range of 0.166–1.66 mm/s to get the WSS in the range known to stimulate ECM mineralization (Zhao et al., 2018). This flow velocity was computed by the CFD model based on the empty scaffolds, which represents the initial phase (e.g. the first week in cell culturing experiments) in BTE in vitro. However, the WSS on cells that are directly exposed to the medium flow increases with the growth of tissue (Guyot et al., 2015, Guyot et al., 2016a). Consequently, for maximizing the mineralization in BTE experiments in vitro, the external medium flow velocity (or flow rate) would have to decrease over the time.

The goal of the present study was to investigate what reduction of the flow applied to the bioreactor is needed to keep the WSS in a range that promotes ECM mineralization. A combined ECM growth model and CFD model was applied to calculate the change in WSS over time. Different scenarios for reducing the fluid flow were applied in the in silico model to investigate their effect of the ECM mineralization. The output from this study may help the BTE field to design in vitro experiments that use a decreasing flow rather than the conventional constant flow rate/velocity over time. Furthermore, this study also provides an in silico tool for real-time loading optimization for maximizing the amount of mineralized ECM.

Methods

In silico model development

In in vitro BTE, appositional tissue growth occurs when the cells are flatly attached to the scaffold surfaces and deposit ECM at this surface, so that ECM production results in a thickening of the scaffold struts. With the production of ECM, some cells will be encapsulated within the matrix, while others will attach on the surface of the matrix as illustrated in Fig. 1. Under perfusion flow, these cells on the surface of the matrix will be exposed to the WSS generated by the medium flow. In this study, we homogenized the mixture of ECM and cells (Fig. 1). It was assumed that the WSS at the surface of the homogenized cell/ECM is representative for the WSS sensed by the cell at the surface (Fig. 1).

Fig. 1

Schematic illustration of homogenization of cell and ECM mixture.

In the in silico model, the growth of cell/ECM followed a diffusion equation (Chaplain and Stuart, 1991; Nava et al., 2013). The ECM formation within the pore space of the scaffold was represented by a chance (or agent) parameter C that could range from 0 at the center of the void to 1 at the scaffold surface. A diffusion equation was applied for the pore space to generate a gradient of the formation chance C at time step 21 (day 21) (Fig. 2a). If the WSS τ ≥ 60 mPa no more ECM formation was allowed as cell death is expected for such high WSS values (McCoy and O'Brien, 2010; Olivares et al., 2009):with, D a diffusivity constant determining the rate. Bone ECM then was modelled in regions, for which the formation chance exceeded a threshold C (Fig. 2b).

Fig. 2

Schematic illustration of parameter exchange between (a) diffusion process and (b) rheological property (i.e. dynamic viscosity) adaptation (ECM/cells growth) in an idealized porous domain.

The local fluid velocity at position was obtained by solving the Navier-Stokes equation in ANSYS CFX solver (ANSYS Inc., PA, USA) using a finite volume method (FVM). In this study, the medium was modelled as incompressible fluid, thus the Navier-Stokes equation was expressed as:where, t the time; η is the fluid dynamic viscosity and p the fluid pressure.

The shear rate was then calculated from the fluid velocity:and the shear stress at the medium – homogenized cell/ECM interface г was calculated using:where, the medium – cell/ECM mixture interface г was tracked by defining the element nodes with the value of C = C0 as shown in Fig. 2(b).

The ECM formation was modelled by changing the rheological property (i.e. viscosity) of the element. For regions representing fluid, the element viscosity was set to a value η = 1 mPa·s representing Dulbecco's Modified Eagle medium (Maisonneuve et al., 2013). For regions representing ECM, the element viscosity was set to η = 2 Pa·s (Köster et al., 2008). In case the WSS exceeded 60 mPa, it was assumed that ECM would be flushed out and no new formation would take place. This resulted in the following equation for the element viscosity:

Eqs. (1), (4), (5) were formulated using CFX Expression Language (CEL), and put into the CFX solver.

In silico simulation and parameter tuning

To apply the algorithm (Section 2.1) in simulating the matrix formation and calculating the WSS on the homogenized ECM/cell surface, we developed a CFD model based on a perfusion bioreactor system (Fig. 3). The bioreactor had a cylindrical shape with a radius of 3 mm, while the scaffold with a lattice microstructure was located at the center and had a length of 6 mm. Because of symmetry, only one quarter section of the bioreactor and scaffold was modelled. Corresponding to experimental conditions (Papantoniou et al., 2014), a constant fluid velocity of 1.2 mm/s was applied to the inlet of the bioreactor chamber as shown in Fig. 3(a), and free outflow was prescribed at the outlet. The outer side surface and scaffold surfaces were defined as non-slip walls. On the two symmetry plane surfaces (Fig. 3a), the fluid velocity followed Eq. (6):with j being the direction perpendicular to the cutting surface.

Fig. 3

(a) One quarter of the fluid domain in the bioreactor for CFD analysis: the up surface and bottom surface are defined as inlet and outlet, respectively; two side cutting surfaces are symmetric boundaries; the outer side surface and scaffold surfaces are defined as non-slip walls, (b) scaffold geometry with pore size and porosity equal in the three orthogonal directions.

The whole fluid domain was meshed by 1.104 × 106 tetrahedral elements with a patch conforming method. The CFD model was solved under transient state. A FVM was employed by ANSYS CFX solver to resolve the model under the convergence criteria of root-mean-square residual of the mass and momentum <1 × 10−4.

To simulate the cell/ECM growth, two parameters (D0 and C0) needed to be determined. Firstly, an arbitrary value of 1.0 × 10−13 m2/s (~1.0 × 10−8 m2/day), which was within the reported range of 3.3 × 10−14 m2/s–8.3 × 10−13 m2/s (Croll et al., 2005), was assigned to D0. This D0 was used to generate the gradient of parameter C in the CFD domain (C = 0–1.00 at day 21). Afterwards, parameter C0 was varied in the range between 1.00 and 0.01 to find the value that resulted in the best fit to the measured ECM volume in experimental results described in literature (Papantoniou et al., 2014).

Applied fluid velocity variation

To investigate what reduction of the flow velocity can generate a WSS in the range for stimulating ECM formation and mineralization, the prescribed fluid velocity was made dependent on time according to:where, V is the inlet fluid velocity used in the experiment (1.2 mm/s) (Papantoniou et al., 2014) and k a coefficient [mm/s day] that determines the decrease in speed over time. Seven different values were chosen for k, ranging from 0 to 0.059 mm/s/day. A value of k = 0 implies a constant fluid velocity, while increasing values represent a linear decrease over time (Fig. 4).

Fig. 4

Seven types of applied inlet fluid velocities (V), where type 1 (k1 = 0 mm/s day) is constant fluid velocity, types 2–7 have a decreasing velocity profile from day 1 to day 21.

Results

According to the results reported in a previous experimental study (Papantoniou et al., 2014), 45% and 65% of the porous volume were filled with ECM/cells at week 2 (day 14) and week 3 (day 21), respectively. The results from parameter (C0) tuning showed that at day 21, the tissue volume fraction increased by approximately 3% with 0.01 increment of parameter C0. The best agreement with the results at day 21 was found when setting parameter C0 = 0.12, which led to a volume filling fraction of 65% in the in silico model. Using this value, a volume filling fraction of 47% was predicted for day 14, which was in good agreement with the experimental result (Papantoniou et al., 2014). A visual comparison of tissue growth throughout the scaffold, which was from our in silico prediction and the previously reported in vitro experiment (Papantoniou et al., 2014) was presented in Supplementary Fig. 1.

Under constant fluid velocity (V1 = 1.2 mm/s) as applied in the experiment (Papantoniou et al., 2014), a 3D visualisation revealed that high WSS were mainly found along the surfaces of the vertical struts while lower values were found for the surfaces of the transversal struts (Fig. 5). The results also showed that the WSS on the cell/ECM surface increased over time (Fig. 5b). At 21 days, most of the ECM/cell surface was subjected to a WSS higher than 30 mPa, which was considered the higher end of the optimal range. The simulations predicted that the homogenized ECM/cell volume kept increasing from day 2 to day 21 (from 7.5% to 65.0%) under constant velocity (V1 = 1.2 mm/s). However, the surface area fraction of ECM/cell that experienced a WSS of 10–30 mPa (for mineralization) had a decreasing trend (Fig. 6) due to the increasing local WSS at the cell/ECM surface.

Fig. 5

(a) Simulated pore volume filling by homogenized ECM/cell at day 7, 14 and 21, the ECM/cell volume increases within 21 days, A-A is a zoomed-in top view that shows the cell/ECM growth within one unit scaffold, (b) increasing wall shear stress distribution at the interface between medium and ECM/cell at days 7, 14 and 21, the geometry is the fluid volume (counterpart of the cell/ECM volume within the pores).

Fig. 6

Evolution of ECM/cell surface area fraction that undergoes the WSS of 10–30 mPa (for stimulating mineralization) under different loading conditions of the perfusion bioreactor.

The cell/ECM surface area fraction subjected to a WSS between 10 and 30 mPa was increased when reducing the fluid flow over time (i.e. with increasing k values) from day 10 onwards according to Fig. 6. The highest cell/ECM surface area fraction within the WSS range of 10–30 mPa was 40.9% (day 21) for fluid velocity profile V3 (decreasing from 1.2–0.4 mm/s from day 1 to day 21). Under constant fluid velocity (V1), the cell/ECM surface area fraction was only 19.3% at day 21. This demonstrated that using a decreasing flow velocity may result in more ECM mineralization than constant velocity. However, the results were highly dependent on the rate of the fluid velocity decrease. For example, for fluid velocity profile V6, only 18.6% of the cell/ECM surface had a WSS in the proper range at day 21, which was even less than the constant loading profile (V1). This was due to the fact that the WSS became too low too fast (0.02 mm/s at day 21). Therefore, to maintain larger portion of cells undergoing the optimal WSS range of 10–30 mPa within 3 weeks, the decreasing speeds of 0.04 mm/s per day and 0.048 mm/s per day are the preferable ones from the initial flow velocity of 1.2 mm/s.

Discussion and conclusion

The goal of the present study was to investigate whether a reduction of the flow applied to a perfusion bioreactor over time can keep a larger amount of WSS in a range to promote ECM mineralization. The first result of this study was that keeping the fluid flow rate constant would result in a considerably increased WSS over time due to filling of pores. This result was in agreement with results of earlier studies that showed that the WSS on the neo-tissue surface increased with tissue growth within the scaffold (Guyot et al., 2015, Guyot et al., 2016a). It was implied that calculations of the WSS using analytical models or CFD analyses based on the scaffold geometry (without tissue in it) could only predict the initial WSS. In our study, when using a constant fluid velocity of 1.2 mm/s, the cell/ECM surface area fraction being stimulated for mineralization (WSS = 10–30 mPa) decreased from 50.11% to 19.26%. It was indicated that with constant fluid flow 30.8% of the cell/ECM surface that initially was exposed to the proper range now has exceeded a WSS of 30 mPa at day 21.

The second result of this study was that a linear reduction of the applied fluid velocity over time resulted in a considerably larger surface fraction being subjected to a WSS in the range known to stimulate ECM/mineralization formation. For the best case investigated here, at day 21 a surface fraction of 40.9% was subjected to a WSS in the range of 10–30 mPa which was considerably higher than the surface fraction of 19.3% found under the constant fluid flow. The rate by which the fluid flow should be reduced, however, needs to be optimized based on the estimated or measured ECM formation. If the fluid flow was reduced too much, the WSS would become too low. For studies that use 3D micro-CT imaging of the scaffold over time (Vetsch et al., 2016, Vetsch et al., 2017), it would be possible to measure at least the mineral formation and use this information to calculate an updated fluid flow rate after scanning. However, non-mineralized ECM formation would not be visible using micro-CT. Hence, estimates of the expected ECM formation rates will likely be needed to get accurate predictions.

There were several limitations to our study. First, the cells were assumed to be flatly attached on the scaffold surface or the earlier formed ECM all the time. A previous in vitro experimental study, in which cells were cultured in pores of different sizes found that cells would detach from the pore surface when the pore size was smaller than 500 μm (Herklotz et al., 2015). In such small pores the cells bridged across the pores and produced ECM that infiltrated within the pores (i.e. interstitial growth), as was observed in other BTE experiments (Li et al., 2009; Melke et al., 2018; Woloszyk et al., 2014). Our in silico model was based on the experimental setup in (Papantoniou et al., 2014), in which the scaffold had a uniform pore size of 755 μm (>500 μm). Thus, appositional growth of cell/ECM was applied in our in silico model, assuming that interstitial growth was absent. Nevertheless, during culturing, the pore size would be reduced by the newly formed tissue and thus bridging may become possible. This will further reduce the porosity, hence the permeability, of the scaffold, and further increase the fluid velocity. It thus was possible, that a further reduction of fluid flow was needed once the pore size gets below 500 μm. In addition, in our model, the cells and ECM were homogenized, and ECM formation was modelled only by changing the element viscosity (Eqs. (6), (7)). Previously, Guyot et al., 2015, Guyot et al., 2016b used another technique, in which the neo-tissue was modelled as a porous media. A limitation of using this technique, however, was the unknown permeability of the neo-tissue. Third, in our study, the prediction of ECM mineralization was based on a previously reported mechano-regulation theory, wherein a WSS of 10–30 mPa on the stem cells surface was optimal for stimulating mineralization within a titanium scaffold (Sikavitsas et al., 2003). It has been proposed that not only the WSS but also the mechanical strain could stimulate stem cells through enhancing osteogenic differentiation and mineral deposition on 2D substrates (Virjula et al., 2017; Yu et al., 2010). However, the findings of 2D mechanobiological experiments were not likely to be directly translated to 3D tissue engineering experiments (McCoy and O'Brien, 2010). Therefore, the influence of such internal stress/strain within cell-embedded matrix on stimulating mineralization of ECM in tissue engineering (e.g. 3D cell culturing environment) is still unclear. Fourth, considering the available experimental data for validation, the WSS-enhanced cell/ECM growth rate was not included in our in silico model. Previously, Chapman et al. (2014) proposed a mathematical expression for WSS-enhanced tissue growth rate, which happened in the WSS range of 0–10 mPa. The aim of this study was to find an optimal decreasing flow rate that could generate a WSS of 10–30 mPa on cells/ECM. However, under this WSS range, the cell/ECM growth rate was constant, according to (Chapman et al., 2014). Nevertheless, the accuracy of simulated tissue growth volume still would be affect. The possible reason could be due to the specific scaffold pore shape, the local WSS at some locations in scaffold still had the possibility to fall within 0–10 mPa under optimal flow rate, even though WSS in majority of regions was in the range of 10–30 mPa. Therefore, to build a more precise in silico model for predicting tissue growth, more factors (e.g. WSS-enhanced tissue growth) need to be considered in the in silico model by comparing to the respective experimental data. Therefore, to build a more precise in silico model for predicting tissue growth, more factors (e.g. WSS-enhanced tissue growth) need to be considered in the in silico model by comparing to the respective experimental data. Finally, as the theory for mineralization used here was based on the WSS, we focused on simulating dynamic culturing condition (e.g. applying WSS on cells) in this study. Therefore, this in silico model is not valid for predicting the ECM mineralization under static condition.

As a conclusion, a constant flow rate, which was commonly applied in BTE experiments could result in a too high WSS on surface cells once a considerable amount of ECM has been formed. To keep the WSS in an optimal stimulation range for the cells located at the tissue surface, the fluid velocity should be reduced over time.

The following is the supplementary data related to this article.

Supplementary Fig. 1

Tissue growth at day 1, 14 and 21, which are from in silico model and the previously reported in vitro experiment (Papantoniou et al., 2014).

CRediT authorship contribution statement

Feihu Zhao:Methodology, Resources, Formal analysis, Writing - original draft, Writing - review & editing.Damien Lacroix:Methodology.Keita Ito:Methodology, Formal analysis.Bert van Rietbergen:Methodology, Formal analysis.Sandra Hofmann:Methodology, Formal analysis.

Declaration of competing interest

The authors declare that there is no conflict of interest in this study.