Numerical Study of Granular Scaffold Efficiency to Convert Fluid Flow into Mechanical Stimulation in Bone Tissue Engineering

Abstract

Controlling the mechanical environment in bioreactors represents a key element in the reactors' optimization. Positive effects of fluid flow in three-dimensional bioreactors have been observed, but local stresses at cell scale remain unknown. These effects led to the development of numerical tools to assess the micromechanical environment of cells in bioreactors. Recently, new possible scaffold geometry has emerged: granular packings. In the present study, the primary goal was to compare the efficiency of such a scaffold to the other ones from literature in terms of wall shear stress levels and distributions. To that aim, three different types of granular packings were generated through discrete element method, and computational fluid dynamics was used to simulate the flow within these packings. Shear stress levels and distributions were determined. A linear relationship between shear stress and inlet velocity was observed, and its slope was similar to published data. The distributions of normalized stress were independent of the inlet velocity and were highly comparable to those of widely used porous scaffolds. Granular packings present similar features to more classical porous scaffolds and have the advantage of being easy to manipulate and seed. The methods of this work are generalizable to the study of other granular packing configurations.

Introduction

Bone tissue engineering strategies for constructs of clinically relevant sizes require dynamic cultures in bioreactors.^1–4 These devices provide nutrients, oxygen, and mechanical stimuli to mesenchymal stem cells (MSCs) seeded on three-dimensional (3D) biomaterial scaffolds. MSCs are known to be shear stress sensitive.^5–14 Bioreactors have been shown to increase MSC proliferation^4,15 and osteogenesis.^2,15–18 In these systems, cell feeding and mechanostimulation are ensured by culture media directly perfused through pores in the 3D scaffold. Therefore, the design of scaffold geometry is a key challenge to regulate shear stress in 3D constructs.

A large diversity of scaffolds have been used and investigated. Highly porous foam scaffolds were extensively used in perfusion bioreactors,^17–22 but fiber mesh scaffolds^21,23–26 and scaffolds with regular architecture^27–29 were also used by several authors. All these scaffolds were massive and porous and cells were seeded by perfusing a cell suspension. In this large variety of 3D porous scaffolds, low inlet flows lead to cell proliferation,^17,21 upregulation of genes involved in the mechanotransduction process,^19,22 and increase of alkaline phosphatase activity.^16,26

However, due to complex internal geometries and flow conditions, mechanical stresses in scaffolds are not homogeneously distributed and cannot be analytically predicted. In this context, numerical tools have been widely used in an attempt to evaluate shear stress levels and distributions within actual scaffolds. Geometries were either determined through microcomputed tomography or created with geometrical objects. Simulations were performed using computational fluid dynamics^{20,23,30–33} or Lattice-Boltzmann methods.^21,34 At efficient inlet flow, the levels of shear stress evaluated with simulation techniques ranged from 10⁻⁵ to 10⁻² Pa. In 2010, Voronov et al.³⁵ studied an extensive set of porous scaffolds and determined a specific three-parameter gamma function to predict in a simple way probability density functions (PDFs) of all scaffolds. This function was also validated for other porous scaffolds of the literature.

Recently, packings of particles were also considered as new possible scaffold configurations.^36–40 Instead of being massive, they are composed of an assembly of several biomaterial particles on which cells are seeded. In a clinical context, such scaffolds would be more adaptable to the surgical vagaries, in terms of size and shape. They would also be easier to seed and manipulate. Perfusion of such granular scaffolds with a culture medium had positive effects on cell proliferation,^36–38,40 bone formation markers,^36,37 and mineralization.^36–38 However, to date and to the best of our knowledge, wall shear stress has not been investigated through 3D simulations.

In the present work, a packed-column bioreactor composed of a vertical cylinder filled with granular particles was studied. Three types of granular packings were investigated: a packing of cubes (previously experimentally studied in the work of David et al.³⁶), a packing of beads to investigate the effect of geometry, and a bidisperse packing of beads to investigate the effect of porosity. To this aim, 3D computational fluid dynamics was used to simulate the flow within the packed-column bioreactor. The objective was to determine numerically wall shear stress levels and distributions inside these scaffolds and to compare these results with porous scaffolds, as depicted in literature. On a biological point of view, the results of this numerical study will represent a key step in the optimization of the mechanical microenvironment of cells in 3D dynamic culture. Distributions of fluid velocity and wall shear stress on the particle surface were analyzed as a function of inlet flow and particle geometry. The mechanical efficiency of fluid flow was evaluated, and heterogeneity was compared with the three-parameter gamma function proposed by Voronov et al.

Materials and Methods

The aim of this study was to assess mechanical stimulation perceived by cells in a bone tissue engineering packed-column bioreactor³⁶ designed as a 13.1-mm-diameter vertical Plexiglas^® tube filled with coral particles. In this bioreactor, cells were seeded on the surface of the particles and a continuous unidirectional flow of alpha-MEM culture medium was imposed.

Generation of granular packings

Simulation of the flow in a bone tissue engineering packed-column bioreactor requires full knowledge of the scaffold geometry composed by different random packings of particles. In the present study, three types of particles were used: (1) 3-mm-side cubes, (2) 2-mm-diameter beads, and (3) 3.5-mm-diameter beads. In each case, the particles were assumed to be perfect on a geometrical point of view and sufficiently large in regard of typical cell size to assume a flat cell configuration. Three granular packings of cubic and spherical particles were generated with these three types of particles. The three configurations were composed of (1) 51 cubes with a side of 3 mm denoted Ccube, (2) 71 beads with a diameter of 3.5 mm denoted C1bead, (3) 36 beads with a diameter of 3.5 mm mixed with 108 beads with a diameter of 2 mm denoted C2bead. In each configuration, the number of particles was determined with respect to (1) a representative volume of cubic particle packing corresponding to a minimal size defect of 15 mm length and (2) an equal overall surface of the particles close to 2740 mm² for every configuration, thus providing the same surface for cell attachment. To create all granular packings, the particles were randomly picked one by one and left to fall under gravity in a tube representing the cylindrical bioreactor. This numerical procedure reproduced in a realistic manner the experimental conditions of the filling of the bioreactor. For that purpose, the LMGC90 software based on the discrete element method was used. This software was mainly developed for granular material simulation and all details of this simulation method can be found in Jean⁴¹ and Radjai and Richefeu⁴² For each configuration, four random assemblies were generated leading to 12 granular packings studied in the present work. An example of each configuration can be seen in Figure 1.

FIG. 1.

Typical velocity streamlines obtained for one assembly of each configuration Ccube (a), C1bead (b), and C2bead (c) at 200 mL/min. The line's color gives the norm of the local velocity vector. The fluid circulated from bottom to top of packing. Color images available online at www.liebertpub.com/tec

As the overall surface of the packing was fixed to be equal for all configurations, the average height of C2bead (h=17.6 mm) was significantly smaller than the one of Ccube (h=22.4 mm) or C1bead (h=23.6 mm). The porosity corresponding to the ratio of fluid volume to total packing volume was 54.2%, 49.8%, and 46.9% for Ccube, C1bead, and C2bead configurations, respectively.

Thanks to the data of positions and orientations obtained with LMGC90, the 12 geometries were created in the numerical model and then meshed with the software Comsol Multiphysics 3.4. Particles were considered to be fixed and nondeformable. In all granular packings, the same density of elements on the surface of the particles was imposed to a value close to 40 triangular elements/mm². This density of surface elements was identified after a preliminary convergence study allowing both accurate simulations and shortened calculation times. Then, the fluid domain composed by the whole cylinder minus the particles was meshed by 3D tetrahedral elements (10-node tetrahedrons, Lagrange P2-P1). In all configurations, the mean tetrahedron density was set to ∼200 tetrahedrons/mm³.

Computational fluid dynamic simulations

In all simulations, the fluid was considered to be homogeneous, Newtonian, and incompressible with a density ρ=1000 kg/m³ and a dynamic viscosity η=1.45×10⁻³ Pa·s.⁴³ The fluid followed a classical Poiseuille flow at the inlet bottom boundary and the inlet velocity was given by the parabolic equation (1): \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}v ( r ) = \frac { 2 \varphi } { \pi R^2 } \ast \left( 1 - \frac { r^2 } { R^2 } \right) \tag { 1 } \end{align*} \end{document}

where r is the radial position of the fluid particle, ϕ is the inlet flow, and R is the radius of the cylinder. In the present study, five inlet flows were tested: 10, 50, 100, 150, and 200 mL/min corresponding to a mean inlet velocity of 1.2, 6.2, 12.4, 18.5, and 24.7 mm/s, respectively. The top boundary was defined as an open boundary where the atmospheric pressure was imposed. A classical no slip condition was used on all other boundaries corresponding to walls. The flow was simulated in three dimensions inside the bioreactor by solving the Navier–Stokes equations in a steady-state analysis.

Analyses

Analyses were focused on velocity and shear stress values and distributions on the surface of particles. Velocity streamlines were displayed to observe the fluid behavior within the packings of particles, while the wall shear stress norms were computed on every node of the surface of the particles with the following equation (2): \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} \tau = \mid \overrightarrow{T} - ( \overrightarrow{T} . \overrightarrow{n} ) \overrightarrow{n} \mid \tag{2}\end{align*} \end{document}

where \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\overrightarrow{T} = \sigma. \overrightarrow{n}$$ \end{document} is the stress vector, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\overrightarrow{n}$$ \end{document} is the normal vector of the surface, and σ is the Cauchy stress tensor given by the classical behavior law for a Newtonian incompressible viscous fluid. In all distributions, every wall shear stress value was balanced with the local size of the mesh. Due to the Reynolds number values (<230) calculated in every configuration, a laminar flow was assumed.

The wall stress PDFs for all computational fluid dynamics (CFD) simulations were also generated. Wall shear stresses were normalized using equation (3): \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}\tau^ { \ast } = \frac { \tau - \bar \tau } { \sigma_ { \tau } } \tag { 3 } \end{align*} \end{document}

where \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\bar \tau$$ \end{document} was the mean wall shear stress and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\sigma_{ \tau}$$ \end{document} was the standard deviation of the stress distribution. The obtained dimensionless wall shear stresses were then distributed into 50 equally sized bins, ranging from the smallest to the largest value using MATLAB software. The three-parameter gamma PDF Γ(2.91,−1.43,0.45) defined by Voronov et al.³⁵ was compared with experimental PDFs.

Heterogeneity

Any reliable tissue engineering process requires a homogeneous local stress field in the whole of the 3D packing. A radial heterogeneity variable was defined to evaluate the variability of stresses on the particle walls. For each calculation, the model was discretized and mean wall shear stresses were calculated for 20 ranges of radial position. The radial heterogeneity H_r was defined as the standard deviation of these 20 values reported to the overall mean value, expressed in percentage.

Statistical analysis

For each configuration and each flow, four initial random packings were generated and analyzed. In the present study, the mean value corresponded to the average value from the four simulations and the error bars plotted were plus and minus one standard deviation. Statistical analysis was performed using R (the R foundation) with a threshold of 5% corresponding to the alpha risk. A nonparametric Kruskal–Wallis test was used to detect any significant differences between the configurations. If at least one of the configurations significantly differed from the others, a post hoc Conover–Iman test was used to compare the configurations two by two. The influence of the inlet flow on the mean wall shear stress was studied using a Spearman correlation test.

Results

Velocity results

A typical result of 3D velocity field for one assembly of each configuration at 200 mL/min is given on Figure 1. The mean velocity values were 0.045, 0.053, and 0.054 m/s for Ccube, C1bead, and C2bead, respectively. These values should be compared to the mean velocity imposed at the bottom boundary of 0.025 m/s. The highest levels of velocity were recorded in the C2bead configuration. Fluid velocity levels within the construct are dependent on the particles geometry and arrangement.

Shear stress levels

The wall shear stresses on the surfaces of particles are shown in Figure 2 for one assembly of each configuration for every inlet flow. Wall shear stress increased when passing from Ccube to C1bead and to C2bead with a same inlet flow. Obviously, wall shear stress also increased when the inlet flow increased. Localizations of maximum wall shear stresses were similar whatever the inlet flow.

FIG. 2.

Wall shear stress obtained for one assembly of each configuration Ccube, C1bead, and C2bead at every inlet flow. The fluid circulated from bottom to top of packing. Color images available online at www.liebertpub.com/tec

Figure 3 shows the mean wall shear stress versus the inlet flow for all configurations. Each bar corresponds to the average value of the four granular packing simulations with the associated standard deviation. As expected, the mean wall shear stress increased with the inlet flow. In the considered range of flows, the relationship between the mean wall shear stress and the inlet flow was linear. The slopes were 16.1 (p<0.0001, R²=0.995), 27.5 (p<0.0001, R²=0.995), and 38.6 (p<0.0001, R²=0.998) Pa·s·m⁻¹ for Ccube, C1bead, and C2bead, respectively. Mechanical efficiency of the inlet flow obtained through this slope depended strongly on the configuration. At the same inlet flow, the mean wall shear stress increased by 75% and 40% switching from Ccube to C1bead and C1bead to C2bead, respectively.

FIG. 3.

Wall shear stress for Ccube, C1bead, and C2bead configurations at 10, 50, 100, 150, and 200 mL/min. Each bar represents the average value calculated for the four assemblies of each configuration.

Shear stress distributions

On a local scale, the inlet flow played a significant role on the distributions of the wall shear stress applied at the surface of the particles. The distributions are shown in Figure 4 where each bar corresponds to the average of the four granular packings with the associated standard deviation. The distribution for an inlet flow of 10 mL/min was not presented because 100% of the shear stress values were included in the first range of wall shear stress (0–0.2 Pa) for all configurations. As expected, the local wall shear stress increased with the inlet flow. At each flow, the distribution for C2bead was more skewed toward higher shear stresses than the distributions for Ccube and C1bead.

FIG. 4.

Wall shear stress distribution for Ccube, C1bead, and C2bead at 50 mL/min (a), 100 mL/min (b), 150 mL/min (c), and 200 mL/min (d). Each bar represents the average value of the four assemblies of each configuration.

Local dispersion analysis

The obtained PDFs for all three geometries of this study are presented in Figure 5. The 1000 dots corresponding to the 50 bins of all four assemblies and five inlet flows are presented in the same figure. All three distributions have a gamma-3p shape. Specifically, Ccube and C1bead both present high frequencies for low shear stresses resulting in a high and narrow peak on the left. This peak is higher in Ccube (0.83) than in C1bead (0.67). The C2bead PDF is wider with a lower maximum value (0.53) and an increased skewness toward higher shear stresses. The dashed line on every graph represents the PDF of stresses Γ(2.91,−1.43,0.45) obtained by Voronov et al. for 36 home-manufactured highly porous scaffolds studied. In the present study, a Kolmogorov–Smirnov goodness-of-fit test showed that the PDF from all configurations followed Voronov's gamma-3P distribution Γ(2.91,−1.43,0.45) with a significance level of 25–30%.

FIG. 5.

Probability density functions followed by the nondimensional wall shear stress obtained for Ccube (a), C1bead (b), and C2bead (c). Fifty bins were created for each calculation resulting in a total of 1000 dots per configuration. The dashed line represents the gamma distribution obtained by Voronov et al.: Γ(2.91,−1.43,0.45).

Shear stress heterogeneity

Figure 6 presents the radial heterogeneities obtained for all configurations and all inlet flows. For every inlet flow, H_r was of the order of 15% in C2bead, lower than in the two other configurations (20–30%). For all geometries, H_r slightly decreased when the inlet flow increased.

FIG. 6.

Radial heterogeneities calculated for Ccube, C1bead, and C2bead configurations at 10, 50, 100, 150 mL/min, and 200 mL/min. Each bar represents the average value calculated for the four assemblies of each configuration.

Discussion

In the present study, three porous scaffolds made of cubic or spherical particles of millimeter size were studied. Local shear stresses were calculated using the CFD method for the three granular packings. This present 3D numerical study on shear stress in granular scaffolds highlighted some interesting trails for the future development of such scaffolds. The particle geometry strongly influences fluid velocity, wall shear stress mean values and distributions, mechanical efficiency of the flow, and plays a role on homogeneity of shear stresses within the constructs.

As expected for each packing, a linear relationship between the mean shear stress at the wall and the mean inlet velocity was highlighted. The slope of this linear relationship, denoted S_w thereafter, provides the capability of the scaffold to convert the effect of the fluid flow into mechanical stimulation. Table 1 gives the slopes S_w calculated for the scaffolds of the present study and different scaffolds analyzed in the literature. In this study, S_w appears to be strongly dependent on the scaffold geometry and on the distribution of particle sizes. In particular, the slope of C2bead configuration is twice greater than the value obtained for Ccube configuration.

Table 1.

Slopes S_w Obtained for Different Scaffolds from Previous Studies and Their Porosity and Pore Size Properties

Reference	Scaffold type	Porosity Φ (%)	Mean diameter pore size d (μm)	Slope S_w (Pa·s/m)
Zermatten et al. (2014)³³	Foam–SF	55.5	160	56
Zermatten et al. (2014)³³	Foam–PCL	37.9	220	46
Maes et al. (2012)³²	Foam–Ha	73	270	41
Maes et al. (2012)³²	Foam–Ti	77	280	32
McCoy et al. (2012)⁴⁴	Foam–CG	90	85–325	75
VanGordon et al. (2011)²¹	Foam–PLA	85	271	79
Lesman et al. (2010)⁴⁵	Foam–PLA	90	500	30
Lesman et al. (2010)⁴⁵	Foam–PLA	90	150	105
Milan et al. (2009)⁴³	Foam–PLA	90	—	30
Jungreuthmayer et al. (2009)³¹	Foam–CG	90.5	96	80
Jungreuthmayer et al. (2009)³¹	Foam–CaPo	58.5	350	29
Cioffi et al. (2008)²⁰	Foam–PEGT	80	180	25
Cioffi et al. (2006)⁴⁶	Foam–Degrapo	77	100	68
Boschetti et al. (2006)²⁷	Foam–PLA	77	100	78
Present study Ccube	Granular cubes	54.2	—	16
Present Study C1bead	Granular beads	49.8	—	28
Present Study C2bead	Granular beads	46.4	—	39
VanGordon et al. (2011)²¹	Nonwoven fiber	85	—	79
Raimondi et al. (2002)²³	Nonwoven fiber	—	—	67

In the literature, a large number of studies focused on porous foam used as scaffold. These foams are composed with different materials, with a porosity Φ varying from 38% to 90% and a pore size diameter d from 100 to 350 μm. Although highly dispersed, all these foams gave the same order of magnitude of the slope with an average value close to 50 Pa·s·m⁻¹. Goldstein et al.⁴⁷ proposed a simplified cylindrical pore model leading to an analytical evaluation of the slope inversely proportional to the pore diameter size and the porosity. This formula roughly predicted the variation observed in the different configurations. Simulations of scaffolds with fiber architecture were scarcer, but the two studies referenced in the table gave an average slope slightly higher than the slope value measured in foam.

The slopes of the granular packings of the present study were smaller compared with literature data, but on the same order of magnitude. These variations cannot be explained by the small porosity decrease observed in granular configurations, suggesting a strong size effect of void space between particles. The typical morphology of granular packing led to a poor definition of pore diameter size. However, the void space was reduced when the number of particles and their contact points increased in a same volume fraction. In this context, the slope of C2bead configuration was naturally higher with more numerous particles. This result suggests that the slope of granular packing could be improved either by using a larger distribution of particle diameters or by using smaller particles. However, granular size diameter governed the rate of dissolution of the construct and needed a minimum size threshold.

Recently, Voronov et al. developed a new method based on PDFs to evaluate a priori the stress distribution of highly porous scaffolds. In particular, these authors showed that three-parameter gamma PDF Γ(2.91,−1.43,0.45) fits PDF of stresses obtained in 36 home-manufactured highly porous scaffolds studied and validated other porous foam scaffolds from different laboratories.^20,48,49 This new method was validated on granular packings used in the present study showing that PDF from all configurations followed the same gamma-3P distribution. This study proved that granular packing scaffolds followed the same PDF as massive porous scaffolds. In Ccube and C1bead, a high frequency peak was observed close to −1 normalized stress suggesting a more ordered microstructure.⁵⁰ To date and in our knowledge, PDF of wall shear stresses in scaffolds with porosity of the order of 45–55% was only studied by Pham et al.⁵⁰ In this work, authors found that for the most random configurations and with a porosity lower than 60%, the PDF still followed the gamma-3P distribution but with a shape modification: a high frequency peak close to −1 normalized stress occurred, the decrease classically observed at the left of this point disappeared, and this effect was more pronounced as porosity was decreased. Based on these results, the absence of any peak in the C2bead configuration revealed, as expected, a more random geometry and suggested that other bead size distributions leading to smaller porosity could be acceptable.

The heterogeneity analysis showed that the C2bead configuration provided more homogeneous stress within the constructs compared to the two other configurations. It is likely that the addition of smaller beads in C2bead compared to C1bead reduced and homogenized pore size and shape in the packing and, consequently, strongly reduced the heterogeneity. This is very important from the point of view of bioreactor design as the homogeneity of the stresses is directly related to the homogeneity of cell adhesion and production of extracellular matrix.

In the present study, the initial scaffold geometry was defined using a numerical Discrete Element Method consisting of randomly picked particles one by one into a set of defined particles. This original procedure was rapid and easy to implement compared with the microcomputed tomography approach. One limitation of the method used in the present study was the perfect characteristics of the particles and the lack of roughness at the surface of particles. However, to our knowledge, the submicrometer roughness effect on wall shear stress was not analyzed in the different porous scaffold simulations due to limits imposed by the microcomputed tomography resolution. Taking the roughness into account was not within the scope of this study, but creating a setup of rough particles in place of perfect particles is technically achievable.

In 3D scaffolds, beneficial fluid flow effects were observed and reported in literature for shear stress values ranging from 10⁻⁵ to 10⁻² Pa. Cell proliferation, production of extracellular matrix, activation of proteins, and expression of some genes related to mechanotransduction were frequently reported in literature with these stress levels.^{19,22,24–26} In the present study, such levels of stress were reached at 10 mL/min in all configurations. However, this level of stimulation is rather different from the values of shear stress in vivo 0.8–3 Pa⁵¹ or in two-dimensional experiments ranging from 0.3 to 2 Pa.^12,52,53 The device studied in the present work, a packed-column bioreactor filled with a granular packing scaffold, enables levels of stimulation of this order of magnitude to be reached. The high inlet flows required to reach these levels indicate that intermittent flows are to be planned. A discontinuous loading composed with alternating sequences of low and high flows would be comparable to in vivo loadings with alternative activity and rest phases.

In conclusion, scaffolds composed with granular packings were equivalent in terms of fluid flow mechanical efficiency and stress distribution, but would be advantageous in a clinical context. The C2bead configuration is the best among the three tested in the present work in terms of efficiency and homogeneity. The numerical method proposed in this work could be directly applied to the distribution of other particles to optimize scaffold regarding these criteria. In future experimental work, effects on cell proliferation, bone formation markers, and mineralization of the three granular packing will be analyzed.

Footnotes

Acknowledgments

This project was supported by the French National Agency. The authors acknowledge financial support from ANR OMBIOS (ANR-11-BS09-0036). They also wish to thank the IVTV platform (ANR-10-EQPX-06-01).

Disclosure Statement

No competing financial interests exist.

References

Bancroft

G.N.

, Sikavitsas

V.I.

, van den Dolder

, Sheffield

T.L.

, Ambrose

C.G.

, Jansen

J.A.

, and Mikos

A.G.

Fluid flow increases mineralized matrix deposition in 3D perfusion culture of marrow stromal osteoblasts in a dose-dependent manner. Proc Natl Acad Sci U S A, 99, 12600, 2002.

Bancroft

G.N.

, Sikavitsas

V.I.

, and Mikos

A.G.

Design of a flow perfusion bioreactor system for bone tissue-engineering applications. Tissue Eng, 9, 549, 2003.

Bilodeau

, and Mantovani

Bioreactors for tissue engineering: focus in mechanical constraints. A comparative review. Tissue Eng, 12, 2367, 2006.

Volkmer

, Drosse

, Otto

, Stangelmayer

, Stengele

, Kallukalam

B.C.

, Mutschler

, and Schieker

Hypoxia in static and dynamic 3D culture systems for tissue engineering of bone. Tissue Eng Part A, 14, 1331, 2008.

Donahue

T.L.H.

, Haut

T.R.

, Yellowley

C.E.

, Donahue

H.J.

, and Jacobs

C.R.

Mechanosensitivity of bone cells to oscillating fluid flow induced shear stress may be modulated by chemotransport. J Biomech, 36, 1363, 2003.

Kapur

, Bayling

D.J.

, and Lau

K.H.W.

Fluid flow shear stress stimulates human osteoblast proliferation and differentiation through multiple interacting and competing signal transduction pathways. Bone, 32, 241, 2003.

Knippenberg

, Helder

M.N.

, Doulabi

B.Z.

, Semeins

C.M.

, Wuisman

P.I.J.M.

, and Klein-Nulend

Adipose tissue-derived mesenchymal stem cells acquire bone cell-like responsiveness to fluid shear stress on osteogenic stimulation. Tissue Eng, 11, 1780, 2005.

Kreke

M.R.

, Huckle

W.R.

, and Goldstein

A.S.

Fluid flow stimulates expression of osteopontin and bone sialoprotein by bone marrow stromal cells in a temporally dependent manner. Bone, 36, 1047, 2005.

Lau

K.H.W.

, Kapur

, Kesavan

, and Baylink

D.J.

Up-regulation of the wnt, estrogen receptor, insulin-like growth factor-I, and bone morphogenetic protein pathways in C57BL/6J osteoblasts as opposed to C3H/HeJ osteoblasts in part contributes to the differential anabolic response to fluid shear. J Biol Chem, 281, 9576, 2006.

10.

Kreke

M.R.

, Sharp

L.A.

, Lee

Y.W.

, and Goldstein

A.S.

Effect of intermittent shear stress on mechanotransductive signaling and osteoblastic differentiation of bone marrow stromal cells. Tissue Eng Part A, 14, 529, 2008.

11.

Scaglione

, Wendt

, Miggino

, Papadimitropoulos

, Fato

, Quarto

, and Martin

Effects of fluid flow and calcium phosphate coating on human bone marrow stromal cells cultured in a defined 2D model system. J Biomed Mater Res A, 86A, 411, 2008.

12.

Grellier

, Bareille

, Bourget

, and Amédée

Responsiveness of human bone marrow stromal cells to shear stress. J Tissue Eng Regen Med, 3, 302, 2009.

13.

Sharp

L.A.

, Lee

, and Goldstein

A.S.

Effect of low-frequency pulsatile flow on expression of osteoblastic genes by bone marrow stromal cells. Ann Biomed Eng, 37, 445, 2009.

14.

Morris

H.L.

, Reed

C.I.

, Haycock

J.W.

, and Reilly

G.C.

Mechanisms of fluid-flow-induced matrix production in bone tissue engineering. Proc Inst Mech Eng H, 224, 1509, 2010.

15.

Gomes

M.E.

, Holtorf

H.L.

, Reis

R.L.

, and Mikos

A.G.

Influence of the porosity of starch-based fiber mesh scaffolds on the proliferation and osteogenic differentiation of bone marrow stromal cells cultured in a flow perfusion bioreactor. Tissue Eng, 12, 801, 2006.

16.

Gomes

M.E.

, Sikavitsas

V.I.

, Behravesh

, Reis

R.L.

, and Mikos

A.G.

Effect of flow perfusion on the osteogenic differentiation of bone marrow stromal cells cultured on starch-based three-dimensional scaffolds. J Biomed Mater Res A, 67A, 87, 2003.

17.

Holtorf

H.L.

, Jansen

J.A.

, and Mikos

A.G.

Flow perfusion culture induces the osteoblastic differentiation of marrow stromal cell-scaffold constructs in the absence of dexamethasone. J Biomed Mater Res A, 72A, 326, 2005.

18.

Datta

, Pham

Q.P.

, Sharma

, Sikavitsas

V.I.

, Jansen

J.A.

, and Mikos

A.G.

In vitro generated extracellular matrix and fluid shear stress synergistically enhance 3D osteoblastic differentiation. Proc Natl Acad Sci U S A, 103, 2488, 2006.

19.

Tanaka

S.M.

, Sun

H.B.

, Roeder

R.K.

, Burr

D.B.

, Turner

C.H.

, and Yokota

Osteoblast responses one hour after load-induced fluid flow in a three-dimensional porous matrix. Calcif Tissue Int, 76, 261, 2005.

20.

Cioffi

, Boschetti

, Raimondi

M.T.

, and Dubini

Modeling evaluation of the fluid-dynamic microenvironment in tissue-engineered constructs: a micro-CT based model. Biotechnol Bioeng, 93, 500, 2006.

21.

VanGordon

S.B.

, Voronov

R.S.

, Blue

T.B.

, Shambaugh

R.L.

, Papavassiliou

D.V.

, and Sikavitsas

V.I.

Effects of scaffold architecture on preosteoblastic cultures under continuous fluid shear. Ind Eng Chem Res, 50, 620, 2011.

22.

Liu

L.Y.

, Yu

, Chen

J.R.

, Tang

Z.H.

, Zong

, Shen

, Zheng

, Tong

X.M.

, Gao

C.Y.

, and Wang

J.F.

Different effects of intermittent and continuous fluid shear stresses on osteogenic differentiation of human mesenchymal stem cells. Biomech Model Mechanobiol, 11, 391, 2012.

23.

Raimondi

M.T.

, Boschetti

, Falcone

, Fiore

G.B.

, Remuzzi

, Marinoni

, Marazzi

, and Pietrabissa

Mechanobiology of engineered cartilage cultured under a quantified fluid-dynamic environment. Biomech Model Mechanobiol, 1, 69, 2002.

24.

Sikavitsas

V.I.

, Bancroft

G.N.

, Holtorf

H.L.

, Jansen

J.A.

, and Mikos

A.G.

Mineralized matrix deposition by marrow stromal osteoblasts in 3D perfusion culture increases with increasing fluid shear forces. Proc Natl Acad Sci U S A, 100, 14683, 2003.

25.

Sikavitsas

V.I.

, Bancroft

G.N.

, Lemoine

J.J.

, Liebschner

M.A.K.

, Dauner

, and Mikos

A.G.

Flow perfusion enhances the calcified matrix deposition of marrow stromal cells in biodegradable nonwoven fiber mesh scaffolds. Ann Biomed Eng, 33, 63, 2005.

26.

Zhao

, Chella

, and Ma

Effects of shear stress on 3-D human mesenchymal stem cell construct development in a perfusion bioreactor system: Experiments and hydrodynamic modelling. Biotechnol Bioeng, 96, 584, 2007.

27.

Boschetti

, Raimondi

M.T.

, Migliavacca

, and Dubini

Prediction of the micro-fluid dynamic environment imposed to three-dimensional engineered cell systems in bioreactors. J Biomech, 39, 418, 2006.

28.

Olivares

A.L.

, Marshal

, Planell

J.A.

, and Lacroix

Finite element study of scaffold architecture design and culture conditions for tissue engineering. Biomaterials, 30, 6142, 2009.

29.

Melchels

F.P.W.

, Tonnarelli

, Olivares

A.L.

, Martin

, Lacroix

, Feijen

, Wendt

D.J.

, and Grijpma

D.W.

The influence of the scaffold design on the distribution of adhering cells after perfusion cell seeding. Biomaterials, 32, 2878, 2011.

30.

Sandino

, Planell

J.A.

, and Lacroix

A finite element study of mechanical stimuli in scaffolds for bone tissue engineering. J Biomech, 41, 1005, 2008.

31.

Jungreuthmayer

, Jaasma

M.J.

, Al-Munajjed

A.A.

, Zanghellini

, Kelly

D.J.

, and O'Brien

F.J.

Deformation simulation of cells seeded on a collagen-GAG scaffold in a flow perfusion bioreactor using a sequential 3D CFD-elastostatics model. Med Eng Phys, 31, 420, 2009.

32.

Maes

, Claessens

, Moesen

, Van Oosterwyck

, Van Ransbeeck

, and Verdonck

Computational models for wall shear stress estimation in scaffolds: A comparative study of two complete geometries. J Biomech, 45, 1586, 2012.

33.

Zermatten

, Vetsch

J.R.

, Ruffoni

, Hofmann

, Müller

, and Steinfeld

Micro-computed tomography based computational fluid dynamics for the determination of shear stresses in scaffolds within a perfusion bioreactor. Ann Biomed Eng, 42, 1085, 2014.

34.

Porter

, Zauel

, Stockman

, Guldberg

, and Fyhrie

3-D computational modeling of media flow through scaffolds in a perfusion bioreactor. J Biomech, 38, 543, 2005.

35.

Voronov

R.S.

, VanGordon

S.B.

, Sikavitsas

V.I.

, and Papavassiliou

D.V.

Distribution of flow-induced stresses in highly porous media. Appl Phys Lett, 97, 024101, 2010.

36.

David

, Bonnefont-Rousselot

, Oudina

, Degat

M.C.

, Deschepper

, Viateau

, Bensidhoum

, Oddou

, and Petite

A perfusion bioreactor for engineering bone constructs: an in vitro and in vivo study. Tissue Eng Part C Methods, 17, 505, 2011.

37.

Yeatts

A.B.

, and Fisher

J.P.

Tubular perfusion system for the long-term dynamic culture of human mesenchymal stem cells. Tissue Eng Part C Methods, 17, 337, 2011.

38.

Yeatts

A.B.

, Gordon

C.N.

, and Fisher

J.P.

Formation of an aggregated alginate construct in a tubular perfusion system. Tissue Eng Part C Methods, 17, 1171, 2011.

39.

Pham

N.H.

, Voronov

R.S.

, Tummala

N.R.

, and Papavassiliou

D.V.

Bulk stress distributions in the pore space of sphere-packed beds under Darcy flow conditions. Phys Rev E, 89, 2014.

40.

Orr

D.E.

, and Burg

K.J.L.

Design of a modular bioreactor to incorporate both perfusion flow and hydrostatic compression for tissue engineering applications. Ann Biomed Eng, 36, 1228, 2008.

41.

Jean

The non-smooth contact dynamics method. Comput Methods Appl Mech Eng, 177, 235, 1999.

42.

Radjai

, and Richefeu

Contact dynamics as a nonsmooth discrete element method. Mech Mater, 41, 715, 2009.

43.

Milan

J.-L.

, Planell

J.A.

, and Lacroix

Computational modelling of the mechanical environment of osteogenesis within a polylactic acid–calcium phosphate glass scaffold. Biomaterials, 30, 4219, 2009.

44.

McCoy

R.J.

, Jungreuthmayer

, and O'Brien

F.J.

Influence of flow rate and scaffold pore size on cell behavior during mechanical stimulation in a flow perfusion bioreactor. Biotechnol Bioeng, 109, 1583, 2012.

45.

Lesman

, Blinder

, and Levenberg

Modeling of flow-induced shear stress applied on 3D cellular scaffolds: implications for vascular tissue engineering. Biotechnol Bioeng, 105, 645, 2010.

46.

Cioffi

, Küffer

, Ströbel

, Dubini

, Martin

, and Wendt

Computational evaluation of oxygen and shear stress distributions in 3D perfusion culture systems: Macro-scale and micro-structured models. J Biomech, 41, 2918, 2008.

47.

Goldstein

A.S.

, Juarez

T.M.

, Helmke

C.D.

, Gustin

M.C.

, and Mikos

A.G.

Effect of convection on osteoblastic cell growth and function in biodegradable polymer foam scaffolds. Biomaterials, 22, 1279, 2001.

48.

Jungreuthmayer

, Donahue

S.W.

, Jaasma

M.J.

, Al-Munajjed

A.A.

, Zanghellini

, Kelly

D.J.

, and O'Brien

F.J.

A comparative study of shear stresses in collagen-glycosaminoglycan and calcium phosphate scaffolds in bone tissue-engineering bioreactors. Tissue Eng Part A, 15, 1141, 2009.

49.

Jia

, Bagnaninchi

P.O.

, Yang

, El Haj

, Hinds

M.T.

, Kirkpatrick

S.J.

, and Wang

R.K.

Doppler optical coherence tomography imaging of local fluid flow and shear stress within microporous scaffolds. J Biomed Opt, 14, 034014, 2009.

50.

Pham

N.H.

, Voronov

R.S.

, VanGordon

S.B.

, Sikavitsas

V.I.

, and Papavassiliou

D.V.

Predicting the stress distribution within scaffolds with ordered architecture. Biorheology, 49, 235, 2012.

51.

Weinbaum

, Cowin

S.C.

, and Zeng

A model for the excitation of osteocytes by mechanical loading-induced bone fluid shear stresses. J Biomech, 27, 339, 1994.

52.

Klein-Nulend

, Helfrich

M.H.

, Sterck

J.G.H.

, MacPherson

, Joldersma

, Ralston

S.H.

, Semeins

C.M.

, and Burger

E.H.

Nitric oxide response to shear stress by human bone cell cultures is endothelial nitric oxide synthase dependent. Biochem Biophys Res Commun, 250, 108, 1998.

53.

Y.J.

, Batra

N.N.

, You

, Meier

S.C.

, Coe

I.A.

, Yellowley

C.E.

, and Jacobs

C.R.

Oscillatory fluid flow affects human marrow stromal cell proliferation and differentiation. J Orthop Res, 22, 1283, 2004.