Structural optimization design of perfusion bioreactors with multilayer circular parallel plates scaffold

Abstract

Perfusion bioreactors with multilayer circular parallel plates scaffold have the potential to provide a considerable scalability in cells’ cultivation for tissue engineering. By now, relatively little attention has been focused on the uniform distribution of fluid shear stress (FSS) among the parallel plates. However, the magnitude of FSS has a significant effect on both proliferation and differentiation of cells. This paper analyses the relationship between the structural parameters of the scaffold and the distribution of FSS among the plates and then presents an approach to optimize the flow field within the bioreactors. By regulating the structural parameters of the scaffold, the uniform distribution of FSS on the 23-layers plates can be achieved. The resolution of FSS can be up to 0.001 Pa and the range of error restrained within ±0.00015 Pa. By means of structural optimization, the discrepancy of the flow field among the plates can be effectively diminished.

Keywords:

• perfusion bioreactor
• multilayer circular parallel plates scaffold
• fluid shear stress
• computational fluid dynamics

Introduction

Perfusion bioreactors with multilayer parallel plates scaffold have a potential in cell cultivation because of their good scalability in tissue engineering bioreactors.[1,2] Numerous researchers have implemented diverse studies based on bioreactors with this structure. Chen et al. [3] adopted the multilayer parallel plate structure to design bioreactors. Then Park et al. [4] and De Bartolo et al. [5] designed a flat membrane bioreactor stacked with 20-layer plates. Afterward, Zhang et al. [6] presented a bioreactor with radial flow within 50 layers of circular plates. Similarly, other authors [7–15] utilized a cylindrical bioreactor composed of stacked double sandwich culture plates. Among these reports, little information is available about the distribution discrepancy of fluid shear stress (FSS) on each plate in different depth within the medium. However, the magnitude of FSS has a significant influence on cell proliferation and cell differentiation.[16,17] Thus, the discrepancy of FSS distribution on the different culture plates within one bioreactor will cause inconsistent cell cultivation, which is undesired in tissue engineering.

This paper is focused on implementing a uniform distribution of FSS on each plate in different depth under the culture medium's surface by optimizing the structural parameters of the bioreactor. In terms of relevant contents from many researches, such as Cherry and Papoutsakis,[18] Bancroft et al.,[19] Godwin et al. [20] and Porter et al.,[21] the level of FSS in their studies was up to 0.001 Pa. In this situation, the distribution uniformity of FSS on every plate had to be analysed. It is necessary to seek for an approach to reduce the discrepancy of the FSS distribution. This paper studies the relation between the FSS distribution and the structural parameters of the scaffold using the simulation technology of computational fluid dynamics (CFD).[22] Based on this relation, a method for optimising the flow field on every layer with 0.001 Pa FSS was presented. As a result, some cells, which are very sensitive to FSS, such as bone marrow mesenchymal stem cells, will obtain more consistent proliferation rate in the perfusion bioreactor, optimized by this approach.

Materials and methods

Cylindrical perfusion bioreactor

A cylindrical perfusion bioreactor was constructed from a circular glass tank, which was used as a perfusion room that contained 23 layers of plates, stacked concentrically through a vertical stainless steel pipe, which served as a draft tube. In each interstitial space between every two adjacent plates, there were six deflector holes pierced through the wall of the steel pipe. A rectangular glass vessel, which served as a buffer chamber adhered to the circular tank, segregated by a baffle. Cultural medium was first pumped into the buffer chamber through a silicon hose and then it was spilled over into the perfusion tank through the window on the top of the baffle. Afterward, the fluid flowed into the draft steel pipe. Finally, driven by a peristaltic pump, the medium flowed into the outlet pipe coupled with the end of the draft steel pipe. The prototype of the configuration is shown in Figure 1.

Figure 1. Prototype of the bioreactor with multilayer plates scaffold.

Note: Perfusion tank (1), support of perfusion tank (2), silicon hose (3), buffer chamber (4), baffle (5), multilayer plate scaffold (6), circular plate (7), draft tube (8), deflector hole (9), outlet of draft tube (10), outlet pipe (11), support of plate (12).

The scale of this type of bioreactor can be easily changed by altering the number of layers.

Simulation in computational fluid dynamics

The geometry of the bioreactor was created by using the commercial software (Solidworks^TM 2012). As shown in Figure 2(a), D1 (100 mm) and H1 (2 mm) represent the diameter and thickness of the circular plates, respectively. W (3 mm) denotes the height of the interstitial space between the neighbouring plates. D2 (11 mm) and D3 (18 mm) are the inner and the outer diameter of the central draft pipe, respectively. The fluid height of 140 mm (H2) and the fluid diameter of 126 mm (D4) were used to simulate the flow. The upper surface of the first plate was located at 17.6 mm under the fluid's surface. These structural parameters were kept constant in all simulations. The meshing and fluid calculations were carried out by the CFD package (ANSYS™ 14.5). As depicted in Figure 2(b), for the symmetrical geometry, simulated with CFD, only quarter of the computational domain along the vertical axis needed consideration. Non-structural tetrahedral meshing, isothermal, incompressible and Newtonian fluid with a dynamic viscosity of 0.0025 kg/(m s) and a density of 1007.4 kg/m³ were adopted in the simulations. At the top the fluid had an inlet velocity set at 0.000381 m/s (i.e. a flow rate of 4.75 mL/s) and the outlet of the central draft pipe was set as outflow. A laminar model was utilized ignoring the energy transfer process. These computational parameters also remained invariable in all simulations.

Figure 2. Cut-open view of multilayer plates configuration and computational grid.

Note: Cut-open view (a), computational grid for simulation (b). Diameter of the circular plates (D1), thickness of the circular plates (H1), height of the interstitial space between the neighbouring plates (W); inner (D2) and the outer (D3) diameter of the central draft pipe; fluid height (H2) and fluid diameter (D4).

Six groups of simulations with a deflector hole aperture of 1, 1.2, 1.4, 1.6, 1.8 and 2 mm were implemented for comparison. In each group the sizes of all deflector holes in the draft tube were the same.

The unit of FSS is Pa, calculated by CFD software FLUENT®, if not stated otherwise.

Statistical analysis

The simulated data of wall shear stress within a radius of 15–45 mm of the circular region on the top of the surface of each plate were exported into the software MATLAB. The software was set to get the mean values of FSS on every plate. The data of FSS on the first plate under the surface of the liquid was precluded because of the complex flow field upon it and because there were no cells cultured on it. Thus, to simplify the description, hereafter ‘the first plate’ represented the plate that was located at 22.6 mm depth under the surface of the liquid.

Six curves composed of six groups of mean values of FSS on 23-layers plates were separately fitted with exponential regression analysis category in software Origin8. Nonlinear curve fitting function, ExpDec1, was adopted. The formulation of ExpDec1 function was (1) $y=y_0+A_1\cdot \exp \left(\frac{x}{t_1}\right),x=22.6+i\times 5,(i=0,1,2,\dots ,22)$(1) where $y_0$, $A_1$ and $t_1$ are variables of regression coefficients and x, y indicate the depth (mm) under the surface of the liquid and FSS (Pa), respectively.

Finally, in terms of the distribution of these FSS regression coefficients, which changed with the deflector hole size, suitable regression modes were employed to fit the curves about coefficients $y_0$ , $A_1$ and $t_1$ . P < 0.05 was considered significant.

Results and discussion

The simulation result of FSS distribution on 23-layer plates is shown in Figure 3(a). Figure 3(b) depicts the correspondent path-lines of the simulated velocity flow.

Figure 3. Distribution of fluid shear stress on every plate and flow velocity in vertical section.

Note: Distribution of FSS, the maximal value of y-coordinate is 5.00e−03 Pa (a); path-lines for flow velocity, the maximal of y-coordinate value is 9.00e−02 Pa (b).

Six groups of mean values ±SD of FSS on 23-layer plates are depicted in Figure 4, in which each colour curve indicates a kind of a deflector hole size.

Figure 4. Distribution of FSS on 23-layer plates, with six kinds of deflector hole sizes.

The results of the regression analysis for the six curves constructed from those mean values mentioned above are shown in Figure 5, in which the x-coordinate represents the deflector hole size and the y-coordinate represents the magnitude of the coefficients from Equation (1)(1) $y=y_0+A_1\cdot \exp \left(\frac{x}{t_1}\right),x=22.6+i\times 5,(i=0,1,2,\dots ,22)$(1) . According to the distribution characteristic of curves in Figure 5, linear, two-order polynomial and exponential categories were performed to fit the coefficient curves y₀, $A_1$ and t₁, respectively. The correspondent fitting functions are represented as follows:

Figure 5. The distribution of FSS regression coefficients y₀, A₁, t_1.changing with deflector hole size.

Linear fitting function for $y_0$: (2.1) $y_0=a+b\cdot h$(2.1)

Two-order polynomial fitting function for $A_1$: (3.1) $A_1=c+B_1\cdot h+B_2\cdot h^2$(3.1)

Exponential fitting function for $t_1$: (4.1) $t_1=t_0+A_2\cdot \exp \left(\frac{h}{t_2}\right)$(4.1) where h is the deflector hole aperture (mm) and $a$, $b$, c, $B_1$, $B_2$, $t_0$, $A_2$ and $t_2$ are regression coefficients. The fitting results are shown in Figure 6 and the formulations are (2.2) $y_0=0.00179-8.51808\times {10}^{-4}\cdot h$(2.2) (3.2) $A_1=-4.93405\times {10}^{-4}+9.04065\times {10}^{-4}\cdot h-3.28002\times {10}^{-4}\cdot h^2$(3.2) (4.2) $t_1=63.96826-1.33435\cdot \exp \left(\frac{h}{0.5666}\right)$(4.2)

Figure 6. Fitting results of curves from Figure 5.

By joining Equations (2.2)(2.2) $y_0=0.00179-8.51808\times {10}^{-4}\cdot h$(2.2) , (3.2(3.2) $A_1=-4.93405\times {10}^{-4}+9.04065\times {10}^{-4}\cdot h-3.28002\times {10}^{-4}\cdot h^2$(3.2) ), (4.2(4.2) $t_1=63.96826-1.33435\cdot \exp \left(\frac{h}{0.5666}\right)$(4.2) ) and (1(1) $y=y_0+A_1\cdot \exp \left(\frac{x}{t_1}\right),x=22.6+i\times 5,(i=0,1,2,\dots ,22)$(1) ), the function about x, h and y can be represented as follows: (5) $\begin{array}{l}y\left(x,h\right)=y_0+A_1\cdot \exp \left(\frac{x}{t_1}\right)\\ \{\begin{array}{l}{y}_0=0.00179-8.51808\times {10}^{-4}\cdot h\\ A_1=-4.93405\times {10}^{-4}+9.04065\times {10}^{-4}\cdot h-3.28002\times {10}^{-4}\cdot h^2\\ t_1=63.96826-1.33435\cdot \exp \left(\frac{h}{0.5666}\right)\end{array}\end{array}x=22.6+i\times 5,(i=0,1,2,\dots ,22),h\in [0.5,2.5]$(5) where y is FSS, and x and h separately represent depth under the surface of the liquid and the deflector hole aperture. The range of h was defined by the parameter W as shown in Figure 2(a).

The CFD method provides effective insights into the optimization and design of scaffolds and a full characterization of flow fields within the bioreactors. This technology allows us to predict the flow without carrying out numerous bioreactor experiments.

Based on the CFD simulation results shown in Figure 3(b), the flow field among each parallel plate exhibited a steady laminar flow. Thus, the FSS on every plate can be calculated as given in the following [23]: (6) $\tau =\frac{6\mu Q}{d{w}^2}$(6) where τ is FSS (Pa), μ is the dynamic viscosity (kg/(m s)), Q is the flow rate (mL/s), d is the outer arc length (cm) of each plate and w is the vertical height (cm) of the fluid among every plate.

In light of Equation (6)(6) $\tau =\frac{6\mu Q}{d{w}^2}$(6) , FSS linearly follows the change of Q, while μ, d and w keep constant. Consequently, considering Equation (5)(5) $\begin{array}{l}y\left(x,h\right)=y_0+A_1\cdot \exp \left(\frac{x}{t_1}\right)\\ \{\begin{array}{l}{y}_0=0.00179-8.51808\times {10}^{-4}\cdot h\\ A_1=-4.93405\times {10}^{-4}+9.04065\times {10}^{-4}\cdot h-3.28002\times {10}^{-4}\cdot h^2\\ t_1=63.96826-1.33435\cdot \exp \left(\frac{h}{0.5666}\right)\end{array}\end{array}x=22.6+i\times 5,(i=0,1,2,\dots ,22),h\in [0.5,2.5]$(5) established at Q = 4.75 mL/s and μ, d, w remained invariable, Equation (5)(5) $\begin{array}{l}y\left(x,h\right)=y_0+A_1\cdot \exp \left(\frac{x}{t_1}\right)\\ \{\begin{array}{l}{y}_0=0.00179-8.51808\times {10}^{-4}\cdot h\\ A_1=-4.93405\times {10}^{-4}+9.04065\times {10}^{-4}\cdot h-3.28002\times {10}^{-4}\cdot h^2\\ t_1=63.96826-1.33435\cdot \exp \left(\frac{h}{0.5666}\right)\end{array}\end{array}x=22.6+i\times 5,(i=0,1,2,\dots ,22),h\in [0.5,2.5]$(5) can be adjusted as (7) $\begin{array}{l}f\left(x,h\right)=\frac{Q}{4.75}\times y\left(x,h\right)\\ y\left(x,h\right)=y_0+A_1\times \exp \left(\frac{x}{t_1}\right)\\ \{\begin{array}{l}{y}_0=0.00179-8.51808\times {10}^{-4}\cdot h\\ A_1=-4.93405\times {10}^{-4}+9.04065\times {10}^{-4}\cdot h-3.28002\times {10}^{-4}\cdot h^2\\ t_1=63.96826-1.33435\times \exp \left(\frac{h}{0.5666}\right)\end{array}\begin{array}{l},x=22.6+i\times 5,(i=0,1,2,\dots ,22)\hfill \\ ,h\in [0.5,2.5]\hfill \end{array}\end{array}$(7) where f is FSS based on the variable Q.

Furthermore, when only d and w remain constant, the formulation about FSS and the variable Q, μ can be similarly established as follows: (8) $\begin{array}{l}\text{fss}\left(x,h\right)=\frac{\mu }{0.0025}\times f\left(x,h\right)\\ f\left(x,h\right)=\frac{Q}{4.75}\times y\left(x,h\right)\\ y\left(x,h\right)=y_0+A_1\times \exp \left(\frac{x}{t_1}\right)\\ \{\begin{array}{l}{y}_0=0.00179-8.51808\times {10}^{-4}\cdot h\\ A_1=-4.93405\times {10}^{-4}+9.04065\times {10}^{-4}\cdot h-3.28002\times {10}^{-4}\cdot h^2\\ t_1=63.96826-1.33435\times \exp \left(\frac{h}{0.5666}\right)\end{array}\begin{array}{l},x=22.6+i\times 5,(i=0,1,2,\dots ,22)\hfill \\ ,h\in [0.5,2.5]\hfill \end{array}\end{array}$(8) where fss is FSS based on variable Q and μ.

According to the transverse changing tendency of FSS curves in Figure 4, the FSS increased with the increase of the layer depth. Comparably, according to the vertical changing trend of FSS values in Figure 4, the varying direction of FSS magnitude was divided at approximately 15th layer (i.e. before the 15th layer, FSS values along the y-coordinate decreased with the increase of deflector hole aperture; after the 15th layer, FSS values increased with the increase of deflector hole aperture). Based on these FSS distribution rules, an optimization method to uniform the FSS distribution among each plate by regulating the deflector hole aperture can be obtained intuitively: the aperture of the deflector holes had to decrease appropriately along both changing directions of the layer number from the 15th layer. The specific adjusting approach for the deflector hole size is depicted in Figure 7.

Figure 7. Diagram of adjusting the deflector hole size.

In Figure 7, because the beginning layer, i.e. the 15th layer (x = 92.6 mm), was obtained intuitively from Figure 4, the hole apertures had to be changed by both increased and decreased steps in the process of adjusting the deflector hole size, in order to improve the accuracy of the outcome. To exemplify the effect of the optimizing method, a smaller value h = 1.2 mm and a corresponding FSS = 0.00145 Pa were selected to perform the adjusting process shown in Figure 7. The adjusted results of the deflector hole size and the corresponding FSS distributions, deduced by Equation (5)(5) $\begin{array}{l}y\left(x,h\right)=y_0+A_1\cdot \exp \left(\frac{x}{t_1}\right)\\ \{\begin{array}{l}{y}_0=0.00179-8.51808\times {10}^{-4}\cdot h\\ A_1=-4.93405\times {10}^{-4}+9.04065\times {10}^{-4}\cdot h-3.28002\times {10}^{-4}\cdot h^2\\ t_1=63.96826-1.33435\cdot \exp \left(\frac{h}{0.5666}\right)\end{array}\end{array}x=22.6+i\times 5,(i=0,1,2,\dots ,22),h\in [0.5,2.5]$(5) , are separately shown in Figures 8(b) and 9(b). For comparison, the original distributions of the hole size and FSS from Figure 4 are separately exhibited in Figures 8(a) and 9(a). To simplify the adjusting process, the optimization strategy was established on the preconditions that the adjusting step n and the variation of h remained small enough so that the cross-correlation among every deflector hole size can be ignored. According to the range of h, n = 0.1 mm, i = +0.5 mm and j = −0.5 mm were adopted, where i and j were the searching range of n (i.e. ∑ n ∈ [j, i]). Based on these preconditions the practical FSS distribution under the optimized structural parameters will have discrepancy to the result obtained via Equation (5)(5) $\begin{array}{l}y\left(x,h\right)=y_0+A_1\cdot \exp \left(\frac{x}{t_1}\right)\\ \{\begin{array}{l}{y}_0=0.00179-8.51808\times {10}^{-4}\cdot h\\ A_1=-4.93405\times {10}^{-4}+9.04065\times {10}^{-4}\cdot h-3.28002\times {10}^{-4}\cdot h^2\\ t_1=63.96826-1.33435\cdot \exp \left(\frac{h}{0.5666}\right)\end{array}\end{array}x=22.6+i\times 5,(i=0,1,2,\dots ,22),h\in [0.5,2.5]$(5) . Thus, the optimized structure was validated in CFD tool. The result of the validated FSS is shown in Figure 9(c).

Figure 8. Distribution of deflector hole size (h) at different depth (x).

Note: Original distribution of deflector hole size (a); optimized distribution of deflector hole size (b).

Figure 9. Distribution of FSS at different depth (x).

Note: Original distribution of FSS (a); optimized distribution of FSS (b); verified distribution of FSS (c). Fluid shear stress (FSS); computational fluid dynamics (CFD).

The validated result (the curve shown in Figure 9(c)) indicated that the uniformity of practical FSS distribution under the optimized structural parameters was the best among the three curves in Figure 9, which can be demonstrated in Table 1.

Table 1. Mean value, variance and peak-to-peak value of curves in Figure 9.

Curve	Mean value (Pa)	Variance (Pa)	Peak-to-peak value (Pa)
A	0.00137	14.9288E−8	0.0012
B	0.00127	4.20158E−8	0.0006
C	0.001278	2.72332E−8	0.0006

Note: a, b and c indicate the curves a, b and c as shown in Figure 9.

As shown in Table 1, curve ‘c’ had the smallest variance and peak-to-peak value, which denoted that its uniformity of distribution was the best.

The verified results indicate that the optimizing approach presented in this paper is efficient to improve the uniformity of FSS distribution among parallel plates. By the optimizing strategy the consistent FSS with 0.001 Pa level can be provided for a large-scale cells cultivating, sensitive to FSS, such as bone marrow mesenchymal stem cells.

Conclusions

The uniformity of the flow field among multilayers parallel plates affects the scalability of perfusion bioreactors for the cultivation of cells. Within bioreactors with multilayers parallel plates scaffold, FSS discrepancy of distribution exists among plates even when the scaffold has perfectly the same structural parameters. The extent of this discrepancy intensifies as the number of layers increases, which diminishes the scalability of bioreactors. The analysis of the flow of the fluid within this structure, using CFD technology, indicated that the FSS distribution and the structural parameters of the deflector hole aperture among the parallel plates exhibited certain relationship. Based on this relationship, the uniformity of FSS distribution among the parallel plates could be improved by adjusting the structural parameters of the size of the deflector holes and the optimization design of the perfusion bioreactors with multilayers parallel plates scaffold could be realized.

Disclosure statement

No conflict of interest exits in this paper, and the paper is approved by all authors for publication.

Additional information

Funding

This work was supported by the National Natural Science Foundation of China [grant number 51205421]; Postdoctoral Science Foundation of China [grant number 2012M521647]; Science and Technology Program of Guangzhou, China [grant number (2013)164; 2060204]; Key Laboratory of Sensing Technology and Biomedical Instruments of Guangdong province, China [grant number 2011A060901013].