A Computational Tool for the Upscaling of Regular Scaffolds During In Vitro Perfusion Culture

Abstract

Inhomogeneous and uncontrolled cellular and tissue responses in bone tissue engineering constructs, as a result of heterogeneous oxygen delivery throughout the scaffold volume, is one of the hurdles hampering clinical transfer of cell–scaffold combinations. This study presents an accurate and computationally efficient one-dimensional model that predicts the oxygen distribution for a regular cell-seeded scaffold in a perfusion bioreactor and the maximum (i.e., critical) scaffold length (L_max) as a function of given oxygen constraints. After validation against computational fluid dynamics models, the one-dimensional model was applied to calculate L_max in the perfusion direction, to ensure appropriate oxygen levels throughout the bone tissue engineering construct during in vitro culture. Both cell-related (cell density and oxygen consumption rate) and bioreactor-related (oxygen constraints and flow rate) culture parameters were varied. Results demonstrated the substantial influence of the culture parameters on L_max. In conclusion, the presented computational tool was able to predict oxygen distribution and maximum scaffold length for regular cell-seeded scaffold. It can be used to design perfusion experiments wherein quantitative knowledge on both oxygen and flow characteristics is needed.

Introduction

Bone tissue engineering (TE) can represent an alternative to autologous or allogeneic bone transplantation for the treatment of large or nonhealing bone defects. Typically, bone TE makes use of osteoprogenitor cells, which are seeded into a suitable scaffold and cultured in a (perfusion) bioreactor, resulting in a TE construct that is finally implanted into the patient's bone defect site.^1,2 The need for cell and tissue (in)growth and for transport of nutrients and waste products to and from the cells poses strict requirements on the scaffold architecture, in terms of pore size, porosity, interconnectivity, and permeability—parameters that are obviously interdependent.^3,4 As to the in vitro culture, bioreactors should provide a suitable environment, which includes efficient nutrient delivery, waste removal, and mechanical stimulation.^5–7

Studies performed in vitro have shown that inhomogeneous delivery of oxygen throughout the TE construct volume can lead to inhomogeneous cellular and tissue responses. Local reduction of cell viability has been attributed to too low levels of oxygen tension, whereas increased cell density at the scaffold borders have been related to high levels of oxygen tension.^8–11 These phenomena can be expected to become even more important with increased bone defect size, thus rendering the so-called upscaling of tissue engineering constructs to clinically relevant sizes a major challenge.^12,13 Oxygen availability has been considered as a limiting factor for scaffold upscaling,¹⁴ because of the low solubility and diffusivity of oxygen in culture medium, opposed to high cellular oxygen consumption.¹⁵ Apart from its effect on cell viability, other studies have demonstrated that cell proliferation and differentiation are influenced by oxygen tension as well.^8,16,17 Limiting ourselves to skeletal applications, low oxygen tension (oxygen levels around 5%) has been found to enhance the proliferation of progenitor cells and to stimulate chondrogenesis.^18–21 Conflicting results have been reported with respect to osteogenic differentiation: some studies found enhanced osteogenesis at low oxygen tension, but this was not confirmed by others.^22–24

All these studies suggest that oxygen tension is a crucial environmental signal that needs to be controlled during in vitro culturing. Perfusion bioreactors can be used to influence and control the oxygen tension inside TE constructs and to ensure adequate oxygen delivery by advective transport. At the same time, construct perfusion will also influence fluid flow-induced shear stresses, another environmental signal that has been shown to influence cell behavior.²⁵ Focusing here on the 3D culture of osteogenic cells in scaffolds, both proliferation and osteogenic differentiation can be enhanced by fluid flow, but conflicting results have been reported with respect to the corresponding values of wall shear stress (WSS).^26,27 Excessive values in the order of 50 mPa have been related to increased cell apoptosis.²⁸

Quantification of both oxygen tension and WSS is crucial when studying the role of scaffold perfusion in tissue engineering. Different approaches have been taken to quantitatively study both the oxygen and the WSS distribution inside porous scaffolds. Regarding the quantification of the shear stress either analytical formulae or computational fluid dynamics (CFD) models were used.^29–31 Experimental techniques such as microparticle imaging velocimetry (μ-PIV) can be used as well, either or not in combination with CFD calculations.³²

Similarly, computational models have been used extensively to quantify the oxygen distribution inside a perfused, cell-seeded scaffold.^33–35 Cioffi et al. compared two CFD modeling approaches to predict the oxygen tension inside a perfused, irregular scaffold for cartilage engineering.³⁵ On the one hand, a microscale model was created, which considered the heterogeneous flow through the microarchitecture. On the other hand, a macroscale model was created as well, considering only the permeability and the porous medium flow through the scaffold. They demonstrated that the macroscale model was sufficiently accurate to predict the oxygen gradient along the flow direction, compared with a computationally more expensive microscale model. It must be noticed that in the microscale model the cells that consume oxygen were assumed to be homogeneously distributed throughout the fluid volume (pore volume). In reality, cells are attached to the scaffold surface, leading to a heterogeneous cell distribution and additional oxygen gradients, perpendicular to the scaffold surface. These gradients may affect the comparison between a microscale and macroscale modeling approach and, therefore, the validity of the macroscale approach. Further, the scaffold considered in this study was an irregular scaffold for which the geometry was reconstructed using microfocus computed tomography (μCT) images. The large computational cost, associated with a microscale model, did not allow modeling the entire scaffold. Hence, only part of it that served as a region of interest (representative volume) was considered and symmetry boundary conditions were assumed on the lateral walls of the region of interest. Maes et al., in another μCT-based CFD study of flow in irregular scaffolds, demonstrated that the size of the considered region of interest as well as the presence of the boundary conditions on the lateral walls can strongly influence the simulation results.³⁶ These results warrant for a careful selection of a region of interest and interpretation of model results, especially for strongly heterogeneous scaffolds.

In this context, it is easy to understand the advantage of using designed, regular scaffolds that are based on a repeatable unit cell geometry, in terms of both computational cost and model accuracy. Owing to the symmetry of the scaffold structure, the flow field will be symmetrical as well. The definition of a smaller region of interest becomes straightforward and the assumption of symmetry boundary conditions is fully valid.

The goals of this study are twofold. First, an accurate but, at the same time, computationally efficient model will be determined to predict the oxygen distribution for a regular cell-seeded scaffold in a perfusion bioreactor for a given cell type and density (i.e., snapshot of the cell-seeded scaffold, not taking into account any temporal changes in the cell density). Related to this first goal, we hypothesize that a computationally efficient one-dimensional (1D) model is sufficiently accurate to capture the oxygen distribution. Second, the 1D model will be applied to predict the maximum (i.e., critical) scaffold length as a function of given oxygen constraints. Both goals will be elaborated and illustrated by making use of a regular cell-seeded Ti6Al4V scaffold in an in-house–developed perfusion bioreactor.³²

Materials and Methods

Experimental setup

The experimental test case in this study was a one unit cell high version of an additive manufacturing (AM) Ti6Al4V scaffold produced by Selective Laser Sintering (SLS). This scaffold represents a simplification of the full AM 3D version and its design allows for a reliable representation of the flow field by considering just one unit cell. The scaffold overall dimensions are 6×10×0.5 mm³. The geometric features of the scaffold are shown in Figure 1 and consist of a regular structure with 0.5 mm by 0.5 mm struts and 0.3-mm-wide channels at a 45° angle with respect to the flow direction. The construct is perfused in a so-called 2D+ bioreactor chamber, which permits the real-time visualization of the fluid flow through the channels of the scaffold. A detailed description of this setup can be found in ref.³²

FIG. 1.

Perfusion 2D+ bioreactor (left) with part of the scaffold design and an enlargement of the unit cell (top right). The red zone represents the modeled region of interest. The white arrows represent the flow direction (left). The geometrical features of the scaffold and the modeled region of interest are summarized in the table (bottom right). Color images available online at www.liebertonline.com/tec

Computational models

Model description

Four different models (further denoted as model A to D) were developed to describe the medium flow and oxygen transport in a regular, cell-seeded scaffold. A description of the model characteristics in terms of the modeling of medium flow, oxygen transport (including oxygen consumption by the cells present in the scaffold), and the distribution of the cells inside the scaffold is given as follows. Each model has its own level of complexity, with model A being the most complex (and therefore most realistic with respect to the experimental conditions) and model D the simplest.

Three of the four models were CFD models of the culture medium flow inside the 2D+ scaffold: two detailed models (models A and B), both solving the Navier-Stokes equation, and one porous model (model C), solving Darcy's law.

The ACIS-based solid modeler Gambit 2.2 (Fluent) was used to build the 3D geometry and the mesh. The finite volume code Fluent 6.3 (Fluent) was used to set up and solve the CFD problem. Because of symmetry and repeatability of the scaffold geometry and flow field, the width and the length of the domain were limited to one and three unit cells, respectively. The fluid domain, corresponding to the region of interest in Figure 1, is depicted in Figure 2 and measured 5.3 mm in length (2 mm for the inlet region and 3.3 mm for the three unit cells), 1 mm in width, and 0.5 mm in height. The culture medium was modeled as an incompressible and homogeneous Newtonian fluid with a viscosity of 10⁻³ Pa s and a density of 10³ kg/m³.

FIG. 2.

Geometry representation of the modeled fluid domain and boundary conditions. The white arrows represent the flow direction.

Detailed models (A and B): The detailed models A and B were built to evaluate the influence of the local microstructure on the velocity field and consequently on the oxygen distribution. The actual scaffold geometry was considered and the flow was described by means of the steady-state Navier-Stokes equation: \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*} \rho {\bf v} \cdot \nabla {\bf v} = - \nabla p + \mu \nabla^2 {\bf v} \tag{1} \end{align*} \end{document}

where ρ and μ are, respectively, the medium density and viscosity; p and v are the pressure and the velocity vector, respectively. A mesh refinement study was conducted to determine the appropriate number of finite volumes, resulting in a mesh of 1,137,708 tetrahedral volumes. The average tetrahedral edge length in the fluid domain was 25 μm, with a progressive refinement down to an edge length of 10 μm close to the scaffold walls. The boundary conditions for the fluid dynamic analysis consisted of a constant inlet velocity (uniform velocity field), an atmospheric pressure (set to zero) applied to the outlet, and no-slip condition on the scaffold walls. Symmetry conditions were applied where appropriate (Fig. 2).

Porous model (C): The fluid-filled scaffold was considered as a porous continuous medium with a permeability K of 1.84×10⁻⁹ m², as calculated by means of the detailed model A, and a porosity Φ of 61% (ratio between fluid volume and total volume). The porous medium flow was described by means of Darcy's law: \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*} {\bf v}_m = - \frac {K} {\mu} \nabla p \tag {2} \end{align*} \end{document}

where v_m is the mean velocity (averaged over cross-sectional area).

To compare this model with the detailed models described in the previous section, the same symmetric fluid domain was taken into account and the same boundary conditions were assigned. The mesh presented 182,384 tetrahedral volumes and the average tetrahedral edge length was 50 μm.

For both the detailed (A and B) and porous (C) models, the oxygen distribution inside the domain was obtained by solving the steady-state conservation law, taking into account advection, diffusion, and consumption: \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*} \rho {\bf v} \cdot \nabla c = \rho D \nabla^2 c - r \tag{3} \end{align*} \end{document}

where c is the oxygen concentration, expressed in terms of mass fraction (ratio between oxygen mass and fluid mass); D is the diffusion coefficient of oxygen (equal to 3×10⁻⁹ m²/s when assuming water as diffusion medium); and r is the volumetric oxygen consumption rate (in kg/[s m³]), which is assumed to follow a Michaelis-Menten (MM) kinetic law: \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*} r = r_{\max} \frac {c} {k_m + c} \tag {4} \end{align*} \end{document}

In the equation (4), r_max is the maximum volumetric consumption rate and k_m corresponds to the oxygen concentration at which the consumption rate is equal to half of r_max. For a given cellular maximum consumption rate r_max,cell, the maximum volumetric consumption rate r_max is given by \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*} r_{\max} = r_{\max , cell} \frac {N} {V} \tag {5} \end{align*} \end{document}

with N being the total number of seeded cells and V the fluid volume to which cells (and therefore consumption) are assigned.

1D model (D): A further model simplification was developed to capture the oxygen distribution along the length of the scaffold (i.e., along the flow direction). In this case, a 1D mass conservation law was considered, assuming a constant (mean) velocity and neglecting the diffusive phenomena: \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*} \rho {\bf v} \frac {dc} {dx} = r_{\max} \frac {c} {k_m + c} \tag {6} \end{align*} \end{document}

where x is the coordinate relative to the scaffold length (i.e., in the flow direction). The Matlab 7.1 code was used to solve this nonlinear equation by means of the finite difference method (FDM). A uniform grid with an element size of 10 μm was created and a backward finite difference scheme having a first-order approximation was adopted.

Cell distribution: Figure 3 shows how, for the three CFD models (A to C), the same number of cells was distributed into the fluid domain. In model A, the cells were considered to be attached (as a cell layer) to the scaffold walls and to the bottom walls of the channels. This was implemented by only assigning a consumption term in a 10-μm-thick fluid layer at the wall surfaces, resulting in a nonuniform oxygen consumption throughout the fluid domain (zero consumption outside the 10-μm-thick layer). In model B, the cells are assumed to be uniformly distributed throughout the entire fluid volume (equal in this case to the scaffold pore volume), similar to that in a previous study,³⁵ thus resulting in uniform oxygen consumption. Finally, as model C and also model D do not take into account the detailed architecture of the scaffold (and, therefore, do not make a distinction between scaffold pore volume and strut volume), the cells are assumed to occupy the entire region of interest (i.e., the red area in Fig. 1). The values for the volume to which cellular consumption was assigned [V in equation (5)] for the different models are summarized in Table 1.

FIG. 3.

Geometry of the scaffold unit cell with representation of the distribution of the oxygen-consuming cells for the detailed (A, B) and porous (C) computational fluid dynamics (CFD) models. Distributions of a given cell number are given in a cross-sectional plane parallel (left figures) and perpendicular (right figures) to the bioreactor bottom plane. For models A and B, fluid is represented in white and scaffold strut in gray. Cross-sectional planes are indicated in the top figure (blue=parallel; red=perpendicular). Color images available online at www.liebertonline.com/tec

Table 1.

Summary of Oxygen Consumption Characteristics for the Four Models Used During Model Evaluation

Model name	Model type	Cell distribution	V (mm ³ )	r_max (mol/s×mm ³ )
A	Detailed	Localized on the walls	8.41×10⁻²	2.8×10⁻¹¹
B	Detailed	Homogenous	1.21	1.9×10⁻¹²
C	Porous	Homogenous	1.91	1.2×10⁻¹²
D	ID	Homogenous	1.91	1.2×10⁻¹²

The value of r_max for the different models was calculated according to equation (5), with N=62×10³ and r_max,cell=3.8×10⁻¹⁷ mol/cell (ES cells³⁷).

Model evaluation

The accuracy in terms of the predicted distribution of oxygen tension pO₂ was evaluated for each model, with model A being the reference model as it represents best the actual scaffold geometry and cell distribution. As an accurate prediction can be expected to be most challenging when strong oxygen gradients are encountered, culture conditions that favor such gradients were selected. This was achieved by selecting high maximum volumetric consumption rates—which in turn was obtained by choosing a high cellular consumption rate and a high cell density—and low bioreactor inlet velocities v_in. Limiting ourselves to cell sources for skeletal engineering applications, cellular consumption rates have been reported in a range from 4×10⁻¹⁸ mol/(cell s) (lower bound for chondrocytes¹⁰) to 3.8×10⁻¹⁷ mol/(cell s) (embryonic stem [ES] cells³⁷). The latter value was adopted here. The number of seeded cells into the modeled region of interest (red area in Fig. 1) was 62×10³, which corresponded to a relatively high cell-seeding density ρ_cell (with respect to total volume, i.e., sum of fluid and strut volume) of 32×10⁶ cells/mL. The resulting maximum volumetric consumption rates r_max for the different models are summarized in Table 1.

For v_in, values of 0.03 and 0.1 mm/s were evaluated, the former being of the same order of magnitude as the lowest values found in literature for skeletal engineering applications.^35,38,39 For our 2D+ bioreactor setup, this leads to flow rates of 0.0054 and 0.018 mL/min, respectively.

Apart from strong oxygen gradients, it was important to assess model accuracy over a broad range of oxygen values, so that the effect of the MM nonlinear consumption rate, embedded in the MM kinetic law, could be determined. To this end, it is important that pO₂ values in the range of k_m are encountered. The value adopted for k_m corresponded to an oxygen tension of 0.68%, a value previously used for ES cells.³⁷ The oxygen tension at the bioreactor inlet (pO_2in) was taken equal to 5%. In combination with the selected volumetric consumption rate and the lowest inlet velocity (0.03 mm/s), it was verified that these conditions led to complete oxygen depletion at the outlet of the modeled region of interest (pO_2out=0) (as will be shown in the Results section).

Model application: Scaffold upscaling

Once the model with the lowest computational cost that was able to accurately predict the oxygen tension along the flow direction was selected, it was applied to scaffold upscaling. To this end, the maximum (i.e., critical) scaffold length L_max, which ensured appropriate oxygen tension levels throughout the cell-seeded scaffold, was determined as a function of both cell- and bioreactor-related parameters. Cell-related parameters included cell density and oxygen consumption rate, and bioreactor-related parameters concerned the oxygen constraints and bioreactor inlet velocity. The selected parameter values are summarized in Table 2. Cellular consumption rates corresponded to the values that were mentioned before for chondrocytes (low consumption rate¹⁰) and ES cells (high consumption rate³⁷). As to the oxygen constraints, two target ranges of oxygen tension were selected: one corresponding to low values ([pO_2in, pO_2out]=[5%, 1%]) and one corresponding to high values ([pO_2in, pO_2out]=[21%, 10%]). Low values may be considered as an oxygen regime to promote the proliferation of progenitor cells (as mentioned earlier, see, e.g., ref.¹⁸), whereas high values may (or may not) enhance osteogenic differentiation (see, e.g., ref.²⁴). As the goal is to merely illustrate the potential of the model for a number of selected oxygen constraints, the biological effect of these oxygen levels is not important in this setting.

Table 2.

Parameter Values (or Range) of Cell- and Bioreactor-Related Parameters Considered for Scaffold Upscaling

r_max,cell (mol/cell)		v _in (mm/s)			(pO_2in, pO_2out) (%, %)
L	H	L	H	ρ_cell (10⁶ cells/mL)	L	H
4×10⁻¹⁸	3.8×10⁻¹⁷	0.03	0.1	1–32	(5, 1)	(21, 10)

L and H refer to low and high values, respectively.

The bioreactor inlet velocity v_in influences not only oxygen transport, but also WSS values in the scaffold. Given the fact that excessive WSS may lead to cell death (as mentioned in the Introduction section), an upper limit for the inlet velocity, corresponding to a critical value of WSS, must be determined. Therefore, WSS distribution and the relation between WSS values and inlet velocity were quantified by means of a detailed CFD model (A and B), to determine an upper limit for v_in. This was done prior to addressing scaffold upscaling. To this end, three different values for v_in were considered: 0.03, 0.1, and 1 mm/s.

Results

Model evaluation

Oxygen distribution was evaluated for all four models in both the longitudinal (along the flow direction, x direction in Fig. 4A) and transverse (perpendicular to flow direction, y direction in Fig. 4A) directions. As to the longitudinal direction, several cross sections were considered perpendicular to it and for each cross section the average oxygen tension was calculated (Fig. 4B, D). For an inlet velocity of 0.03 mm/s, all four models predicted the same trend for the oxygen tension in the longitudinal direction: a linear decrease in the first part of the scaffold followed by a logarithmic decrease, leading to complete oxygen depletion at the outlet (Fig. 4B). In the transverse direction, models A and B presented a nonuniform oxygen distribution, whereas models C and D led to a constant oxygen tension (Fig. 4C; evaluated along a line, perpendicular to the longitudinal direction, at x=1 mm and z=0 mm). For the selected line, the absolute difference in terms of oxygen tension was smaller than 0.2% oxygen between models A and B, whereas it was smaller than 0.3% and 0.4%, respectively, between models A and C and between models A and D.

FIG. 4.

Oxygen distribution in a 2D+ bioreactor for a cell seeding density ρ_cell of 32×10⁶ cells/mL, a maximum cellular consumption rate r_max,cell of 3.8×10⁻¹⁷ mol/(cell s) (ES cells³⁷) and an inlet oxygen tension of 5%. Oxygen values were calculated by means of two detailed CFD models (models A and B), a porous CFD model (model C), and a one-dimensional model (model D). (A) Geometry of models A and B showing the longitudinal (x) and transverse (y) directions and the position of the line that represents the oxygen distribution in the transverse direction. (B, D) Average oxygen tension as a function of the longitudinal distance (for an inlet velocity of 0.03 mm/s (B) and 0.1 mm/s (D), respectively). (C, E) Oxygen tension as a function of the transverse distance in a line at x=1 mm and z=0 mm for an inlet velocity of 0.03 mm/s (C) and 0.1 mm/s (E), respectively. Color images available online at www.liebertonline.com/tec

For an inlet velocity of 0.1 mm/s, the agreement between the different models was confirmed, as shown in Figure 4D. In this case, a linear decrease of oxygen tension as a function of the longitudinal distance was predicted for the entire modeled region. The absolute difference in terms of oxygen tension, evaluated in the line selected in the transverse direction (Fig. 4A), was smaller than 0.09% between models A and B, whereas it was smaller than 0.13% and 0.15%, respectively, between models A and C and between models A and D.

As there were small differences in oxygen tension between all models, model D was preferred for studying scaffold upscaling, because of its low computational cost.

Model application: Scaffold upscaling

Prior to the application of model D to scaffold upscaling, the WSS distribution was studied to determine an upper limit for the inlet velocity v_in. Figure 5 depicts the distribution in the region of interest for an inlet velocity of 0.1 mm/s, obtained by means of the detailed CFD model (models A and B, which produce identical flow results): the maximum value of the WSS was found at the corners of the struts. The WSS distribution was identical for the other two inlet velocities considered here (0.03 and 1 mm/s) as a consequence of the linear relation that exists between the WSS and the inlet velocity for a laminar flow (Reynolds number Re<1 for the range of inlet velocities considered). For the different values of inlet velocity, the mean and maximum values of WSS are summarized in Table 3. An inlet velocity of 1 mm/s resulted in a WSS mean value of 34 mPa, which is in the range of WSS supposed to generate cell apoptosis for osteoblast-like cells, as previously mentioned.²⁸ Hence, this inlet velocity was discarded, and only the values of 0.03 and 0.1 mm/s were further considered as case studies for scaffold upscaling (see also Table 2).

FIG. 5.

Wall shear stress contour plot in the modeled scaffold region for the detailed models A and B for an inlet velocity of 0.1 mm/s. Color images available online at www.liebertonline.com/tec

Table 3.

Wall Shear Stress Results of the Computational Fluid Dynamics Simulation Using the Detailed Model (A and B)

v _in (mm/s)	τ _max (mPa)	τ _mean (mPa)
0.03	6	1.02
0.1	20	3.4
1	200	34

The results for the maximum scaffold length L_max, as iteratively predicted by model D for different culture parameters, are shown in Figure 6. It can be noticed that L_max is inversely proportional to the cell-seeding density. For instance, for the oxygen target range of [5%, 1%], its value decreases from 48.9 mm for 10⁶ cells/mL to 1.5 mm for 32×10⁶ cells/mL in case of an inlet velocity of 0.03 mm/s. Moreover, for this oxygen range and a cell density of 10×10⁶ cells/mL, increasing the inlet velocity v_in from 0.03 to 0.1 mm/s leads to an increase of L_max, from 4.9 to 16.3 mm. Finally, when increasing the length (pO_2in−pO_2out) of the oxygen target range from [5%, 1%] to [21%, 10%], L_max also increases: for a cell density of 10×10⁶ cells/mL and an inlet velocity of 0.03 mm/s, L_max changes from 4.9 to 11 mm.

FIG. 6.

Maximum scaffold length L_max as a function of cell seeding density ρ_cell, as predicted by model D. Results are shown for low ([pO_2in, pO_2out]=[5%, 1%]) (A, B) and high ([pO_2in, pO_2out]=[21%, 10%]) (C, D) target ranges of oxygen tension and an inlet velocity equal to 0.1 mm/s (A, C) and 0.03 mm/s (B, D), respectively. Maximum cellular consumption rate r_max,cell was equal to 3.8×10⁻¹⁷ mol/(cell s) (ES cells³⁷).

The effect of combined changes of all culture parameters on L_max is summarized in Table 4. The high (H) and low (L) values for the inlet velocity v_in, the oxygen target range [pO_2in, pO_2out], and the cellular consumption rate r_max,cell correspond to values in Table 2. As an example, three cell-seeding densities ρ_cell were considered: 1×10⁶, 10×10⁶, and 30×10⁶ cells/mL. The value of L_max gradually increases from the top left corner (having the lowest values equal to 1.6 mm) to the bottom right corner (with the highest value equal to 3448 mm).

Table 4.

L _max Results of One-Dimensional Model (Model D), Expressed in mm, as a Function of Culture Parameters

	v _in−[pO_2in, pO_2out]−r_max,cell
ρ_cell (10⁶ cell/mL)	L - L - H	L - H - H	H - L - H	H - H - H	L - L - L	L - H - L	H - L - L	H - H - L
30	1.6	3.6	5.4	12.2	15.5	34.6	51.7	116.3
10	4.9	11	16.3	36.5	46.6	104	155.2	348
1	49	110	163.2	365	466	1040	1552	3448

Low (L) and high (H) values for v_in, [pO_2in, pO_2out], and r_max,cell correspond to values in Table 2. The cell density ρ_cell equaled 1×10⁶, 10×10⁶, and 30×10⁶ cells/mL.

Discussion

The first part of this study was focused on the determination of the most efficient computational model (in terms of computational cost) that is still sufficiently accurate in capturing the oxygen tension inside a regular cell-seeded scaffold, placed in a perfusion bioreactor. As was evident from Figure 4B and D, all models predicted very similar results for the oxygen gradients in the longitudinal (i.e., flow) direction, with the maximum error between the most complex (detailed CFD model A) and the simplest model (1D model D) being only 0.4% in terms of average oxygen tension. For the oxygen gradients in the transverse direction, the models behaved differently (Fig. 4C, E). Owing to its porous medium description, model C had a constant velocity (uniform velocity field) and a uniform oxygen consumption in each cross-section perpendicular to the flow direction. As a result, oxygen tension was constant in a transverse direction. As model D only took into account the change of oxygen tension in the longitudinal direction [see also equation (6)], it inherently led to a constant value of oxygen tension in a transverse direction as well. On the contrary, for models A and B, a nonuniform oxygen distribution was found in the transverse direction, owing to the nonuniform cell distribution (for model A) and the nonuniform velocity field (for models A and B). Despite this qualitatively different behavior, the absolute differences between models A and D remained low in the transverse direction as well, with the maximum error being located in the regions close to the scaffold walls and being less than 0.4% of oxygen tension. As it is unlikely that such small errors in oxygen tension would be relevant for predicting a biological outcome, it was concluded that model D is the most appropriate (computationally most efficient) model for the application envisioned here, which is the prediction of a maximum scaffold length for a given target range of oxygen tension and cell type. Moreover, given the fact that model comparison was performed for conditions that led to strong (local) oxygen gradients (due to high cell densities and cellular consumption rates) and for oxygen values that covered the entire range of the MM kinetic law (including the nonlinear part), it strengthens the applicability of model D to these more challenging (from a computational point of view) culture conditions. These results confirm the research hypothesis on the accuracy of the 1D model.

The nonlinear consumption rate, expressed by the MM kinetic law, becomes important for oxygen tensions in the range of k_m (which can also be appreciated from Fig. 4). This can be clearly seen in Figure 4B, wherein a nonlinear oxygen gradient was found in the longitudinal direction for the lowest (0.03 mm/s) inlet velocity. In contrast, for the higher (0.1 mm/s) inlet velocity, oxygen tensions were well above k_m, leading to a linear oxygen gradient (Fig. 4D). These results also imply that for sufficiently large oxygen tension (pO₂≫ k_m), equation (4) can be replaced by a constant consumption rate (which is equal to the maximum consumption rate r_max). As a consequence, the 1D mass conservation law of model D [equation (6)] can be further simplified and solved analytically. Therefore, for pO₂≫ k_m, a very simple expression for L_max can be written: \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*} L_{\max} = \frac {{\bf v}_{in} (pO_{2in} - pO_{2out})} {r_{\max , cell} \rho_{cell}} \tag {7} \end{align*} \end{document}

In the second part of this study, model D was applied as a computational tool to tackle scaffold upscaling and to calculate the maximum scaffold length for a given target range of oxygen tension and as a function of cell- and bioreactor-related culture parameters. The approach was illustrated for both high (10%<pO₂<21%) and low (1%<pO₂<5%) oxygen ranges. Although equation (7) is sufficiently accurate for the calculation of L_max for the high oxygen range, it underestimates L_max of about 13% for the low oxygen range. Nevertheless, it is instructive to qualitatively explain the results in Figure 6 and Table 4 by means of this equation. Basically, this equation shows that L_max is proportional to the bioreactor inlet velocity and the oxygen target range Δp_O2 (which equals [pO_2in−pO_2out]) and inversely proportional to the cellular consumption rate and cell density.

In general, one can notice from the results in Table 4 the substantial effect of culture parameters on L_max, which varies over 3 orders of magnitude (from millimeter to meter scale). Obviously, values of L_max of the order of meters are not realistic from a production or clinical point of view. In such cases, the model simply points out that oxygen availability inside the scaffold is not the limiting factor for upscaling. This is, for example, the case for the lowest cell density (10⁶ cells/mL), where even for the most demanding culture conditions (low inlet velocity and oxygen range, and high consumption rate) L_max amounts to 48.9 mm. The values for L_max changed considerably between the low and high oxygen range: for an inlet velocity of 0.1 mm/s and for a seeding density of 30×10⁶ ES cells/mL, L_max equaled 5.44 and 12.16 mm, respectively, for the low and higher oxygen range. This effect follows from equation (7), wherein an increase of the length of the oxygen target range Δp_O2 (from 4% to 11%) leads to an increase of L_max. However, as equation (7) is not valid for the low target range of oxygen tension (due to the nonlinearity of the Michaelis-Menten law for low oxygen tension), it leads to an underestimation of L_max with 13% for this range.

As mentioned earlier, the maximum specific cell consumption rate r_max,cell is inversely proportional to L_max. However, in contrast to, for example, the bioreactor inlet velocity, the value of r_max,cell cannot be varied to obtain the required L_max, as it is intrinsic to a certain cell (pheno)type. Therefore, once the cell of interest is selected, the value of r_max,cell becomes a constant in the analysis. On the other hand, it is interesting to notice how the value of L_max changes with cell type. As an example, Table 4 shows that L_max values were about 10 times higher when chondrocytes were used instead of ES cells, which follows from the fact that r_max,cell is about 10 times lower for chondrocytes.

Finally, and again following equation (7), L_max increases linearly with an increase of the inlet velocity v_in. As this is a bioreactor-related parameter, it can in principle be changed to reach a desired L_max, thus improving the advective transport. An example can be found in Figure 6: for a scaffold seeded with 1×10⁶ ES cells/mL and for the low oxygen target range, L_max linearly increased from 48.96 to 163.2 mm when the inlet velocity increased from 0.03 to 0.1 mm/s (Fig. 6A, B). As mentioned earlier, it cannot be increased limitless because excessive inlet velocities will lead to WSS values that are detrimental for cell viability. A compromise between the WSS and oxygen transport for the determination of the appropriate flow rate must be made. To do so, it is essential to assess an upper limit for the WSS value (which should follow from biological experiments), which is then translated into an upper limit for the inlet velocity. To this end, the detailed CFD model (model A or B), able to capture the 3D velocity gradients and consequently the WSS values, was applied. The results in Table 3 illustrate that, because of laminar flow, the relation between the WSS and the inlet velocity is linear. Further, for regular scaffolds that are composed as a period repetition of a unit cell, the WSS distribution will be periodic as well, as shown in Figure 5. Thus, in principle, a detailed CFD model of the fluid flow in one unit cell gives the necessary information on the WSS distribution and, in combination with the laminar flow regime, provides an upper limit to the inlet velocity. For irregular structures, because of the absence of periodicity and increased heterogeneity, a similar approach for the assessment of the WSS distribution will not be possible, leading to much more expensive computations.^35,36

In summary, for a regular scaffold, based on a repeatable unit cell, an expensive CFD analysis, performed on the entire scaffold, can be fully replaced by a CFD study of the fluid flow in one unit cell (for the determination of the WSS distribution) and a 1D model for the prediction of the oxygen distribution and the determination of L_max. Figure 7 summarizes the different steps of the approach. A CFD model of the representative unit cell is performed by imposing an arbitrary (but reasonable, i.e., sufficiently low to obtain a laminar flow) value for the inlet velocity and the WSS distribution is determined. Experimental input on the effect of WSS on cell behavior is required to have an upper limit for the WSS, which is not detrimental (or beneficial) for cell behavior. Consequently, because of the proportionality of the WSS to the inlet velocity, an upper limit v_max for v_in is established. Constraints on the oxygen range come from biological experiments: depending on the cell type and on the desired cell response, a desired range of oxygen values can be fixed. For the desired oxygen range and within the valid range of v_in (v_in<v_max), L_max can be determined as a function of the cell-related parameters (cell density and cell consumption rate).

FIG. 7.

Flow chart of the computational approach for the determination of L_max as a combination of a detailed CFD analysis of a unit cell and a one-dimensional analysis of oxygen transport.

Alternative applications of the proposed computational approach can be thought of, depending on the problem definition (i.e., what are the dependent and independent variables). Figure 8 illustrates how for a given scaffold length, oxygen target range, and oxygen consumption rate, the minimal inlet velocity to meet these requirements can be calculated. If the inlet velocity exceeds the upper limit v_max (as calculated by means of a CFD model of a unit cell, as previously described), then a feedback loop to the cell seeding density is necessary. The aim of this loop is to decrease proportionally the number of seeded cells so that the oxygen constraints can be met for an acceptable value for the inlet velocity (i.e., v_in<v_max).

FIG. 8.

Flow chart of the computational approach for the determination of the inlet velocity for a given scaffold length L.

The presented computational methodology relies on a number of assumptions, which need to be discussed to better appreciate its applicability.

The study was performed on a thin (0.5-mm-thick) scaffold, which was dictated by the previously developed (2D+) bioreactor. Although it limits the use of this setup to laboratory applications, it does not limit the applicability of the computational approach. As long as the scaffold architecture is regular and periodic, the approach can be extended to several true 3D scaffolds geometries (i.e., thickness of the same order as length and width).

In the detailed model A, cells were assumed to be attached to the scaffold surface, which is more realistic than a homogeneous distribution throughout the entire pore volume (model B, which is similar to that described by Cioffi et al.³⁵). As it was found that the cell distribution only had a minor effect on the oxygen distribution (which followed from the comparison between models A, B, and C, the latter being a porous model), our study confirms and extends the findings of Cioffi et al. with respect to the applicability of a porous (macroscale) model. Model A assumes that the cell surface density (number of cells per unit of scaffold area) is constant. Real cell distributions, as obtained after cell seeding, may be heterogeneous either at a macro- or microscale and may involve phenomena such as cell clustering, such as described by Wendt et al.⁴⁰ At the same time, this study demonstrated that the cell distribution strongly depends on the seeding protocol and that perfusion seeding may result in fairly homogeneous cell distributions, both at the macroscale (entire scaffold) and the microscale (one scaffold pore). Given the fact that these results were obtained with irregular scaffolds and that one can expect an even higher cell uniformity for regular scaffolds, it supports the assumption made for detailed model A.

Nevertheless, in case of a strongly heterogeneous cell distribution, higher oxygen gradients may develop through the scaffold because of highly localized consumption. As a result, diffusive effects may start to play a role as well, rendering model D invalid. Therefore, the validity of a uniform cell surface distribution must be checked experimentally.

The cell dynamics (change of cell density due to, e.g., cell proliferation or death) was not considered in the study. Instead, a static snapshot of the cell density was considered, which can be interpreted as the initial seeding density. This also means that the model cannot be used to predict the long-term (days, weeks) oxygen dynamics during culturing. The modeled cell density may therefore be regarded as the initial seeding density. Only if, during culturing, the effect of cell dynamics and extracellular matrix synthesis on the scaffold pore volume would be negligible (see also discussion below), the modeling approach may provide a static snapshot of the oxygen distribution at an arbitrary time point (provided one knows the cell density for that time point). To verify that the assumption on constant cell density was valid for the goals set out in this study (and particularly the accurate calculation of oxygen distribution for a given time point), a transient mass transport simulation was performed for model A (reference model). The time to reach the steady state for the oxygen distribution was compared with typical cell duplication times. The oxygen distribution inside the region of interest of the scaffold reached a steady state after 80 s. As cell duplication times can vary from 21 to 36 h for adult chondrocytes⁴¹ and ES cells,⁴² respectively, it confirms the validity of the assumption on constant cell density.

In addition, in the detailed model A, the thickness of the oxygen-consuming cell layer near the scaffold surface was kept constant and equal to 10 μm, which corresponds to a volume fraction (with respect to the pore volume) of 7%. This layer only plays a role for the calculation of the oxygen transport. For the flow calculation, it is considered as part of the fluid domain, and therefore, the cell volume is neglected. One can estimate the total cell volume as follows: assuming a value of 850 μm³ for the volume of a single (spherical) cell,⁴³ it leads to a total cell volume fraction (with respect to the pore volume) of 0.1%–4.5% for the considered range of cell densities of 10⁶ to 32×10⁶ cells/mL. As these values are small, neglecting the cell volume seems valid, even for large cell densities.

Apart from cell dynamics, tissue dynamics was neglected as well. In case of extensive extracellular matrix synthesis, this assumption may no longer be valid. The approach can be extended by modeling the tissue volume (and if necessary also the cell volume) more correctly. At different time points, more narrow channels (for the detailed models A and B) and lower permeability (for porous model C) can be implemented. This will lead to larger pressures and WSS for the same inlet velocity, therefore affecting the upper limit for the inlet velocity. Further, for thick layers of cells and tissues, the diffusion through these layers may be the dominant transport mechanism, which must be taken into account for an accurate prediction of the local oxygen values.

The approach developed in this study focuses on the quantification of important environmental signals, such as oxygen tension and WSS, for a given cell density. As mentioned in the Introduction section, these signals are important mediators of cell viability, proliferation, differentiation, and matrix synthesis. As our computational approach was developed to address the scaffold upscaling for a given time point, no attempt was made here to perform dynamic simulations of the interaction between these signals and cell and tissue growth. Examples of such simulation studies can be found in the literature, wherein cell and tissue growth (and differentiation) in tissue engineering constructs are made dependent on local oxygen tension^16,44 or mechanical stimuli, such as stress, strain, or fluid flow.^45–47 In the future, our computational approach could be embedded in a larger simulation platform that takes into account cell and tissue dynamics and their regulation by means of local WSS and oxygen tension.

Conclusion

In conclusion, a computationally efficient methodology that enables to predict the oxygen distribution inside perfused, regular scaffolds has been developed. Further, by combining a detailed CFD model of a scaffold unit cell and a 1D model of oxygen transport, the maximum scaffold length can be calculated as a function of given oxygen and WSS constraints. The computational tool can be used to design perfusion experiments wherein quantitative knowledge on both oxygen and flow characteristics is needed. In this way, it avoids any trial and error concerning the role of these environmental signals.

Footnotes

Acknowledgments

This work was funded by the K.U. Leuven IDO project 05/009 (QuEST) and is part of Prometheus, the Leuven Research and Development Division of Skeletal Tissue Engineering of the Katholieke Universiteit Leuven (www.kuleuven.be/prometheus).

Disclosure Statement

No competing financial interests exist.

Part of this work has been presented at the TERMIS-EU conference in June 2010 and at the 17th congress of the European Society of Biomechanics in July 2010.

References

Porter

J.R.

, Ruckh

T.T.

, Popat

K.C.

Bone tissue engineering: a review in bone biomimetics and drug delivery strategies. Biotechnol Prog, 25:1539. 2009.

Grayson

W.L.

, Frohlich

, Yeager

, Bhumiratana

, Chan

M.E.

, Cannizzaro

, Wan

L.Q.

, Liu

X.S.

, Guoc

X.E.

, Vunjak-Novakovic

Engineering anatomically shaped human bone grafts. Proc Natl Acad Sci U S A, 107:3299. 2010.

Hutmacher

D.W.

, Schantz

J.T.

, Lam

C.X.F.

, Tan

K.C.

, Lim

T.C.

State of the art and future directions of scaffold-based bone engineering from a biomaterials perspective. J Tissue Eng Regen Med, 1:245. 2007.

Ikada

Challenges in tissue engineering. J R Soc Interface, 3:589. 2006.

Chen

H.-C.

, Hu

Y.C.

Bioreactors for tissue engineering. Biotechnol Lett, 28:1415. 2006.

Martin

, Vermette

Bioreactors for tissue mass culture: design, characterization, and recent advances. Biomaterials, 26:7481. 2005.

Butler

D.L.

, Goldstein

S.A.

, Guilak

Functional tissue engineering: the role of biomechanics. J Biomech Eng, 122:570. 2000.

Radisic

, Malda

, Epping

, Geng

, Langer

, Vunjak-Novakovic

Oxygen gradients correlate with cell density and cell viability in engineered cardiac tissue. Biotechnol Bioeng, 93:332. 2006.

Lewis

M.C.

, Macarthur

B.D.

, Malda

, Pettet

, Please

C.P.

Heterogeneous proliferation within engineered cartilaginous tissue: the role of oxygen tension. Biotechnol Bioeng, 91:607. 2005.

10.

Malda

, Woodfield

T.B.

, Van der Vloodt

, Kooy

F.K.

, Martens

D.E.

, Tramper

, Van Blitterswijk

T.B.F.

, Riesle

The effect of PEGT/PBT scaffold architecture on oxygen gradients in tissue engineered cartilaginous constructs. Biomaterials, 25:5773. 2004.

11.

Wendt

, Stroebel

, Jakob

, John

G.T.

, Martin

Uniform tissues engineered by seeding and culturing cells in 3D scaffolds under perfusion at defined oxygen tensions. Biorheology, 43:481. 2006.

12.

Naughton

G.K.

From lab bench to market critical issues in tissue engineering. Ann N Y Acad Sci, 961:372. 2002.

13.

Wendt

, Riboldi

S.A.

, Cioffi

, Martin

Bioreactors in tissue engineering: scientific challenges and clinical perspectives. Adv Biochem Eng Biotechnol, 112:1. 2009.

14.

Martin

, Wendt

, Heberer

The role of bioreactors in tissue engineering. Trends Biotechnol, 22:80. 2004.

15.

Lovett

, Lee

, Edwards

, Kaplan

D.L.

Vascularization strategies for tissue engineering. Tissue Eng Part B Rev, 15:353. 2009.

16.

Malda

, van den Brink

, Meeuwse

, Grojec

D.E.

, Martens , Tramper

et al. Effect of oxygen tension on adult articular chondrocytes in micro-carrier bioreactor culture. Tissue Eng, 10:987. 2004.

17.

Das

, Jahr

, van Osch

G.J.

, Farrell

The role of hypoxia in MSCs: considerations for regenerative medicine approaches. Tissue Eng Part B Rev, 16:159. 2010.

18.

Grayson

W.L.

, Zhao

, Bunnell

, Ma

Hypoxia enhances proliferation and tissue formation of human mesenchymal stem cells. Biochem Biophys Res Commun, 358:948. 2007.

19.

Lennon

D.P.

, Edmison

J.M.

, Caplan

A.I.

Cultivation of rat marrow-derived mesenchymal stem cells in reduced oxygen tension: effects on in vitro and in vivo osteochondrogenesis. J Cell Physiol, 187:345. 2001.

20.

Henderson

J.H.

, Ginley

N.M.

, Caplan

A.I.

, Niyibizi

, Dennis

J.E.

Low oxygen tension during incubation periods of chondrocyte expansion is sufficient to enhance postexpansion chondrogenesis. Tissue Eng Part A, 16:1585. 2010.

21.

et al. In vitro expansion of adipose-derived adult stromal cells in hypoxia enhances early chondrogenesis. Tissue Eng, 13:2981. 2007.

22.

, Grayson

W.L.

, Frohlich

, Vunjak-Novakovic

Hypoxia and stem cell-based engineering of mesenchymal tissues. Biotechnol Prog, 25:32. 2009.

23.

Volkmer

, Kallukalam

B.C.

, Maertz

, Otto

, Drosse

, Polzer

, Bocker

, Stengele

, Denitsa

, Mutschler

, Schieker

Hypoxic preconditioning of human mesenchymal stem cells overcomes hypoxia-induced inhibition of osteogenic differentiation. Tissue Eng Part A, 16:153. 2010.

24.

Malladi

, Xu

, Chiou

, Giaccia

A.J.

, Longaker

M.T.

Effect of reduced oxygen tension on chondrogenesis and osteogenesis in adipose-derived mesenchymal cells. Am J Physiol Cell Physiol, 290:C1139. 2006.

25.

Yeatts

A.B.

, Fisher

J.P.

Bone tissue engineering bioreactors: dynamic culture and the influence of shear stress. Bone, 2010In press doi:10.1016/j.bone.2010.09.138.

26.

Grayson

W.L.

, Bhumiratana

, Cannizzaro

, Chao

P.H.G.

, Lennon

D.P.

, Caplan

A.I.

, Vunjak-Novakovic

Effects of initial seeding density and fluid perfusion rate on formation of tissue engineered bone. Tissue Eng, 14:1809. 2008.

27.

Bancroft

G.N.

, Sikavitsast

V.I.

, van den Dolder

, Sheffield

T.L.

, Ambrose

C.G.

, Jansen

J.A.

, Mikos

A.G.

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

28.

Cartmell

S.H.

, Porter

B.D.

, Garcia

A.J.

, Guldberg

R.E.

Effects of medium perfusion rate on cell-seeded three-dimensional bone constructs in vitro. Tissue Eng, 9:1197. 2003.

29.

Brown

, Meenan

B.J.

Investigating the effects of fluid shear forces on cellular responses to profiled surfaces in-vitro: a computational and experimental investigation. Conf Proc IEEE Eng Med Biol Soc, 2007:5387. 2007.

30.

Cioffi

, Boschetti

, Raimondi

M.T.

, Dubini

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

31.

Botchwey

E.A.

, Pollack

S.R.

, El-Amin

, Levine

E.M.

, Tuan

R.S.

, Laurencin

C.T.

Human osteoblast-like cells in three-dimensional culture with fluid flow. Biorheology, 40:299. 2003.

32.

De Boodt

, Truscello

, Özcan

S.E.

, Leroy

, Van Oosterwyck

, Berckmans

, Schrooten

Bi-modular flow characterisation in tissue engineering scaffolds using computational fluid dynamics and particle imaging velocimetry. Tissue Eng Part C, 2010Epub ahead of print10.1089/ten.tec.2010.0107.

33.

Sengers

B.G.

, van Donkelaar

C.C.

, Oomens

C.W.J.

, Baaijens

F.P.T.

Computational study of culture conditions and nutrient supply in cartilage tissue engineering. Biotechnol Prog, 21:1252. 2005.

34.

Galbusera

, Cioffi

, Raimondi

M.T.

, Pietrabissa

Computational modeling of combined cell population dynamics and oxygen transport in engineered tissue subject to interstitial perfusion. Comput Methods Biomech Biomed Eng, 10:279. 2007.

35.

Cioffi

, Kuffer

, Strobel

, Dubini

, Martin

, 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.

36.

Maes

, Van Ransbeeck

, Van Oosterwyck

, Verdonck

Modeling fluid flow through irregular scaffolds for perfusion bioreactors. Biotechnol Bioeng, 103:621. 2009.

37.

Cochran

D.M.

, Fukumura

, Ancukiewicz

, Carmeliet

, Jain

R.K.

Evolution of oxygen and glucose concentration profiles in a tissue-mimetic culture system of embryonic stem cells. Ann Biomed Eng, 34:1247. 2006.

38.

Meinel

et al. Bone tissue engineering using human mesenchymal stem cells: effects of scaffold material and medium flow. Ann Biomed Eng, 32:112. 2004.

39.

Davisson

, Sah

R.L.

, Ratcliffe

Perfusion increases cell content and matrix synthesis in chondrocyte three-dimensional cultures. Tissue Eng, 8:807. 2002.

40.

Wendt

, Marsano

, Jakob

, Heberer

, Martin

Oscillating perfusion of cell suspensions through three-dimensional scaffolds enhances cell seeding efficiency and uniformity. Biotechnol Bioeng, 84:205. 2003.

41.

Saha

, Kirkham

, Wood

, Curran

, Yang

Comparative study of the chondrogenic potential of human bone marrow stromal cells, neonatal chondrocytes and adult chondrocytes. Biochem Biophys Res Commun, 401:333. 2010.

42.

Stojkovic

, Lako

, Strachan

, Murdoch

Derivation, growth and applications of human embryonic stem cells. Reproduction, 128:259. 2004.

43.

Simpson

, McGrath

, Savion

Myocyte hypertrophy in neonatal rat heart cultures and its regulation by serum and by catecholamines. Circ Res, 51:787. 1982.

44.

Demol

, Lambrechts

, Geris

, Schrooten

, Van Oosterwyck

Towards a quantitative understanding of oxygen tension and cell density evolution in fibrin hydrogels. Biomaterials, 2010, DOI:10.1016/j.biomaterials.2010.08.093.

45.

Byrne

D.P.

, Lacroix

, Plannel

J.A.

, Kelly

D.J.

, Prendergast

P.J.

Simulation of tissue differentiation in a scaffold as a function of porosity, Young's modulus and dissolution rate: application of mechanobiological models in tissue engineering. Biomaterials, 28:5544. 2007.

46.

Milan

J.L.

, Planell

, Lacroix

Simulation of bone tissue formation within a porous scaffold under dynamic compression. Biomech Model Mechanobiol, 9:583. 2010.

47.

Kelly

D.J.

, Prendergast

P.J.

Prediction of the optimal mechanical properties for a scaffold used in osteochondral defect repair. Tissue Eng, 12:2509. 2006.