Distributed and Lumped Parameter Models for the Characterization of High Throughput Bioreactors

Laura Iannetti; Giovanna D’Urso; Gioacchino Conoscenti; Elena Cutrì; Rocky S. Tuan; Manuela T. Raimondi; Riccardo Gottardi; Paolo Zunino

doi:10.1371/journal.pone.0162774

Abstract

Next generation bioreactors are being developed to generate multiple human cell-based tissue analogs within the same fluidic system, to better recapitulate the complexity and interconnection of human physiology [1, 2]. The effective development of these devices requires a solid understanding of their interconnected fluidics, to predict the transport of nutrients and waste through the constructs and improve the design accordingly. In this work, we focus on a specific model of bioreactor, with multiple input/outputs, aimed at generating osteochondral constructs, i.e., a biphasic construct in which one side is cartilaginous in nature, while the other is osseous. We next develop a general computational approach to model the microfluidics of a multi-chamber, interconnected system that may be applied to human-on-chip devices. This objective requires overcoming several challenges at the level of computational modeling. The main one consists of addressing the multi-physics nature of the problem that combines free flow in channels with hindered flow in porous media. Fluid dynamics is also coupled with advection-diffusion-reaction equations that model the transport of biomolecules throughout the system and their interaction with living tissues and C constructs. Ultimately, we aim at providing a predictive approach useful for the general organ-on-chip community. To this end, we have developed a lumped parameter approach that allows us to analyze the behavior of multi-unit bioreactor systems with modest computational effort, provided that the behavior of a single unit can be fully characterized.

Citation: Iannetti L, D’Urso G, Conoscenti G, Cutrì E, Tuan RS, Raimondi MT, et al. (2016) Distributed and Lumped Parameter Models for the Characterization of High Throughput Bioreactors. PLoS ONE 11(9): e0162774. https://doi.org/10.1371/journal.pone.0162774

Editor: Feng Zhao, Michigan Technological University, UNITED STATES

Received: May 9, 2016; Accepted: August 29, 2016; Published: September 26, 2016

Copyright: © 2016 Iannetti et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Data Availability: The manuscript is based on computational models and numerical simulations that are fully reproducible to anyone who has access to the commercial computational platform used for the calculations and named in full details in the manuscript. The authors confirm that there isn't any additional data that is not already provided within the paper.

Funding: This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement No. 646990 - NICHOID) awarded to MTR. These results reflect only the authors’ view and the Agency is not responsible for any use that may be made of the information contained. The work of RG and RST was supported by grants from the Commonwealth of Pennsylvania, Department of Health, the Environmental Protection Agency (83573601), the National Institutes of Health (1U18 TR000532), the U.S. Department of Defense (W81XWH-08-2-0032 and W81XWH-10-1-0850), and the Ri.MED Foundation.

Competing interests: The authors have declared that no competing interests exist.

1 Introduction

A number of in vitro approaches have been used over time for high throughput drug screening or toxicology testing. However, most currently available systems are only partial approximations of human biology and their predictive capacity is consequently limited. In fact, such systems are either based on human cell cultures, not capturing the complexity of cell behavior in a three dimensional (3D) environment, or they are based on animal tissues fragments, 3D in nature but only partially biosimilar to human tissues and unable to account for interactions with other organs. To overcome these limitations, next generation bioreactors are being developed to generate multiple human cell-based tissue analogs within the same fluidic system to better recapitulate the complexity and interconnection of human physiology. These efforts aim at creating multi-tissue organ systems (cardiovascular, gastro-intestinal, musculoskeletal, etc.) that ultimately can be joined in an interconnected human-on-chip device capable of providing a veritable representation of the body complex response to diseases and potential drug treatments [3–5].

The effective development of these devices requires a solid understanding of their interconnected fluidics, to predict the transport of nutrients and waste through the constructs and improve the design accordingly. In this work, we have focused on a specific bioreactor with multiple input/output aimed at generating osteochondral constructs, i.e., a biphasic constructs in which one side is cartilaginous in nature, while the other is osseous. This bioreactor [1, 6, 7] represented in Fig 1 has been chosen since it comprises both a dual chamber system to host a single biphasic tissue construct with distinct fluidics (Fig 1, top), and a set of interconnected chambers with common fluidics (Fig 1, bottom). Starting from this specific bioreactor, we have developed a general approach to model the microfluidics of a multi-chamber, interconnected system that may be applied to human-on-chip devices.

Download:

Fig 1. Different bioreactor configurations.

1 cell (top left), 1-unit in cross section (top right), 4-units (bottom left) and 96-units. (bottom right).

https://doi.org/10.1371/journal.pone.0162774.g001

The microphysiological osteochondral bioreactor analyzed in this work is aimed at the study of osteoarthritis (OA), a major pathology of articular joints, affecting over 33% of the population over the age of 65 [8]. The hallmark of this disease that affects all tissues in the joint, is the progressive degeneration of cartilage which begins well before clinical symptoms manifest, ultimately requiring joint replacement surgery [9]. The high incidence of this painful and disabling pathology begs for the understanding of the causes and mechanisms of its development, in order to identify reparative drug therapies to arrest or even regenerate the damaged tissues and ultimately avoid surgery. A novel strategy in this respect adopts a tissue engineering approach and the use of bioreactors [1, 7] to generate a high number of identical in vitro constructs that can replicate the pathogenesis of joint diseases for the identification of therapeutic targets and for drug screening [1, 10–12]. Critical in this respect is the development of a representative model of the interactions between cartilage and other joint tissues and, in particular, with the subchondral bone. In fact, there is growing evidence of the exchange of nutrients, cytokines, and hormones in vivo between bone and cartilage. The osteochondral (OC) unit is then conceived as the main target of OA, to reflect the dynamic cartilage/bone interplay in both health and disease [1, 13–18]. The medium to high throughput system studied in this work, which we call high-throughput bioreactor (HTB) hereon, is the first of its kind. It hosts in a single chamber a biphasic construct, with separate fluidics for its cartilaginous and osseous components, effectively creating a dual-chamber setup (Fig 1) [6, 7]. In this way, cartilage and bone will be in contact and able to signal to each other, while each is exposed to its ideal culture medium. Furthermore, the HTB allows the generation and culture of a high number of identical OC constructs similar in dimensions to native tissue biopsies [1, 6, 7]. It must be noted that the physiological functions of the examined tissue are primarily load bearing and force transduction, which imply a key role for the extracellular matrix (ECM), also an essential player in the regulation of cell differentiation, physiology and response to insults [1, 19, 20]. Consequently, a bioreactor that accommodates a significant ECM tissue component to recapitulate at least some of the physiological aspects of the osteochondral complex requires a relatively larger volume, in the order of millimeters rather than the hundreds of micrometers more common in microfluidic systems. To generate a construct that mimics tissue physiology, the bioreactor chamber is filled with a cell-laden porous polymeric scaffold. Hence, the larger size and the presence of porous scaffold within the insert makes nutrient perfusion within the device a potential challenge, since to avoid cellular hypoxia and to obtain adequate tissue development, nutrients must travel a longer path to reach the inner regions within the bioreactor. In this context, we use computational fluid dynamics to assess the hydrodynamic properties of the system. Previous works [4, 21–23] evaluated the fluid mixing and transport of nutrients between chambers in the same unit of a forced perfusion setup, but to our knowledge there are no similar studies about the interaction of fluid and porous constructs in a design with more effective fluidics as the one in Fig 1.

Furthermore, to achieve a high-throughput drug screening system, single bioreactor dual-chambers (bioreactor unit) have been connected and combined in a multi-unit system, organized in sequential and parallel rows (Fig 1). In the 96 wells design presented in Fig 1, individual units are connected only in series, 8 at a time as this design is best suited for drug or toxicological screening; to asses for instance a dose response, each array of 8 units can be subjected to a different concentration of the compound under examination. “In parallel” connection, although possible, has not been envisioned. The constructs in each row are meant as replicates for multiple endpoint testing (e.g., histology, PCR, etc.). A further challenge is then to guarantee that the tissue constructs in the downstream chambers receive the appropriate amount of nutrients from the fluid that has perfused the units upstream. In other words, not only a dual-chamber bioreactor, but also a multi-unit array shall be analyzed.

The specific objective of this work is to develop a methodology to characterize the flow and transport in a HTB by means of a computational modeling approach, combining distributed and lumped parameter models. In particular, we have assessed the degree of perfusion and mixing of nutrients in each region of the device, evaluating the effect of different scaffold types. The computational model was then used to compare two different engineered constructs, a hydrogel (methacrylated gelatin, GelMA [6, 24]) and a porous polymeric scaffold (poly-L-lactate, PLLA)[25]. The first one features very small pore size and is solute permeable, the second one shows larger pore size and is impenetrable to fluid and nutrients.

Performing such simulations requires overcoming several challenges at the level of computational modeling. The main one consists of addressing the multi-physics nature of the problem that combines free flow in channels with hindered flow in porous media. Fluid dynamics is then coupled with advection-diffusion-reaction equations that model the transport of biomolecules throughout the system and their interaction with living tissue. Besides these modeling challenges, the complex configuration of the bioreactor poses significant difficulties in building the CAD model and discretizing its parts with a computational mesh suitable for the application of a numerical scheme. These issues can be solved using an in-house-made software that incorporates state-of-the-art efficient algorithms for the approximation of partial differential equations. Although this approach is viable, it entails significant costs in terms of man-hours for the implementation and validation of the new software. For this reason, we have adopted here a commercial platform, ANSYS (ANSYS Inc., Canonsburg, PA), which features advanced multi-physics simulation capabilities. Another challenging aspect of this work is then to stretch the limits of the ANSYS platform to address the complex problem at hand. Ultimately, our aim is to provide a predictive approach useful for the general organ-on-chip community. To this end, we have developed a lumped parameter approach that allows us to analyze the behavior of multi-unit bioreactor systems with a modest computational effort, provided that the behavior of a single unit could be fully characterized. If the linearity conditions are satisfied, this computational methodology is independent from the specific osteochondral nature of the biological system being studied. Our approach simply describes a network of interconnected multi-chamber units. Consequently, we believe that our approach can be directly applied to predict the flow and transport of a generic human-on-chip setup, even those comprising multiple physiological systems (e.g., a liver model connected to a kidney model, connected to a bone model, etc.) with single or multi-chamber units.

2 Models and Methods

Exploiting the commercial platform ANSYS (ANSYS Inc., Canonsburg, PA), we have developed a CAD model of the bioreactor and we have used it to simulate flow and transport phenomena in the system. The steps to achieve a realistic simulation of the fluid and transport within the bioreactor are detailed below.

2.1 CAD model

The 3D CAD model of the bioreactor was created using ANSYS ICEM CFD v.15.0 (ANSYS Inc.) CAD modeler. We have considered a row of 4-units connected in series (see Fig 1). Each unit has the same configuration, specifically designed to grow a construct that combines cartilage and bone, and comprises the following parts: two inlets and two outlets consisting of cylindrical channels, to guarantee the circulation of fluid from the upstream units to the downstream ones. Each inlet/outlet channel is characterized by a length (L) of 5.3 mm and an inner diameter (d) of 1 mm. The perforated cylindrical insert that holds the scaffold in place is 8.5 mm high and 3.75 mm wide. Each bioreactor chamber is sealed by an upper cap and by two O-rings (see Fig 1). Forthcoming extensions of this study will consider rows of 8 bioreactor units. By aligning 12 parallel lines of these rows, one obtains a plate of 96-units, which is a realistic prototype of high-throughput bioreactor for drug screening.

2.2 Flow

The bioreactor features the combination of free flow for the inlets, outlets, and the outer chambers with porous media flow for the inner culture chamber (insert). In each region, we assume that the flow is incompressible. For momentum balance, our approach employs a general equation that encompasses the nature of both types of flow, and we will switch between them by suitably tuning the problem parameters in each region. This equation has the structure of Brinkman equation for flow in porous media, because it combines viscous terms, such as in Stokes, with friction terms, such as in Darcy. To model free flow, a convective term, which plays a significant role in case of high Reynolds regimes, was added. Static conditions are also assumed. Then, the momentum balance equation reads as follows: (1) where denotes the velocity vector field ( and denote the restriction of the velocity field to the free fluid and porous medium, respectively), p the hydrostatic pressure, ρ e μ are the fluid viscosity and density respectively, and K_perm the hydraulic conductivity of the porous medium (for the free flow regions we set K_perm → ∞). For the partition of the bioreactor into sub-regions, we refer to Fig 2. We assume that the culture medium that perfuses the bioreactor is comparable to water (ρ = 999,97 kg/m³ μ = 0,001 Pa s) since the dissolved nutrients and other chemical species are relatively dilute.

Download:

Fig 2. Representation of the bioreactor.

Free fluid regions are visualized in grey, the porous medium is red (top). For the localization of boundary surfaces, Ω and Г indicate volume and surface, respectively (bottom).

https://doi.org/10.1371/journal.pone.0162774.g002

For the definition of boundary conditions, we partition the bioreactor surface as illustrated in Fig 2. At the bioreactor inlet, (Γ_{c_up,in} e Γ_{c_down,in}), a given flow rate is applied through the enforcement of a flat velocity profile on the inflow sections; a no-slip condition is adopted on the surfaces that separate the free fluid and the porous medium from the bioreactor walls (), which have been assumed to be rigid walls. At the outlet, (Γ_{c_up,out} e Γ_{c_down,out}) we have set a uniform normal stress field equal to the atmospheric pressure, namely σ_f ∙ n = 0, where is the Cauchy stress in the fluid. Given the previous modeling choices, the flow problem becomes (2)

2.3 Mass transport

An important part of this study consists of modeling the transport of bio-molecules dissolved in the culture media that perfuse the bioreactor. In particular, we focus on oxygen, fundamental to guarantee cell survival. However, the model is general and has been used to describe the transport of glucose and proteins, as it will be reported in forthcoming works.

Since all solutes are diluted, they are modeled as passively transported by the culture media. Their governing equations have been formulated in terms of volumetric concentrations measured in [mg/ml]. The symbol C denotes the solute concentration, D the diffusion coefficient for the specific biomolecule and the subscripts f and s indicate the fluid and the porous medium (scaffold), respectively. Therefore the equation describing the biomolecules’ transport in the fluid phase is: (3)

For the porous medium, namely the scaffold region, we assume that fluid and solid phases coexist. We denote with C_s,s and C_s,f the volumetric concentration of biomolecules in the solid and in the fluid phase of the scaffold, respectively. Denoting with γ the porosity of the scaffold (complement to unity of the solid phase, i.e. for the free flow regions we set γ = 1), the volumetric concentration of biomolecules in the porous medium is given by the following weighted average C_s = γC_s,f + (1 − γ)C_s,s. Then, following the theory of mixtures, the governing equations for biomolecules concentration in the porous medium read as follows: (4)

This model assumes that both the fluid and the solid phases in the porous medium are permeable to biomolecules. The mass transfer coefficient from the fluid to the solid phase in the porous medium is τ, while IAD is the interface area density of the surface separating the two phases. As a result, the term τIAD(C_s,s − C_s,f) represents the flux exchanged between the two phases of the porous medium. The symbol S denotes the source term representing the consumption of nutrients by living cells disseminated into the scaffold. For this reason, it is usually a function (linear or nonlinear) of the nutrient concentration. We will discuss the constitutive models for the parameters S,τ,IAD in the next section.

At the inlet boundaries (Γ_{c_up,in} e Γ_{c_down,in}) a known concentration has been imposed, using independent values on each inlet section. A homogeneous Neumann condition ∇C_f ∙ n = 0 has been adopted on the bioreactor wall and outlets (). In fact, the wall is considered impermeable to nourishments and their flux in the direction normal to the outlets is assumed equal to zero. Moreover, conservation of concentrations C_f = C_s and of biomolecules flux −D_f∇ ∙ C_fn = −D_s∇ ∙ C_sn have been applied at the interface between fluid and porous media (Γ_{fluid–porous}). As a result, the concentration of oxygen is determined by the following problem: (5)

2.4 Model parameters and constitutive laws

2.4.1 Model parameters for the flow model.

First, the characteristic Reynolds number of the flow in the bioreactor was determined from the following definition, (6) where D is the inlet diameter of 1mm, ρ = 999,97 Kg/m³ and μ = 0,001 Pa ∙ s are the fluid density and dynamic viscosity, respectively, is the inlet flow rate into each chamber, equal to 1 ml/day. A Re ≪ 0,01, was found thus confirming that the assumption of laminar flow is accurately verified. As a consequence, the inertial (and nonlinear) term in the momentum equation, namely , can be neglected and the flow model turns out to be a set of linear equations. This will be the key property for the later derivation of a surrogate of the flow model, which is only based on algebraic equations consequently featuring a negligible computational cost.

Another parameter, essential to determining the flow in the porous medium is the (intrinsic) permeability K_perm that is determined by the microscopic structure of the scaffold, quantified by the porosity (γ), the tortuosity, etc. In the case of materials featuring an anisotropic structure, permeability is a tensor quantity. Here, since the scaffolds under consideration are isotropic, it becomes a scalar parameter. In what follows, we will consider two types of scaffolds, one made out of methacrylated gelatin (GelMA) and the other consisting of a poly-L-lactate (PLLA) foam. The porosity and permeability of the latter have been estimated via Boyle’s pycnometer and scanning electron microscopy (SEM) analysis. Data for GelMA are scarce in literature. However, for tissue engineering it is used as a surrogate material to mimic the extracellular matrix of cartilage; hence, we initialized the model for the bioreactor configuration using data that have been previously measured for native cartilage [26]. In both cases, the values for porosity and permeability are reported in Table 1.

Download:

Table 1. Porosity and permeability values used for GelMA and PLLA scaffolds.

https://doi.org/10.1371/journal.pone.0162774.t001

2.4.2 Model Parameters and constitutive laws for mass transport.

Inlet concentrations for oxygen are 3.15 e-3 [mg/ml] and 7.2 e-3 [mg/ml] for the upper and lower chamber, respectively. We observe that the oxygen supply of the upper chamber falls within the range of hypoxic conditions, compatible with the biological need of the chondral tissue, while the lower chamber, where bone is developed, is kept under normoxic conditions. These different environments are aimed at supporting stem cell differentiation into a chondral and osseous phenotype, respectively [27]. The diffusion coefficient was obtained from previously published studies [28].

For the exchange of biomolecules between fluid and solid phases within the scaffold, the coefficients τ,IAD must be calculated. To this purpose, we model the porous medium as a periodic structure whose unit can be idealized as a cube containing a hollow sphere, namely the pore, as illustrated in Fig 3.

Download:

Fig 3. Representative SEM micrograph of the PLLA scaffold and microscopic model of the scaffold pores for the quantification of the exchange between fluid and solid constituents of the porous matrix.

https://doi.org/10.1371/journal.pone.0162774.g003

Although this configuration is incompatible with the flow through the pore, as it is completely closed, it is adequate for modeling mass transfer between the solid and the fluid phases of the porous medium. According to this model, we estimate the value of the interface area density (IAD), which only depends on the configuration of the unit. Let S_{e_s} = 4πR² and S_{i_s} = 4π(R−δ)² be the external and internal pore surface, respectively, and let V_c be the total volume of the unit. Then the interface area density is defined as: (7)

To estimate the mass transfer coefficient, we assume that at the pore scale mass transfer is dominated by diffusion in the solid phase. As a consequence, the Sherwood number magnitude turns out to be in the range of unity. Exploiting this assumption, we have (8) where τ is the mass transfer coefficient and d is the pore diameter. As a result, we obtain, (9)

We observe that GelMA and PLLA have different behaviors with respect to mass transfer and interface area density. GelMA scaffold has homogeneous properties, namely the pore radius is uniform everywhere and equal to R = 9.77205 e-6[m] with a thickness δ = 10%R. The GelMA matrix is permeable to solutes, as shown by the positive diffusion coefficients D_s,s reported in Table 2. The PLLA scaffold is substantially different because it is impermeable to solutes. As a result, the mass transfer coefficient is necessarily null. Since the exchange between solid and fluid phases in the porous medium is modeled by terms τIAD(C_s,s−C_s,f), we notice that the interface area density does not affect the model.

Download:

Table 2. Oxygen parameters adopted for GelMA and PLLA scaffolds.

https://doi.org/10.1371/journal.pone.0162774.t002

In order to complete the mass transport model, we introduced the term S, to account for both catabolite production and metabolite consumption in cell metabolism. Given the importance of maintaining cell viability by ensuring sufficient nutrients supply, we focus in particular on metabolite consumption, for which studying transport of oxygen is ideal. Cells are assumed to be confined in the porous scaffold and consumption of nutrients, S(C_s), is expected to be proportional to their availability, namely S(C_s) = S(γC_s,f + (1−γ)C_s,s). Different models can be adopted for this function, either linear or nonlinear. In the former case we set S(C_s) = rC_s, where r is a constant parameter determined according to the following balance law: (10) where is a reference concentration for each solute, measured in [mol/ml], V_max is the maximal consumption rate for the considered nutrient and for a specific cell phenotype, quantified in [mol/cell s], and N_v is the average volumetric cell density in the scaffold, measured in [cells/ml]. The main limitation of this model is that it does not guarantee any upper bound for nutrient consumption rate. The more nutrients are available, the more they are metabolized. This approach can be improved using a Michaelis-Menten description of cell metabolism [29], which introduces saturation of the consumption rate, according to the following function: (11) where K_m is the Michaelis-Menten constant, equal to the concentration at which the consumption rate reaches 50% of the maximal value. As a result, the consumption term turns out to be a nonlinear function, namely (12)

We observe that for small nutrient concentrations the linear and the Michaelis-Menten models behave similarly, whereas the latter provides a better estimate of metabolic consumption in case of abundance of nutrients.

2.5 Computational solvers

The commercial code ANSYS CFX v.13.0 was used to carry out the fluid dynamic and mass transport simulations. The spatial discretization consists of a cell based finite volume method.

From the computational standpoint, the main challenge of this study consists in solving a fluid-porous interaction problem that involves coupled flow and mass transport. A fully coupled strategy has been adopted, namely all the equations are solved simultaneously through a monolithic linear system that embraces all the degrees of freedom.

More precisely, the Laplace operator in the fluid momentum and oxygen transport equations is approximated by a centered scheme, while the convective terms have been discretized by means of an upwind method. The convective term in the Navier-Stokes equations is linearized by Picard iterations (equivalent to a fictitious time stepping method with semi-implicit treatment of ) (“ANSYS CFX-Solver Theory Guide”, ANSYS Inc., 2010). The pressure variable in the Navier-Stokes equations is evaluated at the same nodes of the velocity field.

The system is then solved using an algebraic multigrid method exploiting incomplete LU factorization as smoother. Numerical simulations have been performed on parallel CPUs using a quad-socket 12-Core AMD Magny Cours CPU, 128 GB RAM at University of Pittsburgh. Convergence criteria were set to 10⁻⁶ for the normalized residuals of the global linear system of equations.

To ease the convergence of the algebraic solver, it turned out to be extremely helpful to neglect the contribution of streamline diffusion in the mass transport model, accounting only for the cross-wind component of the diffusion operator. From the modeling standpoint, this approximation is justified since the Péclet number characterizing mass transport in the ducts and in the scaffold of the bioreactor is larger than unity. More precisely, we define the Péclet number as follows (13) where a is the characteristic length of diffusion, is the characteristic fluid velocity and D is the diffusion coefficient of the nutrient in the fluid (water). The Péclet number has been calculated for two sets of parameters, the first one identifying flow and mass transport in the pores of the insert (a = 9.77205 e-6 m, = 1.546e-3 m/s, D = 2.9e-9 m²/s) and the second one the flow in the chambers that will hold the scaffold (a = 5 e-4 [m], = 1.473e-5 [m/s], D = 2.9e-9 [m²/s]). For the insert we obtained Pe = 2.5, while for the chambers Pe = 5.2.

Domain discretization is a crucial phase in the computational model set up to ensure an accurate description of the investigated phenomena as well as reasonable computational time and costs. The geometrical features of the bioreactors span from 8.5 mm (height of the scaffold), to 1 mm (inlet/outlet channel inner diameter), to 0.25 mm (radius of the pores). The final mesh consists of 735658 and 550226 tetrahedral elements for the GelMA and the PLLA case, respectively, with a minimum dimension of the elements of 0.1 mm and a maximum of 0.25 mm. This discretization is suitable for the fluid dynamics model, because, as previously stated, the Reynold’s number results smaller than 0.01, and consequently the boundary layers can be considered fully developed. The fluid dynamics simulations in single array are performed with moderate computational effort (about 7 minutes on CPUs using a quad-socket 12-Core AMD Magny Cours CPU, 128 GB RAM). A numerical test that uses a coarser mesh consisting of 443740 and 242236 elements, respectively, confirms that the results obtained with the finer discretization are insensitive to the mesh size.

3 Lumped Parameter Models of HTB

Although in-silico analysis is rightfully considered a cost efficient approach with respect to experimental investigation, section 2.3 illustrates that the development of a computational model of the bioreactor is a challenging task, because of the significant amount of work-hours required to define a detailed CAD model and the considerable computational efforts involved with the definition of a computational mesh and with the solution of the discrete equations.

When using numerical tools in the design or optimization of the bioreactor configuration and working conditions, it is essential to minimize the cost of running simulations for different sets of design parameters. The scientific computing community is well aware of this critical aspect of the approach and has recently made great progress in developing strategies to synthesize surrogate models that replace the brute force simulation approach with much less computational costs. We have mentioned a list of a few examples related to bioengineering [30–36], among many others.

Surrogate or reduced models are based on much simpler mathematical operators than partial differential equations. For steady problems, they may consist of algebraic equations, or ordinary differential equations to capture time dependent phenomena. Such models are often called lumped parameter models, because they synthesize into a small number of coefficients the behavior of spatially dependent functions, solutions of partial differential equations, a.k.a. distributed parameter models.

The aim of this section is to derive a set of lumped parameter models describing flow and mass transport in the bioreactor fulfilling two objectives:

To determine the change of quantitative outputs when the input data are varied, for a fixed single or multi-chamber configuration,
To determine the change of quantitative outputs when the number of chambers in the array is varied.

3.1 Lumped parameter model for a fixed HTB configuration

We aim to develop an input-output relation between parameters of the model and observed quantities of interest. Because of the linearity of the flow model, motivated by low Reynolds numbers, this relation is a linear operator that can be characterized by a limited number of simulations. The number of required simulations depends on the dimension of the input/output parameter space.

To illustrate the derivation of a lumped parameter model, we consider an example that will be later used for the bioreactor design. In particular, we analyze the flow split at the outlet of the bioreactor chambers for prescribed values of the inlet flow rates.

Let us consider the velocity fields defined by fixing unit flow rates at each inlet of the bioreactor, (14)

Since the flow model is linear, the velocity and pressure fields corresponding to any combination of the inlet flow rates, denoted as Q_in1 and Q_in2 respectively, can be represented as a linear combination of solutions (15)

Since we are interested in the quantification of the outflow rates, we calculate (16)

As a result, we have identified the following input-output algebraic relation between inlet and outlet flow rates (17) that represents the lumped parameter model we were looking for. We note that the operator (matrix) M depends on the bioreactor geometric design.

This approach can be extended to the mass transport problem, provided that the model adopted for consumption of nutrients is linear, namely S(C_s) = rC_s. In this case, we denote with d_i the solution of Eq 5 obtained setting and . Then, any solution C_f of the mass transport problem can be expressed as (18)

Let C_out,1,C_out,2 be the nutrient concentration on the upper and lower outlets respectively and for simplicity of notation let us define (19)

Then, because of the linearity of the mass transport model we obtain (20) that can be translated in the following vector form, (21)

3.2 Lumped parameter model for variable bioreactor configurations

Here we focus on the problem of determining a lumped parameter model for a sequence of bioreactor units, when the solution for 1-unit is known. From the methodological standpoint, this problem is more challenging than the one of characterizing the lumped parameter model for one bioreactor unit, because partial differential equations are not linear with respect to the configuration of the domain. In other words, the solution of an n-unit bioreactor is not the superposition of n solutions of a single unit configuration.

Another strategy for determining a lumped parameter model of a multi-unit configuration emerges observing that units are combined in sequence (see Fig 4). Consequently, we conjecture that the behavior of the n-unit bioreactor is the composition of n-unit models. As an example, for a sequence of two units we posit that the input/output relation for flow rates is (22)

Download:

Fig 4. Distribute vs lumped parameter models

Top: A 8-unit bioreactor configuration, showing details of a 2-unit example used for the development of the lumped parameter model (top panels). Bottom: A sketch of a multi-unit bioreactor configuration with heterogeneous unit design in a generic sequence of units, where different unit designs are denoted with letters A, B, C.

https://doi.org/10.1371/journal.pone.0162774.g004

Owing to the similar design of the upper and lower chambers, the resistance to flow of the fluid entering from the upper and lower inlets is comparable. As a result, the following property is valid at any junction between two adjacent bioreactor units, (23)

It shows that equal normal stresses are applied at the intermediate section of a 2-unit bioreactor. Since these are the boundary conditions applied at the outlet of our model for an individual unit it means that any unit in a row functions as an individual one. As a result, we conclude that (24) and consequently (25)

This example can be easily generalized to the case of a row of n-units. More precisely, we infer that the lumped parameter model for an n-unit bioreactor, denoted by M_n is the multiplicative composition of n single unit models, namely (26) where the latter expression denotes the n-th power of the operator M.

This approach can be applied to flow (as illustrated above) as well as to mass transport. In this way, the lumped parameter models M,D, derived in section 3.1 for single unit configurations, can be extended to multi-unit configurations made of units combined in a row. Using direct numerical simulations of multi-cell configurations, we will demonstrate in the next sections the good accuracy of these reduced models.

We finally observe that the model composition rule is also applicable in the case of combination of different unit designs (schematized in Fig 4 with letters A, B, C). In particular, the input/output relation (Y = M ∙ X) for a row of 3-units of generic type A, B, C of which we know the individual lumped parameter models, M_A,M_B,M_C respectively, is given by M = M_A ∙ M_B ∙ M_C. Following the ambitious vision of building a human-on-chip model, any pattern of bioreactors organized in a row can be characterized using this approach, provided that the properties of each individual unit are known.

4 Numerical Simulations

4.1 Numerical simulation of flow

In this study, simulations of flow are performed to compare flow patterns in the GelMA and PLLA scaffold when inlet flow rates are varied. More precisely, the following different flow pairs were simulated: (a) 1 and 1, (b) 1 and 2 and (c) 10 and 10 ml/day for the upper and lower inlet, respectively.

We observe that for all the configurations, the fluid is driven by the pressure gradient to move toward the upper chamber (Fig 5). The flow split obtained by applying the different flow pairs are reported in Tables 3 and 4 for the GelMA and PLLA case, respectively. The comparison of the outlet flow rates for the two scaffolds highlighted opposite outcomes in terms of flow mixing. Indeed, while not significant flow mixing was found for the GelMA scaffold, a significant mixing occurs in the PLLA case. As expected, the maximum mixing (that is 42.9%) occurs with different input fluid flow rates (1 and 2 ml/day at the upper and lower inlet, respectively).

Download:

Fig 5. Streamlines in the 1- unit model with the GelMA scaffold.

https://doi.org/10.1371/journal.pone.0162774.g005

Download:

Table 3. Results obtained by simulating different flow split in the one unit model with the GelMA and PLLA scaffold.

https://doi.org/10.1371/journal.pone.0162774.t003

Download:

Table 4. Percentage of oxygen consumption for the GelMA and PLLA scaffold.

https://doi.org/10.1371/journal.pone.0162774.t004

For the sake of brevity, the results of the 4-units array are not reported since they are the qualitatively equivalent to the single unit configuration.

4.2 Numerical simulation of transport

Simulations of oxygen transport were performed to compare mass transfer in the GelMA and PLLA scaffolds.

Concentrations equal to 3.15 and 7.2 μg/l were applied at the upper and lower inlet, respectively. As in the previous case, the following flow pairs were simulated: (a) 1 and 1, (b) 1 and 2 and (c) 10 and 10 ml/day at the upper and lower inlet. Two configurations of the bioreactor were considered, namely 1-unit and a 4-unit array. The results of 1-unit model are reported in Fig 6 and Fig 7 for the GelMA and PLLA scaffold, respectively. The analysis of the mass transport simulations obtained for the GelMA and the PLLA scaffolds allows us to draw general considerations, which are valid for both single and 4-unit arrays.

Download:

Fig 6. Oxygen concentration with GelMA scaffold.

From left to right, flow pair of 1–1 [ml/day], 1–2 [ml/day], 10–10 [ml/day].

https://doi.org/10.1371/journal.pone.0162774.g006

Download:

Fig 7. Oxygen concentration with PLLA scaffold.

From left to right, flow pair of 1–1 [ml/day], 1–2 [ml/day], 10–10 [ml/day].

https://doi.org/10.1371/journal.pone.0162774.g007

Firstly, as explained in section 2.3.2, we see that axial advection is dominant with respect to the cross-wind diffusion. Therefore, the higher the flow rates and fluid velocity, the more the inlet and outlet oxygen concentrations look similar due to a reduced oxygen drop (Fig 6B and Fig 7B). However, the diffusion of oxygen from the lower chamber to the upper one is not negligible, because different inlet concentrations promote the formation of concentration gradients that trigger transport.

For both the GelMA and PLLA cases, the oxygen concentration in the top region of the scaffold is higher in the case of low flow rate, (a, inlet flow equal to 1 ml/day) than in the case of high flow rate (c, inlet flow equal to 10 ml/day). Concerning case (b), the mix of the two chambers’ flow is greater and a contribution of convective transport is added to the diffusive flux from the bottom towards the top of the bioreactor chamber. For this reason, the oxygen concentration in the top region of the scaffold is greater in case (b) than in cases (a) and (c).

Finally, the simulations suggest that the scaffold porosity and permeability play a relevant role on mass transport. Indeed, while the GelMA is permeable to oxygen, the PLLA is not. This implies that the aforementioned phenomena are more evident with a polymeric scaffold impervious to mass transport through the solid phase, such as PLLA.

4.3 Oxygen consumption

The simulations of oxygen consumption were performed for the two different scaffolds (GelMA and PLLA) for an array of 4-units, in order to study the depletion of nutrients in the culture medium. The flow split is the one of case (a) (1 and 1 ml/day) and the inlets concentrations are equal to 3.15 and 7.2 μg/l at the upper and lower inlet, which correspond to the normoxic levels of the different types of tissue grown in the upper and lower chambers.

Since we consider a 4-unit array, we observe that diffusion develops more easily along the bioreactor axis (longer fluid path with respect to the 1-unit case) and as a consequence, the oxygen concentration tends to become more uniform. More precisely, enhanced diffusion combined with different inlet concentrations causes a decrease of the oxygen level in the lower chamber and an increase in the top one. This trend is heightened by cellular oxygen consumption, which further leads to a diminishing of the oxygen concentration in the lower chamber (Fig 8).

Download:

Fig 8. Oxygen concentration in the 4 cells array.

GelMA (left) and PLLA (right) scaffolds are used when with active consumption rate.

https://doi.org/10.1371/journal.pone.0162774.g008

The two types of scaffold show the same trend of oxygen consumption, but the computations highlighted different percentage of consumed oxygen. Indeed, a higher percentage of oxygen consumption was found for the PLLA scaffolds with respect to GelMA. This effect is likely a result of the different cell density used for the two cases. In fact, cell density is assumed to be equal to 1 x 10⁶ cells/ml for GelMA and to 2.12304019 x 10⁶ cells/ml in the case of PLLA.

4.4 Comparison of distributed and lumped parameter models

In this section, the results of the lumped and the distributed parameter models are presented and compared in terms of fluid dynamics and mass transport. The fluid dynamics results for 1-unit and 4-units array are first presented, then, the mass transport results of both configurations are studied. For the sake of brevity, we present only the results obtained by simulating the GelMA scaffold.

4.4.1 Fluid dynamics.

Two computational fluid dynamics simulations were performed for the single unit configuration to determine the lumped parameter model (LPM). In particular, two inlet flow pairs are applied as reported in Table 5. The resulting LPM matrix M is:

Download:

Table 5. Simulation settings to identify the fluid dynamics characteristics of one-unit bioreactor.

https://doi.org/10.1371/journal.pone.0162774.t005

Then, the results of the 1-unit and 4-unit LPMs are compared to those of the distributed parameter model, see Tables 6 and 7, and in two test cases the error was lower than 1%.

Download:

Table 6. Comparison of the 1-unit fluid dynamics results provided by the distributed (distr) and the lumped (lump) parameter models.

https://doi.org/10.1371/journal.pone.0162774.t006

Download:

Table 7. Comparison of the 4-unit array fluid dynamics results provided by the distributed (distr) and the lumped (lump) parameter models.

https://doi.org/10.1371/journal.pone.0162774.t007

4.4.2 Mass transport.

For the LPM model of mass transport we have adopted the parameters of Table 2 and inlet concentrations summarized in Table 8.

Download:

Table 8. Simulation settings to identify the mass transport input-output characteristics of one-unit bioreactor.

https://doi.org/10.1371/journal.pone.0162774.t008

To start with, we analyze the mass transport model without cell metabolism, that is the case S(C_s) = 0 in Eq 5. The LPM model for the corresponding mass transport simulations is the following matrix:

The results of the 1-unit LPM are compared with those of the distributed parameter model in two simulations with different inlets concentrations, reported in Table 9, whose values are set according to ongoing experimental tests. The results from the LPM model differ from those of the distributed parameters model by less than the 1%.

Download:

Table 9. Comparison of the one-unit oxygen concentration results provided by the distributed (distr) and the lumped (lump) parameter models.

https://doi.org/10.1371/journal.pone.0162774.t009

We also calculate the LPM model for mass transport with active cell metabolism. For the linear model, S(C_s) = rC_s, the LPM matrix for 1-unit is the following while for the Michaelis-Menten case, namely Eqs 11 and 12, the LPM model becomes

The inspection of the matrices D,D_l,D_mm informs about the characteristics of the different consumption models compared here. We observe that the diagonal entries of D_l are the smallest, confirming that the linear model is the one with the highest oxygen consumption rate. The extra-diagonal coefficients correspond to the oxygen exchange between the upper and lower chambers. Their magnitude is similar in all cases, because they depend on the diffusion parameters solely. For the linear case, the theory at the basis of the LPM derivation is satisfied, while it does not rigorously hold true for the Michaelis-Menten model, because the mass transport equation becomes nonlinear. Once again, numerical simulations based on the full model applied to the 8-unit array confirm that the LPM model with linear consumption rate, namely D_l, predicts outlet concentrations with less than 1% error. The corresponding results are reported in Table 10 and visualized in Fig 9. In Table 11 we report the error of the LPM based on the Michaelis-Menten nonlinear consumption rate. Despite the nonlinear nature of the problem, in conflict with the principles at the basis of the LPM derivation, the LPM model is fairly accurate in predicting the concentration split and decay at the outlet also with a Michaelis-Menten consumption rate, with a maximum error of about 10% for an array of 4-units, located on the bottom outlet of the bioreactor.

Download:

Fig 9. Variation of the outlet concentration of oxygen with respect to the number of units (unit #0 denotes the inlet value) for the mass transport model without cell metabolism (dashed line), with linear consumption rate (dotted line) and with Michaelis-Menten consumption model (solid line).

Data calculated using the full 3D model are reported in red.

https://doi.org/10.1371/journal.pone.0162774.g009

Download:

Table 10. Comparison of the 8-unit array oxygen concentration results provided by the distributed (distr) and the lumped (lump) parameter models with linear consumption rate.

https://doi.org/10.1371/journal.pone.0162774.t010

Download:

Table 11. Comparison of the 4-unit array oxygen concentration results provided by the distributed (distr) and the lumped (lump) parameter models with Michaelis-Menten consumption rate.

https://doi.org/10.1371/journal.pone.0162774.t011

The LPM model for mass transport is particularly interesting because it allows us to estimate the decay of nutrient concentrations due to cell metabolism along an arbitrarily long array of units, using the formula . Considering for example the inlet concentrations of Table 9, test case #2 for we estimate the outlet concentration decay for the transport model without oxygen consumption. The same calculation is then repeated for the linear and the Michaelis-Menten models for cell metabolism and the results are compared in Fig 9, where also the outlet concentrations determined using the fully 3D simulations are shown for a qualitative visualization of the LPM error.

5 Discussion

From the engineering standpoint, our study shed lights on important aspects of the bioreactor behavior. We observe that the flow is dominated by viscous effects and by pressure gradients, while inertial effects are negligible. Differences in inlet velocities between upper and lower chamber generate a vertical pressure gradient inside the bioreactor chambers, which promotes mixing of nutrient fluid flowing through the osteochondral construct. Furthermore, we have observed that the magnitude of vertical pressure gradients depends highly on the permeability of the scaffold. Between the two materials tested here, it appears that the most permeable one favors the mixing of fluid among the upper and lower chambers.

Concerning mass transfer, our simulations suggest that it is dominated by convection. Diffusion effects are however non-negligible, but their (relative) intensity varies according to the inlet flow rate and the scaffold properties. More precisely, Fig 6 and Fig 7 show that high flow rates decrease the transport of biochemical species between the two chambers. From the analysis of these plots we also observe that the concentration in the bioreactor top chamber is greater than the one at the upper outlet. This means that the exchange between the chamber and the supplying channels is not sufficient to remove all the chemical species that accumulate in this region, because of combined diffusion and convection. This effect is observable for both types of scaffold, but is more evident for GelMA, suggesting that this type of material hinders flow and mass transport more than PLLA does. When nutrient (or oxygen) consumption is switched on in the simulation, concentration gradients are quickly smoothed out when traveling along multiple bioreactor units. At the same time, concentration levels significantly decrease. The computational model thus serves as a valuable tool to estimate whether the final units of the row receive enough nutrients, as illustrated in the example presented above.

Finally, we have developed a surrogate, inexpensive approach to characterize the output of the bioreactor without the burden of running many computer simulations. It consists of a lumped parameter model, derived exploiting the linearity of the full model. The LPM has proven to be very accurate in capturing the effect of sequentially combining multiple units. A natural application of this model is studying the concentration decay along a sequence of bioreactor units. For example, Fig 9 shows the concentration decay at the bioreactor outlets when the number of units is varied from 1 to 16. Three sets of curves outline the behavior of different cell metabolism models. When cell metabolism is switched off (dashed lines), the upper and lower concentrations equilibrate very quickly, confirming that diffusion effects of oxygen between the two chambers are non negligible. We recall that large oxygen diffusion and transport between the upper and lower chambers is not necessarily desirable, when different types of tissue are grown. Indeed, in our case, cartilage natural environment should be hypoxic, while bone better develops in normoxic conditions. For constant consumption rate, the concentration decay is the largest. As a consequence after 16 bioreactor units, almost all the nutrient concentration has been consumed. The Michaelis-Menten metabolic model is the most realistic of the three options, because it accounts for a saturation effect that limits the consumption rate. According to our preliminary data on cell viability in the bioreactor, obtained by Live/Dead assays (data not shown), the oxygenation computed after 16-units appears to be still at a sufficient level.

The computational approach proposed here is subject to some limitations. One is the approximation of the fluid dynamic and mass transport through steady model. A key challenge in the engineering of three-dimensional tissue is maintenance of cell viability when the volumetric cell density increases. In this study, we assumed a constant cell density equal to the initial culture conditions that occur after distributing cells homogenously throughout the volume of the scaffolds. However, variations in cell density with time could be easily incorporated in both our models, to predict oxygen drops in long-term culture. Secondly, as literature data are lacking, we assumed the GelMA properties (i.e., porosity and permeability) equal to those of native cartilage. Experimental tests will be performed in future work to assess these properties. Finally, we have not accounted for transport along capillaries. This could be acceptable for many engineered constructs that are approximation of native tissues, frequently obtained from single cell types, e.g., mesenchymal stem cells, within a hydrogel or a porous scaffold. If the HTB were to be used with native tissues, we expect our approach to hold true with the necessary adjustments to account for the different tissues types. The avascular components of cartilage would be modeled adjusting the parameters we currently used for GelMA, whereas the for the vascularized bone, the more porous structure we described for the PLLA scaffold could offer a good starting model to approximate the cavities and capillaries present in subchondral bone.

Another improvement of our study would be to validate the oxygen concentration drops predicted by our models with actual measurements performed when the bioreactor is operated with cell-seeded constructs. This validation would be technically challenging, only feasible using oxygen sensors incorporated in the perfusion circuit, at the inlet and outlets of each bioreactor unit or even inserted directly in the chambers, in direct contact with the living cells [37]. Detecting larger molecules, even at low concentrations provides a more simple and reliable quantification. On this basis, extensive validation of the ability of our models to predict the flow-dynamics and mass transport in the bioreactor will be the subject of future work.

6 Conclusions and Perspectives

From the methodological standpoint, we have overcome the challenge of developing a complex multi-physics model of the bioreactor. We have also succeeded in implementing the model into a commercial computational platform, showing the significant potential of computational tools on biomedical research, including analytical cases integrating quantitative biology and translational medicine. Future developments of this study consist of experimental validation of the models and their application to explore different bioreactor configurations. Such findings will allow optimization of the model by incorporating the multi-faceted factors that affect its behavior and functionality.

Author Contributions

Conceptualization: LI GDU GC EC RST MTR RG PZ.
Data curation: LI GDU EC PZ.
Formal analysis: LI GDU EC PZ.
Funding acquisition: RST MTR RG PZ.
Investigation: LI GDU GC EC RST MTR RG PZ.
Methodology: LI GDU GC EC RST MTR RG PZ.
Project administration: RG PZ.
Resources: EC PZ.
Software: LI GDU EC PZ.
Supervision: RST MTR RG PZ.
Validation: LI GDU GC RST RG.
Visualization: LI GDU GC EC RG PZ.
Writing – original draft: LI GDU GC EC RST MTR RG PZ.
Writing – review & editing: LI GDU GC EC RST MTR RG PZ.

References

1. Alexander PG, Gottardi R, Lin H, Lozito TP, Tuan RS. Three-dimensional osteogenic and chondrogenic systems to model osteochondral physiology and degenerative joint diseases. Experimental Biology and Medicine. 2014;239(9):1080–95. pmid:24994814
- View Article
- PubMed/NCBI
- Google Scholar
2. Eghbali H, Nava MM, Mohebbi-Kalhori D, Raimondi MT. Hollow fiber bioreactor technology for tissue engineering applications. The International journal of artificial organs. 2016:0. Epub 2016/02/27. pmid:26916757.
- View Article
- PubMed/NCBI
- Google Scholar
3. Perestrelo A, Águas A, Rainer A, Forte G. Microfluidic Organ/Body-on-a-Chip Devices at the Convergence of Biology and Microengineering. Sensors. 2015;15(12):29848.
- View Article
- Google Scholar
4. Giusti S, Sbrana T, La Marca M, Di Patria V, Martinucci V, Tirella A, et al. A novel dual-flow bioreactor simulates increased fluorescein permeability in epithelial tissue barriers. Biotechnology Journal. 2014;9(9):1175–84. pmid:24756869
- View Article
- PubMed/NCBI
- Google Scholar
5. Loskill P, Marcus SG, Mathur A, Reese WM, Healy KE. μorgano: A Lego®-like plug & play system for modular multi-organ-chips. PLoS ONE. 2015;10(10).
- View Article
- Google Scholar
6. Lin H, Lozito TP, Alexander PG, Gottardi R, Tuan RS. Stem cell-based microphysiological osteochondral system to model tissue response to interleukin-1Β. Molecular Pharmaceutics. 2014;11(7):2203–12. pmid:24830762
- View Article
- PubMed/NCBI
- Google Scholar
7. Lozito TP, Alexander PG, Lin H, Gottardi R, Cheng AWM, Tuan RS. Three-dimensional osteochondral microtissue to model pathogenesis of osteoarthritis. Stem Cell Research and Therapy. 2013;4(SUPPL.1).
- View Article
- Google Scholar
8. Lawrence RC, Felson DT, Helmick CG, Arnold LM, Choi H, Deyo RA, et al. Estimates of the prevalence of arthritis and other rheumatic conditions in the United States. Part II. Arthritis and Rheumatism. 2008;58(1):26–35. pmid:18163497
- View Article
- PubMed/NCBI
- Google Scholar
9. Stolz M, Gottardi R, Raiteri R, Miot S, Martin I, Imer R, et al. Early detection of aging cartilage and osteoarthritis in mice and patient samples using atomic force microscopy. Nature Nanotechnology. 2009;4(3):186–92. pmid:19265849
- View Article
- PubMed/NCBI
- Google Scholar
10. Elder SH, Goldstein SA, Kimura JH, Soslowsky LJ, Spengler DM. Chondrocyte differentiation is modulated by frequency and duration of cyclic compressive loading. Annals of Biomedical Engineering. 2001;29(6):476–82. pmid:11459341
- View Article
- PubMed/NCBI
- Google Scholar
11. Wikswo JP. The relevance and potential roles of microphysiological systems in biology and medicine. Experimental Biology and Medicine. 2014;239(9):1061–72. pmid:25187571
- View Article
- PubMed/NCBI
- Google Scholar
12. Rouwkema J, Gibbs S, Lutolf MP, Martin I, Vunjak-Novakovic G, Malda J. In vitro platforms for tissue engineering: Implications for basic research and clinical translation. Journal of Tissue Engineering and Regenerative Medicine. 2011;5(8):e164–e7. pmid:21774080
- View Article
- PubMed/NCBI
- Google Scholar
13. Goldring MB, Goldring SR. Articular cartilage and subchondral bone in the pathogenesis of osteoarthritis. Annals of the New York Academy of Sciences2010. p. 230–7.
- View Article
- Google Scholar
14. Suri S, Walsh DA. Osteochondral alterations in osteoarthritis. Bone. 2012;51(2):204–11. pmid:22023932
- View Article
- PubMed/NCBI
- Google Scholar
15. Yuan XL, Meng HY, Wang YC, Peng J, Guo QY, Wang AY, et al. Bone-cartilage interface crosstalk in osteoarthritis: Potential pathways and future therapeutic strategies. Osteoarthritis and Cartilage. 2014;22(8):1077–89. pmid:24928319
- View Article
- PubMed/NCBI
- Google Scholar
16. Martin I, Miot S, Barbero A, Jakob M, Wendt D. Osteochondral tissue engineering. Journal of biomechanics. 2007;40(4):750–65. Epub 2006/05/30. pmid:16730354.
- View Article
- PubMed/NCBI
- Google Scholar
17. Giannoni P, Lazzarini E, Ceseracciu L, Barone AC, Quarto R, Scaglione S. Design and characterization of a tissue-engineered bilayer scaffold for osteochondral tissue repair. J Tissue Eng Regen Med. 2015;9(10):1182–92. Epub 2012/11/23. pmid:23172816.
- View Article
- PubMed/NCBI
- Google Scholar
18. Tampieri A, Sandri M, Landi E, Pressato D, Francioli S, Quarto R, et al. Design of graded biomimetic osteochondral composite scaffolds. Biomaterials. 2008;29(26):3539–46. Epub 2008/06/10. pmid:18538387.
- View Article
- PubMed/NCBI
- Google Scholar
19. Yang G, Rothrauff BB, Lin H, Gottardi R, Alexander PG, Tuan RS. Enhancement of tenogenic differentiation of human adipose stem cells by tendon-derived extracellular matrix. Biomaterials. 2013;34(37):9295–306. pmid:24044998
- View Article
- PubMed/NCBI
- Google Scholar
20. Candiani G, Raimondi MT, Aurora R, Lagana K, Dubini G. Chondrocyte response to high regimens of cyclic hydrostatic pressure in 3-dimensional engineered constructs. The International journal of artificial organs. 2008;31(6):490–9. Epub 2008/07/09. pmid:18609501.
- View Article
- PubMed/NCBI
- Google Scholar
21. Mazzei D, Guzzardi MA, Giusti S, Ahluwalia A. A low shear stress modular bioreactor for connected cell culture under high flow rates. Biotechnology and Bioengineering. 2010;106(1):127–37. pmid:20091740
- View Article
- PubMed/NCBI
- Google Scholar
22. Sbrana T, Ucciferri N, Favrè M, Ahmed S, Collnot EM, Lehr CM, et al. Dual flow bioreactor with ultrathin microporous TEER sensing membrane for evaluation of nanoparticle toxicity. Sensors and Actuators, B: Chemical. 2016;223:440–6.
- View Article
- Google Scholar
23. Vozzi F, Mazzei D, Vinci B, Vozzi G, Sbrana T, Ricotti L, et al. A flexible bioreactor system for constructing in vitro tissue and organ models. Biotechnology and Bioengineering. 2011;108(9):2129–40. pmid:21495015
- View Article
- PubMed/NCBI
- Google Scholar
24. Lin H, Zhang D, Alexander PG, Yang G, Tan J, Cheng AWM, et al. Application of visible light-based projection stereolithography for live cell-scaffold fabrication with designed architecture. Biomaterials. 2013;34(2):331–9. pmid:23092861
- View Article
- PubMed/NCBI
- Google Scholar
25. Mannella GA, Conoscenti G, Carfì Pavia F, La Carrubba V, Brucato V. Preparation of polymeric foams with a pore size gradient via Thermally Induced Phase Separation (TIPS). Materials Letters. 2015;160:31–3.
- View Article
- Google Scholar
26. Taffetani M, Gottardi R, Gastaldi D, Raiteri R, Vena P. Poroelastic response of articular cartilage by nanoindentation creep tests at different characteristic lengths. Medical Engineering and Physics. 2014;36(7):850–8. pmid:24814573
- View Article
- PubMed/NCBI
- Google Scholar
27. Nava MM, Raimondi MT, Pietrabissa R. Controlling self-renewal and differentiation of stem cells via mechanical cues. Journal of Biomedicine and Biotechnology. 2012;2012.
- View Article
- Google Scholar
28. Sacco R, Causin P, Zunino P, Raimondi MT. A multiphysics/multiscale 2D numerical simulation of scaffold-based cartilage regeneration under interstitial perfusion in a bioreactor. Biomechanics and Modeling in Mechanobiology. 2011;10(4):577–89. pmid:20865436
- View Article
- PubMed/NCBI
- Google Scholar
29. Michaelis L, Menten MML. The kinetics of invertin action: Translated by T.R.C. Boyde Submitted 4 February 1913. FEBS Letters. 2013;587(17):2712–20. pmid:23867202
- View Article
- PubMed/NCBI
- Google Scholar
30. Astorino M, Hamers J, Shadden SC, Gerbeau JF. A robust and efficient valve model based on resistive immersed surfaces. International Journal for Numerical Methods in Biomedical Engineering. 2012;28(9):937–59. pmid:22941924
- View Article
- PubMed/NCBI
- Google Scholar
31. Braakman R, Sipkema P, Westerhof N. A dynamic nonlinear lumped parameter model for skeletal muscle circulation. Annals of Biomedical Engineering. 1989;17(6):593–616. pmid:2589694
- View Article
- PubMed/NCBI
- Google Scholar
32. Formaggia L, Lamponi D, Tuveri M, Veneziani A. Numerical modeling of 1D arterial networks coupled with a lumped parameters description of the heart. Computer Methods in Biomechanics and Biomedical Engineering. 2006;9(5):273–88. pmid:17132614
- View Article
- PubMed/NCBI
- Google Scholar
33. Kim HJ, Vignon-Clementel IE, Figueroa CA, Ladisa JF, Jansen KE, Feinstein JA, et al. On coupling a lumped parameter heart model and a three-dimensional finite element aorta model. Annals of Biomedical Engineering. 2009;37(11):2153–69. pmid:19609676
- View Article
- PubMed/NCBI
- Google Scholar
34. Lee J, Smith NP. The multi-scale modelling of coronary blood flow. Annals of Biomedical Engineering. 2012;40(11):2399–413. pmid:22565815
- View Article
- PubMed/NCBI
- Google Scholar
35. Pant S, Fabrèges B, Gerbeau JF, Vignon-Clementel IE. A methodological paradigm for patient-specific multi-scale CFD simulations: From clinical measurements to parameter estimates for individual analysis. International Journal for Numerical Methods in Biomedical Engineering. 2014;30(12):1614–48. pmid:25345820
- View Article
- PubMed/NCBI
- Google Scholar
36. Shi Y, Lawford P, Hose R. Review of Zero-D and 1-D Models of Blood Flow in the Cardiovascular System. BioMedical Engineering Online. 2011;10.
- View Article
- Google Scholar
37. Raimondi MT, Giordano C, Pietrabissa R. Oxygen measurement in interstitially perfused cellularized constructs cultured in a miniaturized bioreactor. Journal of Applied Biomaterials and Functional Materials. 2015;13(4):e313–e9. pmid:26450634
- View Article
- PubMed/NCBI
- Google Scholar

[ref1] 1. Alexander PG, Gottardi R, Lin H, Lozito TP, Tuan RS. Three-dimensional osteogenic and chondrogenic systems to model osteochondral physiology and degenerative joint diseases. Experimental Biology and Medicine. 2014;239(9):1080–95. pmid:24994814
View Article
PubMed/NCBI
Google Scholar

[2] View Article

[3] PubMed/NCBI

[4] Google Scholar

[ref2] 2. Eghbali H, Nava MM, Mohebbi-Kalhori D, Raimondi MT. Hollow fiber bioreactor technology for tissue engineering applications. The International journal of artificial organs. 2016:0. Epub 2016/02/27. pmid:26916757.
View Article
PubMed/NCBI
Google Scholar

[6] View Article

[7] PubMed/NCBI

[8] Google Scholar

[ref3] 3. Perestrelo A, Águas A, Rainer A, Forte G. Microfluidic Organ/Body-on-a-Chip Devices at the Convergence of Biology and Microengineering. Sensors. 2015;15(12):29848.
View Article
Google Scholar

[10] View Article

[11] Google Scholar

[ref4] 4. Giusti S, Sbrana T, La Marca M, Di Patria V, Martinucci V, Tirella A, et al. A novel dual-flow bioreactor simulates increased fluorescein permeability in epithelial tissue barriers. Biotechnology Journal. 2014;9(9):1175–84. pmid:24756869
View Article
PubMed/NCBI
Google Scholar

[13] View Article

[14] PubMed/NCBI

[15] Google Scholar

[ref5] 5. Loskill P, Marcus SG, Mathur A, Reese WM, Healy KE. μorgano: A Lego®-like plug & play system for modular multi-organ-chips. PLoS ONE. 2015;10(10).
View Article
Google Scholar

[17] View Article

[18] Google Scholar

[ref6] 6. Lin H, Lozito TP, Alexander PG, Gottardi R, Tuan RS. Stem cell-based microphysiological osteochondral system to model tissue response to interleukin-1Β. Molecular Pharmaceutics. 2014;11(7):2203–12. pmid:24830762
View Article
PubMed/NCBI
Google Scholar

[20] View Article

[21] PubMed/NCBI

[22] Google Scholar

[ref7] 7. Lozito TP, Alexander PG, Lin H, Gottardi R, Cheng AWM, Tuan RS. Three-dimensional osteochondral microtissue to model pathogenesis of osteoarthritis. Stem Cell Research and Therapy. 2013;4(SUPPL.1).
View Article
Google Scholar

[24] View Article

[25] Google Scholar

[ref8] 8. Lawrence RC, Felson DT, Helmick CG, Arnold LM, Choi H, Deyo RA, et al. Estimates of the prevalence of arthritis and other rheumatic conditions in the United States. Part II. Arthritis and Rheumatism. 2008;58(1):26–35. pmid:18163497
View Article
PubMed/NCBI
Google Scholar

[27] View Article

[28] PubMed/NCBI

[29] Google Scholar

[ref9] 9. Stolz M, Gottardi R, Raiteri R, Miot S, Martin I, Imer R, et al. Early detection of aging cartilage and osteoarthritis in mice and patient samples using atomic force microscopy. Nature Nanotechnology. 2009;4(3):186–92. pmid:19265849
View Article
PubMed/NCBI
Google Scholar

[31] View Article

[32] PubMed/NCBI

[33] Google Scholar

[ref10] 10. Elder SH, Goldstein SA, Kimura JH, Soslowsky LJ, Spengler DM. Chondrocyte differentiation is modulated by frequency and duration of cyclic compressive loading. Annals of Biomedical Engineering. 2001;29(6):476–82. pmid:11459341
View Article
PubMed/NCBI
Google Scholar

[35] View Article

[36] PubMed/NCBI

[37] Google Scholar

[ref11] 11. Wikswo JP. The relevance and potential roles of microphysiological systems in biology and medicine. Experimental Biology and Medicine. 2014;239(9):1061–72. pmid:25187571
View Article
PubMed/NCBI
Google Scholar

[39] View Article

[40] PubMed/NCBI

[41] Google Scholar

[ref12] 12. Rouwkema J, Gibbs S, Lutolf MP, Martin I, Vunjak-Novakovic G, Malda J. In vitro platforms for tissue engineering: Implications for basic research and clinical translation. Journal of Tissue Engineering and Regenerative Medicine. 2011;5(8):e164–e7. pmid:21774080
View Article
PubMed/NCBI
Google Scholar

[43] View Article

[44] PubMed/NCBI

[45] Google Scholar

[ref13] 13. Goldring MB, Goldring SR. Articular cartilage and subchondral bone in the pathogenesis of osteoarthritis. Annals of the New York Academy of Sciences2010. p. 230–7.
View Article
Google Scholar

[47] View Article

[48] Google Scholar

[ref14] 14. Suri S, Walsh DA. Osteochondral alterations in osteoarthritis. Bone. 2012;51(2):204–11. pmid:22023932
View Article
PubMed/NCBI
Google Scholar

[50] View Article

[51] PubMed/NCBI

[52] Google Scholar

[ref15] 15. Yuan XL, Meng HY, Wang YC, Peng J, Guo QY, Wang AY, et al. Bone-cartilage interface crosstalk in osteoarthritis: Potential pathways and future therapeutic strategies. Osteoarthritis and Cartilage. 2014;22(8):1077–89. pmid:24928319
View Article
PubMed/NCBI
Google Scholar

[54] View Article

[55] PubMed/NCBI

[56] Google Scholar

[ref16] 16. Martin I, Miot S, Barbero A, Jakob M, Wendt D. Osteochondral tissue engineering. Journal of biomechanics. 2007;40(4):750–65. Epub 2006/05/30. pmid:16730354.
View Article
PubMed/NCBI
Google Scholar

[58] View Article

[59] PubMed/NCBI

[60] Google Scholar

[ref17] 17. Giannoni P, Lazzarini E, Ceseracciu L, Barone AC, Quarto R, Scaglione S. Design and characterization of a tissue-engineered bilayer scaffold for osteochondral tissue repair. J Tissue Eng Regen Med. 2015;9(10):1182–92. Epub 2012/11/23. pmid:23172816.
View Article
PubMed/NCBI
Google Scholar

[62] View Article

[63] PubMed/NCBI

[64] Google Scholar

[ref18] 18. Tampieri A, Sandri M, Landi E, Pressato D, Francioli S, Quarto R, et al. Design of graded biomimetic osteochondral composite scaffolds. Biomaterials. 2008;29(26):3539–46. Epub 2008/06/10. pmid:18538387.
View Article
PubMed/NCBI
Google Scholar

[66] View Article

[67] PubMed/NCBI

[68] Google Scholar

[ref19] 19. Yang G, Rothrauff BB, Lin H, Gottardi R, Alexander PG, Tuan RS. Enhancement of tenogenic differentiation of human adipose stem cells by tendon-derived extracellular matrix. Biomaterials. 2013;34(37):9295–306. pmid:24044998
View Article
PubMed/NCBI
Google Scholar

[70] View Article

[71] PubMed/NCBI

[72] Google Scholar

[ref20] 20. Candiani G, Raimondi MT, Aurora R, Lagana K, Dubini G. Chondrocyte response to high regimens of cyclic hydrostatic pressure in 3-dimensional engineered constructs. The International journal of artificial organs. 2008;31(6):490–9. Epub 2008/07/09. pmid:18609501.
View Article
PubMed/NCBI
Google Scholar

[74] View Article

[75] PubMed/NCBI

[76] Google Scholar

[ref21] 21. Mazzei D, Guzzardi MA, Giusti S, Ahluwalia A. A low shear stress modular bioreactor for connected cell culture under high flow rates. Biotechnology and Bioengineering. 2010;106(1):127–37. pmid:20091740
View Article
PubMed/NCBI
Google Scholar

[78] View Article

[79] PubMed/NCBI

[80] Google Scholar

[ref22] 22. Sbrana T, Ucciferri N, Favrè M, Ahmed S, Collnot EM, Lehr CM, et al. Dual flow bioreactor with ultrathin microporous TEER sensing membrane for evaluation of nanoparticle toxicity. Sensors and Actuators, B: Chemical. 2016;223:440–6.
View Article
Google Scholar

[82] View Article

[83] Google Scholar

[ref23] 23. Vozzi F, Mazzei D, Vinci B, Vozzi G, Sbrana T, Ricotti L, et al. A flexible bioreactor system for constructing in vitro tissue and organ models. Biotechnology and Bioengineering. 2011;108(9):2129–40. pmid:21495015
View Article
PubMed/NCBI
Google Scholar

[85] View Article

[86] PubMed/NCBI

[87] Google Scholar

[ref24] 24. Lin H, Zhang D, Alexander PG, Yang G, Tan J, Cheng AWM, et al. Application of visible light-based projection stereolithography for live cell-scaffold fabrication with designed architecture. Biomaterials. 2013;34(2):331–9. pmid:23092861
View Article
PubMed/NCBI
Google Scholar

[89] View Article

[90] PubMed/NCBI

[91] Google Scholar

[ref25] 25. Mannella GA, Conoscenti G, Carfì Pavia F, La Carrubba V, Brucato V. Preparation of polymeric foams with a pore size gradient via Thermally Induced Phase Separation (TIPS). Materials Letters. 2015;160:31–3.
View Article
Google Scholar

[93] View Article

[94] Google Scholar

[ref26] 26. Taffetani M, Gottardi R, Gastaldi D, Raiteri R, Vena P. Poroelastic response of articular cartilage by nanoindentation creep tests at different characteristic lengths. Medical Engineering and Physics. 2014;36(7):850–8. pmid:24814573
View Article
PubMed/NCBI
Google Scholar

[96] View Article

[97] PubMed/NCBI

[98] Google Scholar

[ref27] 27. Nava MM, Raimondi MT, Pietrabissa R. Controlling self-renewal and differentiation of stem cells via mechanical cues. Journal of Biomedicine and Biotechnology. 2012;2012.
View Article
Google Scholar

[100] View Article

[101] Google Scholar

[ref28] 28. Sacco R, Causin P, Zunino P, Raimondi MT. A multiphysics/multiscale 2D numerical simulation of scaffold-based cartilage regeneration under interstitial perfusion in a bioreactor. Biomechanics and Modeling in Mechanobiology. 2011;10(4):577–89. pmid:20865436
View Article
PubMed/NCBI
Google Scholar

[103] View Article

[104] PubMed/NCBI

[105] Google Scholar

[ref29] 29. Michaelis L, Menten MML. The kinetics of invertin action: Translated by T.R.C. Boyde Submitted 4 February 1913. FEBS Letters. 2013;587(17):2712–20. pmid:23867202
View Article
PubMed/NCBI
Google Scholar

[107] View Article

[108] PubMed/NCBI

[109] Google Scholar

[ref30] 30. Astorino M, Hamers J, Shadden SC, Gerbeau JF. A robust and efficient valve model based on resistive immersed surfaces. International Journal for Numerical Methods in Biomedical Engineering. 2012;28(9):937–59. pmid:22941924
View Article
PubMed/NCBI
Google Scholar

[111] View Article

[112] PubMed/NCBI

[113] Google Scholar

[ref31] 31. Braakman R, Sipkema P, Westerhof N. A dynamic nonlinear lumped parameter model for skeletal muscle circulation. Annals of Biomedical Engineering. 1989;17(6):593–616. pmid:2589694
View Article
PubMed/NCBI
Google Scholar

[115] View Article

[116] PubMed/NCBI

[117] Google Scholar

[ref32] 32. Formaggia L, Lamponi D, Tuveri M, Veneziani A. Numerical modeling of 1D arterial networks coupled with a lumped parameters description of the heart. Computer Methods in Biomechanics and Biomedical Engineering. 2006;9(5):273–88. pmid:17132614
View Article
PubMed/NCBI
Google Scholar

[119] View Article

[120] PubMed/NCBI

[121] Google Scholar

[ref33] 33. Kim HJ, Vignon-Clementel IE, Figueroa CA, Ladisa JF, Jansen KE, Feinstein JA, et al. On coupling a lumped parameter heart model and a three-dimensional finite element aorta model. Annals of Biomedical Engineering. 2009;37(11):2153–69. pmid:19609676
View Article
PubMed/NCBI
Google Scholar

[123] View Article

[124] PubMed/NCBI

[125] Google Scholar

[ref34] 34. Lee J, Smith NP. The multi-scale modelling of coronary blood flow. Annals of Biomedical Engineering. 2012;40(11):2399–413. pmid:22565815
View Article
PubMed/NCBI
Google Scholar

[127] View Article

[128] PubMed/NCBI

[129] Google Scholar

[ref35] 35. Pant S, Fabrèges B, Gerbeau JF, Vignon-Clementel IE. A methodological paradigm for patient-specific multi-scale CFD simulations: From clinical measurements to parameter estimates for individual analysis. International Journal for Numerical Methods in Biomedical Engineering. 2014;30(12):1614–48. pmid:25345820
View Article
PubMed/NCBI
Google Scholar

[131] View Article

[132] PubMed/NCBI

[133] Google Scholar

[ref36] 36. Shi Y, Lawford P, Hose R. Review of Zero-D and 1-D Models of Blood Flow in the Cardiovascular System. BioMedical Engineering Online. 2011;10.
View Article
Google Scholar

[135] View Article

[136] Google Scholar

[ref37] 37. Raimondi MT, Giordano C, Pietrabissa R. Oxygen measurement in interstitially perfused cellularized constructs cultured in a miniaturized bioreactor. Journal of Applied Biomaterials and Functional Materials. 2015;13(4):e313–e9. pmid:26450634
View Article
PubMed/NCBI
Google Scholar

[138] View Article

[139] PubMed/NCBI

[140] Google Scholar

Figures

Abstract

1 Introduction

2 Models and Methods

2.1 CAD model

2.2 Flow

2.3 Mass transport

2.4 Model parameters and constitutive laws

2.4.1 Model parameters for the flow model.

2.4.2 Model Parameters and constitutive laws for mass transport.

2.5 Computational solvers

3 Lumped Parameter Models of HTB

3.1 Lumped parameter model for a fixed HTB configuration

3.2 Lumped parameter model for variable bioreactor configurations

4 Numerical Simulations

4.1 Numerical simulation of flow

4.2 Numerical simulation of transport

4.3 Oxygen consumption

4.4 Comparison of distributed and lumped parameter models

4.4.1 Fluid dynamics.

4.4.2 Mass transport.

5 Discussion

6 Conclusions and Perspectives

Author Contributions

References