Design , characterization , and modeling of microcirculation systems with integrated oxygenators

Here, we describe a microfluidic system for hypoxia assays on human cell culture models. These systems are developed to replace or reduce animal testing in biomedical basic research. The presented system uses a gas-permeable membrane as a gas–liquid interface and a micropump for media actuation to influence the oxygen content in two cell culture chambers. To apply well-defined hypoxic conditions to the cells, a good understanding of the mass transport phenomena is necessary. Therefore, a complete network model of the microfluidic system is presented. This model is validated by means of micro-particle image velocimetry (μPIV) and optical oxygen measurement with fluorescence lifetime detection. Finally, the impact of several process parameters, e.g., the gas permeability of the pump, is discussed using the developed model.


Introduction
In the human body, oxygen supply is a key factor influencing cell vitality and growth.An oxygen restriction is called hypoxia and induces several changes in cell metabolism.For example, the release of hypoxia inducible factors (HIF) during hypoxia influences inflammatory processes, e.g., affecting tumor growth (Eltzschig and Carmeliet, 2011) or wound healing (Kalucka et al., 2013).A deeper understanding of those processes is therefore of particular interest in medical basic research.Currently, two types of hypoxia assays are widely used: animal models with genetically modified mice (McColm et al., 2004) or static human cell culture models (Richter et al., 2013).Both have several advantages and disadvantages.Animal models, on one hand, feature a complete living system with several organs interacting with each other.Unfortunately, recent publications showed that, for example, substances tested as non-toxic for animals may be harmful for humans in the clinical trials (Knight, 2008).Human cell culture systems, on the other hand, give the opportunity to observe those processes in human tissue, but the interactions between several organs are not considered.Perfused microfluidic 3-D cell cultivation systems may overcome these issues and should mimic in vivo-like conditions much better than static cell cultivation platforms or animal models by integrating different human micro-organoids in a closed, artificial microvascular system (Baker, 2011).Based on a microfluidic platform for cell cultivation with integrated micropumps, valves, and reservoirs (Marx et al., 2012) a multilayer-based system was recently developed (Sonntag et al., 2015).The aim is to transfer pharmaceutical and cosmetic substance testing from the animal to this so-called multiorgan chip (MOC).Several human tissue were successfully cultivated within this system, e.g., liver (Maschmeyer et al., 2015), skin (Ataç et al., 2013), and kidney (Schimek et al., 2015).In the present study, an additional oxygenator element is implemented in the microfluidic platform, allowing automated and reproducible hypoxia assays on the chip.

Design
One of the main benefits of the presented microfluidic system is the possibility to integrate microfluidic actuators to produce completely closed and perfused microfluidic circuits and microvascular systems with very small volumes in the range of several microliters.External pumping devices are not needed, as it is mostly the case in microfluidic cell cultivation devices presented in the literature (Bhatia and Ingber, 2014).The microfluidic system is produced with a layer-by- layer manufacturing technology of laser-cut polymer foils (Sonntag et al., 2015).The mentioned fluidic actuators are pneumatically driven systems.Therefore, a flexible membrane has to be integrated in the microfluidic system, which can be easily displaced with pressure rates up to 100 kPa.A commercially available silicone foil (SILPURAN ® FILM 2030, Wacker Chemie AG) with a thickness of 200 microns and high gas permeability is perfectly suited for this application.In Fig. 1a an exploded view of the microfluidic system with the pneumatic part and mounted reservoirs on top, the elastomeric membrane in the middle and the fluidic part at the bottom is shown.
The complete system is mounted in a cell culture support with two connection elements which include the pneumatic ports.In Fig. 1b the developed microfluidic layout is shown.It features a micropump for fluid actuation and several valves to control the liquid stream.The system can be operated as a flow-through system by closing the bypass valve and opening the inlet and outlet valves.Furthermore, a closed circuit operation is possible by opening the bypass valve and closing the inlet and outlet valves, which limits the circulating volume to around 75 µL.Two cell cultivation chambers, A and B (with oxygen sensing spots), are integrated; each is equipped with an upstream and downstream valve (α, β) to control the perfusion to the segment.In previous fluidic designs, the oxygen entry was only performed by the gas-permeable elastomeric pump membrane.So the pump acted as a combined flow source and oxygenator (Busek et al., 2015a).This may be problematic especially when gas bubbles are produced, as reported by Godolsky and Knapp (2013).To enhance the oxygen transport capability an additional serpentine-shaped oxygenator element was placed at the center of the microfluidic system, utilizing the same gas-permeable membrane also used for fluid actuation but providing a much higher gas exchange area A O 2 .The process gas is flushed above the mem-brane in counterflow to the liquid stream with a flow of several L min −1 to prevent oxygen depletion over the length of the oxygenator.Finally, in Fig. 1c a picture of the bonded system with mounted reservoirs is shown.A controlling system based on an embedded Linux PC was developed switching up to 24 pneumatic outlets, which allows the parallel actuation of 8 pneumatic micropumps.Figure 2 shows a block diagram according to ISO-1912-1 (ISO, 2012).
Furthermore, the controller provides up to four electrical heating channels to ensure stable cell cultivation conditions.Several interfaces (USB, CAN, and Ethernet) allow connection to laboratory automation and information systems or direct sensor readout.As shown, an additional oxygenator output is provided which mixes two process gases, e.g., oxygen and nitrogen, to provide a controlled media oxygenation at integrated oxygenator elements.
It is therefore well suited to characterize the transient flow behavior in fully closed microfluidic systems.Small particles flowing with the media stream are utilized to observe the fluid movement.The experimental setup was described in detail earlier (Busek et al., 2012).For image processing an optimized cross-correlation algorithm is used, which only works when the flow velocity vector has only one component.This is, for example, the case in the straight channel segments at the outlet of the oxygenator device.Therefore, µPIV measurements were applied there.In Fig. 3 the velocity-time curve at the channel center for two pumping frequencies f (1.2 and 2.4 Hz) as well as their respective time-averaged values are shown.The pump was operated with a pressure of 50 kPa during pump cycles and a pressure of −50 kPa during filling cycles.
It can be seen, that the flow is pulsatile with peak velocities up to 400 mm s −1 and that there are time periods where no fluid movement can be observed at all.A detailed description and modeling of the micropump was published earlier (Busek et al., 2013).To obtain the time-averaged volumetric flow rate the measured velocity has to be integrated over the time and the channel cross section A: where h is the channel height, b the channel width, T the duration of one pump cycle (= f −1 ), v z the time-averaged flow velocity at the channel center, and v norm (y, z) the normalized flow pattern in both spatial directions.In a tubular channel a parabolic flow pattern can be estimated and the surface integral of v norm (y, z) is 0.5.In case of a rectangular channel (b = 1 mm, h = 0.25 mm) this factor changes to 0.48 because the flow pattern is getting flatter.The flow pattern of a rectangular channel can be calculated using Fourier series expansion (Sommer, 2014).With Eq. ( 1) the mean flow rate for the velocity-time curve shown in

Oxygen measurement with fluorescence lifetime detection
In presence of oxygen, the fluorescence lifetime of several fluorescent dyes is decreased due to quenching (Schmälzlin et al., 2005).This decay in the fluorescence lifetime can be utilized to measure the oxygen content in the chambers.The experimental setup was described earlier in detail (Schmieder et al., 2014).An oxygen sensing spot was placed in each chamber (A and B) and the fluorescence lifetime was measured for several minutes.For calibration purpose the microfluidic system was flushed with three different process gases with varying oxygen contents w O 2 compared to the volume fraction of oxygen in the atmosphere under standard conditions (φ O 2 = 20.9vol %).Those gases were air: w O 2 = 100 %; reference gas: w O 2 = 47 %; and nitrogen: w O 2 = 0 %.Afterwards, the microfluidic system was filled with deionized water for the oxygenation and deoxygenation experiments.During all experiments, the oxygenation and deoxygenation behavior of the system was characterized in the same manner.First pure nitrogen was applied to the membrane oxygenator, until the oxygen content at the measurement spots reached a lower boundary value.Then the process gas was changed to compressed air until the oxygen content remained constant.As mentioned previously the oxygenator uses the same membrane as the micropump.Therefore, an additional gas exchange at the membrane can be estimated as described earlier (Busek et al., 2015b).The oxygenation-deoxygenation behavior is thus not only influenced by the process gas applied to the oxygenator but also by the oxygen content of the pumping gas.This effect is shown in Fig. 4. One can see the measured oxygen content plotted against pumping time for different pump gases (compressed air, nitrogen).The mass transfer takes place in several minutes.Furthermore, the used pump gas has a significant influence on the oxygenation-deoxygenation behavior of the system.

Modeling
The mass transport under pulsatile flow can be described with the dimensionless Reynolds number Re and the Strouhal number St (Nishimura et al., 2000).For the described rectan-gular channel both numbers can be calculated using following formula: With a kinematic viscosity ϑ of 1 mm 2 s −1 , a frequency f of 1 Hz, and the channel dimensions described in Sect.3.1, Re is 6 and St is 0.08.Both values are much lower than the critical values (Re crit = 1400 and St n = 1.92)where turbulence, especially in cavities, can occur, which will enhance mass transport there (Nishimura et al., 2000).Therefore, the flow is laminar and can be seen as quasistatic.One can replace the pulsatile micropump by a constant flow source in the mass transport model.

Oxygenation coefficient
By keeping the process gas flow at the top of the membrane in the range of several L min −1 oxygen depletion over the length of the oxygenator does not take place because the liquid flow is with several µL s −1 much lower.Therefore, the oxygen source concentration is constant and the oxygenation coefficient K O 2 can be used to characterize the oxygenator K O 2 is a function of the oxygen concentration at the oxygenator outlet c(L) and the saturation concentration c s, O 2 which depends on the solubility α L of oxygen in the liquid and can be calculated as follows: A high oxygenation coefficient is useful to allow a fast and broad oxygenation and deoxygenation of the culture media but this should be realizable with the presented manufacturing technology and with pumping rates of several µL s −1 .Therefore, a mathematic model describing the membranebased oxygenator should be found.Mass transport phenomena in perfused systems are generally described by the convection-diffusion equation: where D describes the diffusion coefficient, v is the velocity vector, and R is the reactive term.Reaction means here the oxygen uptake of the cells.Because of the fact that there are no cells present in the oxygenator, R can be neglected.Furthermore, only the steady state should be considered, meaning that δc δt can be set to zero.Equation ( 5) can be further simplified when we have a look at the geometrical boundary conditions of the oxygenator as shown in Fig. 5.With a constant oxygen flux J D through the membrane, the concentration gradient in x direction δc δx is zero.Therefore, a convective mass transport occurs in z direction and diffusion mainly occurs in y direction.Now the partial differential equation (PDE) (Eq.5) only depends on the two Cartesian coordinates y and z.For a fully developed laminar flow one can formulate the following PDE: Equation ( 6) is similar to a first-order wall reaction (Pallares and Grau, 2014).If the liquid flow is simplified to a plug flow with the mean velocity v z , the velocity does not depend on the y coordinate and the error function is a practicable solution of this PDE (Bird et al., 2007): The mean oxygen concentration at the oxygenator outlet (c m (z = L)) can be obtained by integrating this equation over the channel height h.Furthermore, a numerical fit for the error function with the mass transport term f (y, z) and the coefficients a 1 ≈ 0.27839, a 2 ≈ 0.23039, a 3 ≈ 0.00097, a 4 ≈ 0.07810 (Abramowitz and Stegun, 1970) can be used: Now the oxygenation coefficient can be calculated as follows: for an oxygenator with a length L = 110 mm, a channel height h = 0.25 mm, a diffusivity D ≈ 2×10 −3 mm 2 s −1 , and a mean velocity v z ≈ 30 mm s −1 the oxygenation coefficient is K O 2 ≈ 0.38, in alignment with hollow-fiber oxygenators widely described in the literature (Kreulen et al., 1993;Dindore et al., 2005).

Transient oxygen transport through the membrane
Besides the perfused channel area, the oxygenator membrane has a significant influence on the oxygenation and deoxygenation behavior of the system.When the process gas at the top of the membrane is switched, diffusion processes in the membrane will delay the oxygen exchange at the membraneliquid interface.Starting from Eq. ( 5), following PDE can be assumed describing the mass transport in the membrane (no convection): with the membrane thickness d and the overall oxygen transport coefficient k.This factor is composed of the gas-liquid mass transfer rate k G, L and the membrane mass transfer coefficient k M and can be calculated as follows (Zhan et al., 2009): where D denotes the oxygen diffusivity in the silicone.For a 200 µm thick membrane and a diffusivity of D ≈ 3 × 10 −3 mm 2 s −1 and a k i, L ≈ 0.04 mm s −1 (Shiku et al., 2006) k is about 0.01 mm s −1 .A solution for Eq. ( 11) at the membrane bottom (y = d) is analogous to Eq. ( 7) given with the error function In case of deoxygenation, the source concentration c s, O 2 falls to zero and following formula can be assumed: www.j-sens-sens-syst.net/5/221/2016/ 9 Beside the described membrane oxygenator (Eq. 10, Eq. 13, Eq. 14) it includes a constant flow source (Eq. 1) with an additional oxygenator element mimicking the gas exchange at the pump membranes.The mathematical description of the additional oxygen entry is difficult, because the flow pattern in the micro pump is unknown.Therefore the oxygenation coefficient of the micro pump K O2, pump was obtained from experimental results.Furthermore fluidic resistances representing the micro channels and additional volume elements are implemented.The flowing liquid is water with the shown parameters 5 (diffusivity D, dynamic viscosity ν and saturation concentration cs).This model can now be used to compare simulated and measured flow and oxygen values and afterwards an optimization of the fluidic layout can be obtained e.g. by coupling SimulationX to the optimization tool OptiY (Abid et al., 2015).

Model validation
In Fig. 7 (left), the oxygen content w O2 (relatively to the atmospheric oxygen volume fraction 20.9 . %) at spot B 10 is shown.The plot shows the calculated values (with the network model and the parameters shown in Fig. 6) and the measured w O2 values.Similar to the experimental setup described earlier the oxygenator was initially operated with pure nitrogen for deoxygenation and after around 300 seconds the process gas was changed to air.The modelled results fit the experimental data with a root mean square deviation of 0.86 %.Therefore the presented model is suitable to simulate the oxygen transport by the oxygenator and the micro pump within the system.15

Impact of pump oxygenation coefficient
The impact of different process parameters on the oxygenation and deoxygenation curves is of great interest, because it 20 shows the best way to manipulate the oxygen content in the cell culture chambers.As mentioned earlier, the parasitic oxygen entry by the pump is unknown.It is therefore important to examine the influence of this parameter on the oxygenationdeoxygenation-behaviour of the system.In Fig. 7 (right) this effect is shown for four different pump oxygenation coefficients K O2, pump calculated with the developed network model of the microfluidic system.

Network model of the complete microfluidic system
All presented parts of the modeling are included in the microfluidic library for the network simulation tool Simula-tionX (Busek et al., 2015a) using the model description language Modelica (Fritzson and Engelson, 1998).A complete model of the microfluidic system can be developed as shown in Fig. 6.Besides the described membrane oxygenator (Eqs.10, 13, 14), it includes a constant flow source (Eq. 1) with an additional oxygenator element mimicking the gas exchange at the pump membranes.The mathematical description of the additional oxygen entry is difficult, because the flow pattern in the micropump is unknown.Therefore, the oxygenation coefficient of the micropump K O 2 , pump was obtained from experimental results.Furthermore, fluidic resistances representing the microchannels and additional volume elements are implemented.The flowing liquid is water with the shown parameters (diffusivity D, dynamic viscosity ν, and saturation concentration cs).This model can now be used to compare simulated and measured flow and oxygen values and afterwards an optimization of the fluidic layout can be obtained, e.g., by coupling SimulationX to the optimization tool OptiY (Abid et al., 2015).

Model validation
In Fig. 7 (left), the oxygen content w O 2 (relative to the atmospheric oxygen volume fraction φ O 2 = 20.9vol %) at spot B is shown.The plot shows the calculated values (with the network model and the parameters shown in Fig. 6) and the measured w O 2 values.Similar to the experimental setup described earlier, the oxygenator was initially operated with pure nitrogen for deoxygenation and after around 300 s the process gas was changed to air.The modeled results fit the experimental data with a root mean square deviation of 0.86 %.Therefore, the presented model is suitable to simulate the oxygen transport by the oxygenator and the micropump within the system.

Impact of pump oxygenation coefficient
The impact of different process parameters on the oxygenation and deoxygenation curves is of great interest, because it shows the best way to manipulate the oxygen content in the cell culture chambers.As mentioned earlier, the parasitic oxygen entry by the pump is unknown.It is therefore important to examine the influence of this parameter on the oxygenation-deoxygenation behavior of the system.In Fig. 7 (right) this effect is shown for four different pump oxygenation coefficients K O 2 , pump calculated with the developed network model of the microfluidic system.

Conclusion and outlook
The aim of this work was to develop a microfluidic system for hypoxia assays in human cell cultures to replace or reduce animal testing in basic medical research.The presented system uses a gas-permeable membrane as a gas-liquid interface and a micropump for media actuation to influence the oxygen content in two cell culture chambers.To apply welldefined hypoxic conditions to the cells, a good understanding of the mass transport phenomena is necessary which is why a complete model of the microfluidic system was developed.The flow in the microchannels was measured with a µPIV setup and used as constant flow source in the model.Afterwards the oxygenation and deoxygenation behavior of the system was characterized by applying different process gases to the oxygenator and observing the fluorescence decay time in both measurement chambers.A good correlation between measured and calculated oxygen concentrations was achieved, which means that the developed model is suited to calculate the oxygen transport in the microfluidic system.Finally, the impact of the parasitic oxygen entry by the pump was followed up using the developed model.In further investigations the influence of other process parameters, e.g., the process gas pressure should be examined.Moreover, the model can be coupled to an optimization tool to refine the design of the microfluidic system in terms of needed oxygenation coefficients, etc.Later on an oxygen consumer, e.g., a high-density 3-D cell culture model, has to be implemented in the microfluidic system.A mathematical model with different reaction kinetics is already integrated in the simulation tool.Together with the here-presented oxygenator model it could be used to calculate the steady-state oxygen concentration for different oxygen consumption rates or even cell growth as presented in the literature (Ma et al., 2007).

Figure 1 .
Figure 1.(a) Exploded view of the microfluidic system; (b) microfluidic layout indicating micropump, oxygenator, and measurement spots (A, B) as well as upstream and downstream valves (α, β); (c) picture of microfluidic system with mounted reservoirs.

Figure 6 .
Figure 6.SimulationX model of oxygen transport in the microfluidic system.

Fig. 7 :
Fig. 7: left: Comparison between modelled and measured oxygen content for deoxygenation with nitrogen and oxygenation with compressed air; right: Impact of pump oxygenation coefficient on oxygenation-deoxygenation behaviour (calculated)

Figure 7 .
Figure 7. Left: comparison between modeled and measured oxygen content for deoxygenation with nitrogen and oxygenation with compressed air; right: impact of pump oxygenation coefficient on oxygenation-deoxygenation behavior (calculated).