Received date: June 15, 2013; Accepted date: July 17, 2013; Published date: July 22, 2013
Citation: Tanaka S, Iber D (2013) Simulation of a Microfluidic Gradient Generator using Lattice Boltzmann Methods. Curr Synthetic Sys Biol 1:102. doi: 10.4172/2332-0737.1000102
Copyright: © 2013 Tanaka S, 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.
Visit for more related articles at Current Synthetic and Systems Biology
Microfluidics provides a powerful and versatile technology to accurately control spatial and temporal conditions for cell culturing and can therefore be used to study cellular responses to gradients. Here we use Lattice Boltzmann methods (LBM) to solve both the Navier-Stokes equation (NSE) for the fluid and the coupled convection-diffusion equation (CDE) for the compounds that form the diffusion-based gradient. The design of a microfluidic chamber for diffusion-based gradients must avoid flow through the cell chamber. This can be achieved by alternately opening the source and the sink channels. The fast toggling of microfluidic valves requires switching between different boundary conditions. We demonstrate that the LBM is a powerful method for handling complex geometries, high Péclet number conditions, discontinuities in the boundary conditions, and multiphysics coupling.
The spatial organization of complex organisms requires cells to read out concentration differences and respond accordingly. Graded responses have been documented in different contexts, including developmental processes and immunological responses. To study these processes in greater detail it is important to study cells in a well-controlled setting where concentration gradients can be applied. Microfluidics provides a powerful and versatile technology to accurately control spatial and temporal conditions. A standard layout of a microfluidic set-up is shown in Figure 1. The chamber is connected to two channels that act as source and sink of the compound of interest. The compound is transported through these channels mainly by advection. Within the chamber, however, flow must be avoided as the resulting hydrodynamic stress would impact on the cells. Transport of the compound in the chamber must therefore be diffusion-dominated. The transition from advection- to diffusion-dominated regime at the right place is a challenging engineering problem. In the design phase, the use of adequate simulation techniques can provide insight into the governing mechanisms.
Figure 1: Simulation Setup. (a) The channel and chamber geometry (all numbers given in [μm]). The cell chamber of interest, in-between the two channels, has dimensions 100[μm] by 250[μm]. The dead end valves are always closed. (b) In reality the distance from the valves to the chamber would be longer, but the channels are pruned in the simulation setup in order to save simulation time. This has no impact on the solution since the diffusion length is smaller than the length of the pruned channels. (c) The left and right channel are opened alternately for 1[s], respectively. Between the flushings, all valves are closed for 119[s].
To describe the process, we need to simulate both the fluid dynamics and the distribution of the compound in the fluid. An incompressible Newtonian fluid is described by the Navier-Stokes equation (NSE):
where denotes the velocity field, ρ the fluid density, µ the dynamic viscosity, and an external body force, which is zero in our case. The passive advection and diffusion of a diluted compound is described by the advection-diffusion equation:
with the diffusion coefficient D, and the local reaction term R.
Several numerical schemes have been used to simulate fluid flow, including the Lattice Boltzmann method (LBM) . It has been shown that the Navier-Stokes equation is recovered in the nearly incompressible hydrodynamic limit. The LBM was also successfully applied to solve diffusion  and advection-diffusion equations . Furthermore, the LBM has been shown to be capable to simulate reaction-diffusion problems by applying it to Turing type mechanisms [4–6]. The coupling of a passively advected and diffusing scalar to a fluid was shown in [7,8].
The Lattice Boltzmann method has previously been used to simulate a non-flow-free microfluidic gradient generator , whose design was inspired by Dertinger et al. . However, the Pe´clet number in these simulations is about fivefold lower than in the experimental data that is used to validate the simulations and well within the numerically noncritical regime of the standard Lattice Boltzmann methods. Moreover, since the microfluidic gradient generator is not flow-free it would not be suitable for experiments with cells that are sensitive to flow.
Recently, a microfluidic setup was published that prevents flow through the cell chamber by using alternate flushing of the channels . We use LBM-based simulations to computationally explore a similar design. Because of its algorithmic locality, complex geometries and boundary conditions (such as the fast switching of boundary conditions to imitate the valve dynamics), can easily be integrated into the LBM. In the following we will describe the LBM-based simulations of gradient formation in a microfluidic chamber.
Experimental layout and conditions
The implemented layout of the microfluidic chip is shown in Figure 1a. It consists of two parallel channels of width 100[µm] and length 1000[µm], whose inlet and outlet are controlled by microfluidic valves. The two channels are connected by the diffusion chamber, consisting of a pre-chamber region, and the actual cell chamber of size 1000[µm] by 250[µm]. A fence consisting of five cylindrical elements disturbs the flow and thus reduces its advective influence in the chamber. Additional dead-end channels, which are used for cell handling, are explicitly modeled since they act as reservoirs. The height of the chamber is between 40[µm] and 100[µm], such that the processes are well approximated by a two dimensional simulation.
The distance from the valves to the chambers may be larger in reality (Figure 1b); however, they can be pruned such that they are in the range of the effective diffusion length , defined by the diffusion coefficient D and the time period T between two consecutive flushings.
Every two minutes, the left channel opens both its inlet and outlet valve for one second. As shown in Figure 1c, the sequence for the right channel is shifted by one minute, such that the chip is completely closed for 59 seconds. The technique to avoid flow through the cell chamber using alternate flushing of the channels has been described in .
In the beginning, the left channel is flushed with medium with a compound of interest, and the right channel with plain medium. Therefore a concentration gradient emerges across the chamber. After two hours, the fluids are interchanged for another two hours, enforcing an inversion of the gradient.
Upon opening the valve, the flow is virtually immediately fully developed due to the very low Reynolds number and the incompressibility. The high maximal flow speed in the channel Uphys=5000[µm/s] and the low diffusion coefficient of the compound Dphys= 122[µm/s] result in an advection-dominated transport in the channels. Using a characteristic length scale L = 100[µm/s] (the channel width), the Pe´clet number Pe = UL/D≈4100 can be computed. While entering the diffusion chamber (of after closing all valves), the advective transport reduces to zero and the diffusive transport dominates.
A summary of the used physical parameter values can be found in Tables 1 and 2.
|Physical units||LB units|
Table 1: Parameter values for the fluid solver, both in physical and LB units
|Physical units||LB units|
Table 2: Parameter values for the CDE solver, both in physical and LB units.
For the fluid, we employed the standard D2Q9 solver proposed in Chen and Doolen . The D2Q9 lattice is defined by the lattice directions
The lattice Bhatnagar-Gross-Krook (BGK) equation, the equilibrium distribution function in the ith direction, and the corresponding weights wi are taken as:
where Δt denotes the time step, Δx the spatial discretization and c= Δx/Δt the lattice speed. The simplest choice to guarantee consistency with the lattice (Eq. (3)) is Δt=Δt =c=1. The relaxation time τ is related to the kinematic viscosity . For isothermal flow, the speed of sound is defined as
Equation (4) implies a two step algorithm: first, local collision relaxes the populations towards the local equilibrium (right hand side), and second, the populations perform a free flight to the next lattice point (left hand side). The density, momentum density and pressure are computed as:
For solving the advection-diffusion equation of a compound C (Equation (2)), a multi-distribution function (MDF) approach was used. It was first used in Bartoloni et al.  to model the temperature field of a thermal flow (Boussinesq approximation) as a passively advected scalar field, and further improved in Guo et al. . In contrary to the fluid solver,  use a D2Q4 lattice and stencil. The lattice BGK equation, the equilibrium distribution and the zeroth order moment read:
where gi denote the particle distribution functions. The velocity field is transferred from the fluid solver (Equation (7)). The relaxation time time τg is related to the diffusion coefficient as D = (2τg -1)/4.
Note that all variables are measured in tthe LB units δt and δx. The conversion to physical quantities is described below.
For the wall boundary condition of the fluid, the missing incoming populations are approximated by equilibrium distributions. The momentum needed to compute the equilibrium is spatially first order interpolated between the fluid in direction of the missing population, and the known zero momentum at the wall. The density is spatially first order extrapolated from the fluid.
The pressure boundary condition for the fluid is realized by applying do-nothing, corresponding to zeroth order extrapolation in time. Then all populations are rescaled such that the desired pressure is obtained. The velocity field is not affected, since the rescaling factor γ cancels
This technique was used in Zhang and Kwok , although for specifying the pressure difference in periodically closed channels.
The same approach is used for the Dirichlet boundary condition of the diffusing compound at the channel inlet. Since the prescribeddensity boundary condition cannot be applied to the outlet, the wellknown do-nothing boundary condition is applied . For the no-flux boundary condition, the boundary condition presented in Guo et al. [8,15] is used.
Choice of parameter values and conversion
For both the fluid and the advection-diffusion simulations, the flow channels are resolved by 20 lattice points. The lattice spacing thus derives as δx =[100[µm]/20=5[µm].
For the fluid simulation, the dimensionless Reynolds-number can be used to determine the missing LB parameters:
With the choice of a numerically reasonable LB kinematic viscosity vLB = 0.05, the maximal fluid speed ULB can be computed (Table 1). The LB time unit can be computed by equating the physical and LB maximal fluid speed, which results in δt = 1.25 ×10-6[s]. The pressure difference is chosen such that ULB is achieved with less than 1% error.
The procedure for the advection-diffusion simulation is similar. In order to allow for large time steps, the maximal fluid speed is chosen as ULB = 0.5, which is close to the speed of sound . Again, the LB time unit δt = 1.25 ×10-4[s] is determined by comparing the maximal fluid speeds. To use the fluid field from the fluid simulation in the advection-diffusion simulation, the field was scaled according to the change in the LB time unit. By equating the dimensionless Pe´clet number
the the diffusion coefficient DLB can be computed (Table 2).
The diffusing compound is passively advected by the fluid, but does not feedback on the fluid dynamics. The fluid dynamics can therefore be solved independently. After having computed the velocity fields, they can be used as an input to solve the advection-diffusion equation of the compound. The results of these steps are analyzed and discussed in the following.
In a first step, only the fluid dynamics is solved and analyzed. The geometry as shown in Figure 1a was implemented using the wellknown D2Q9 BGK LB scheme. Upon opening the valves, a pressure boundary condition with prescribed pressure was applied. The fluid in the domain, initially at rest, is brought to a steady state flow field in a short time span in a weakly compressible manner. The streamlines for opened left and right channel are shown in Figure 2a and 2b, respectively. Although the majority of streamlines follows a straight line, particles close to the boundaries flow around the fence and into the chamber. It can be observed that the direction of flow is locally reversed. However, the magnitude of the velocity is negligible in the chamber itself (Figure 2c and 2d), which is an important requirement for culturing cells. In the channel, the flow approximates a Poiseuille flow velocity profile, and corruptions as a consequence of the inlet and outlet boundary conditions are negligible. The fluid pressure, shown in Figure 2e and 2f, drops, according to the Poiseuille solution, linearly in the channels. Close to the chamber, the effective flow cross-section is increased, leading to a lower pressure gradient. The pressure field in the cell chamber is homogeneous. The chosen parameter values for the fluid dynamics simulation are given in Table 1.
Figure 2: The flow field. (a-b) Streamlines when the left (a) or right (b) channel are open, respectively. (c-d) The magnitude of velocity. In the channel, a Poiseuille-flow-like velocity profile is developed. The units of the color code is [μm/s]. (e-f) The pressure field drops linearly in the channels, and flattens in the chamber-transition-zone, where the effective flow-crosssection is increased. The color code is given in LB units.
Since the flow virtually immediately reaches state (after about 0.025[s]), and because the time step for the fluid solver has to be chosen considerably smaller (about tenfold) than the desired time step for the advection-diffusion solver, the velocity fields are pre-computed (as described in the a foregoing section), stored and then loaded into the advection-diffusion solver. The velocity fields of three different conditions are needed: left channel open and right channel closed, left channel closed and right channel open, and all channels closed. This approach prevents the advection-diffusion solver from solving the same fluid flow repeatedly.
In a second step, the advection-diffusion solver, which is largely consistent with the fluid solver, is implemented to solve (Equation (2)). Other than in the fluid solver, the velocity to compute the equilibrium distributions (Equation (9)) is taken from the pre-computed fluid velocity field. The zeroth order moment (density) can be interpreted as the compound concentration. In order to allow for maximal temporal step size, the fluid velocity has to be chosen as high as possible. We chose ULB = 0.5 in order to comply with the limit defined by the speed of sound cs. It has been shown that, for advection-diffusion applications, the advection velocity can be chosen comparably high [16,17], although the error grows. However, the choice of a new ULB = 0.5 (and thus a new time step δt) requires the rescaling of the pre-computed velocity fields.
In order to simulate the opening and closing of the valves, the inlet and outlet boundary conditions are changed periodically for the advection-diffusion solver. For closed valves, an ordinary wall (no-flux) boundary condition is applied. An open inlet is realized by a prescribed-density boundary condition, which corresponds to the pressure boundary condition in the fluid solver. The inlet density is set to 1 or 0, depending on whether the channel is fed with the concentrated or the plain medium.
When opening the channel carrying the concentrated medium for the first time, the advected concentration profile has a step-like shape. The advection of discontinuities (or extremely sharp gradients) is numerically challenging. To test the capability of the advectiondiffusion solver to advect a step profile with very low diffusion, the initial opening of the left channel is considered. At time t = 0, the channel is at rest, and the compound concentration is zero. Upon opening the channel, the fluid velocity immediately accelerates to the stationary (pre-computed) velocity field. The compound concentration, being 1 at the inlet boundary, is carried with the flow into the channel. This situation is similar to the advection of a step profile, with two differences: firstly, the velocity is not constant perpendicular to the direction of flow; and, secondly, we have diffusive transport. However, the Péclet number Pe ≈ 4100 is very high, meaning that the advective transport dominates the diffusive transport, and the smoothing of the initially step-like concentration profile is minor. Figure 3a shows the concentration profile along the channel. Severe numerical instabilities are induced close to the high gradient region, which do not decay with the front being advected. A tenfold decrease of the time step δt does not affect the instabilities (Figure 3a). Decreasing, instead, the Péclet number to Pe ≈ 800 by decreasing the flow speed to Uphys = 100[µm/s] (and thus ULB = 0.1[δx/δt]) leads to significant smoothing of the initially present numerical instabilities, as shown in Figure 3c. Again, decreasing the time step tenfold does not lead to an improvement of the solution (Figure 3d). However, since the numerical instabilities only form at the very first opening on a massive scale, and only on a minor scale at later openings, they are acceptable when considering the long time solution. The chosen parameters for the step advection test are summarized in Table 3.
Figure 3: Propagation of a step profile. The concentration profiles along the midline of the left channel is shown, shortly after opening the valve the first time. The top row shows simulations for a Péclet number Pe ≈ 4100, and the bottom row for Pe ≈ 800 after decreasing the fluid speed to Uphys = 1000[μm/s]. In the second column, the temporal resolution is tenfold increased as compared to the first columns. The chosen parameters are summarized in Table 3a For low temporal resolution and high Pe number, high numerical instabilities are triggered, which do not decrease as time passes. (b) Increased temporal resolution does not mitigate the instabilities. (c) At a five-fold lower Pe number, numerical instabilities can be initially found close to the high gradient, but they diminish as the step profile propagates and smoothens. (d) Higher temporal resolution does not lead to significantly different results.
Table 3: Parameter values for the step profile advection test.
As last step, the full microfluidic setup is simulated. The left channel is fed with concentrated medium for two hours, leading to a gradient from the left to the right. For another two hours, the gradient is reversed by changing the medium of the channels. In order to test the effect of decay of the compound, the simulations are once carried out with, and once without degradation. In the latter case, the gradient reaches almost its steady-state after 60 [min] (shown in Figure 4a). According to the theory, a non-linear steady-state concentration profile is obtained. Thirty minutes after reversing the media, the gradient already reaches its mirrored steady-state profile. However, the deadend channels act as reservoirs, revealing the importance of their inclusion into the simulation. When setting the compound degradation to zero, it takes approximately two hours to reach the linear steadystate concentration profile (shown in Figure 4b). The reversion is much faster and the steady-state is already reached after one hour. This can be explained by the fact that a particle only has to diffuse through half the chamber, whereas it had to diffuse through the entire chamber to form the initial gradient. To analyze the dynamics of gradient formation, the time evolution of the concentration at fixed points (shown in Figure 1 as blue, green and red points) is shown in Figure 4c for a decaying compound, and in Figure 4d for a non-decaying compound. Besides the time to reach steady-state, also the reached maximal concentrations in the cell chamber (from -500[µm] to + 500[µm]) are considerably lower in the case when the compound decays.
Figure 4: Evolution of Concentration. (a) The concentration profiles at different time points along the midline of the chamber. The compound is decaying, leading to a non-linear steady-state profile. (b) For a stable compound, a linear concentration profile is obtained in the chamber. (c) Time evolution of the concentration at three different locations in the chamber (Figure 1). The steady state is reached quickly for a decaying compound, as opposed to a stable compound (d).
In order to gain insight into the governing processes in a microfluidic gradient generator, and to support its design and development, computer simulations based on Lattice Boltzmann methods were developed. We have shown that the LBM is capable to simulate both the fluid flow and the coupled advective and diffusive processes. Although the high Péclet number regime (Pe > 4000) leads to severe numerical instabilities when advecting steep gradients, we report that, on the long time scale, the LBM leads to stable solutions. Using simulation it was shown that the presented microfluidic chip design is capable of forming and maintaining spatially and temporally controlled, diffusion based gradients. The valve-switching strategy and the flow-disturbing fence reliably prevent the flow from entering the cell culture chamber.
The benefit of the availability of such a gradient generator is manifold. Firstly, the concentration levels at points of interest can be accurately predicted and controlled. With the steadily increasing demand for quantitative data, microfluidic gradient generators offer the possibility to efficiently conduct multiple experiments on one chip. Secondly, the compound sources and sinks are spatially close to the experimental chamber and frequently renewed. This allows to use biochemical compounds with very high (self-) degradation, since the diffusive transport only has to cover comparably small distances.
The authors acknowledge funding from the SNF Sinergia grant Developmental engineering of endochondral ossification from mesenchymal stem cells and from the SNF SystemsX RTD NeurostemX.