Computational analysis of integrated biosensing and shear flow in a microfluidic vascular model

The forces exerted by flowing blood on vascular cells in vivo play critical roles in regulating vascular cell biology.¹ Traditional static in vitro cell culture models fail to replicate fluid flow-induced shear stresses, motivating the development of microfluidic platforms that better mimic the vascular microenvironment.^2,3 Indeed, microfluidic platforms have been extensively applied to study the effects of flow-induced shear stress on various cell types.⁴ 

Among the many cell behaviours affected by flow-induced shear stress is the secretion of biomolecules that impact cellular signaling and function. For example, nitric oxide, an important mediator of vasodilation, is released by endothelial cells in a shear-dependent manner.⁵ However, in situ microfluidic biosensing of cell-secreted biomolecules is typically performed under static fluid conditions.⁶ Assay integration is an oft-cited advantage of microfluidic cell culture devices,^7–9 but in situ detection of secreted biomolecules under flow conditions presents significant design challenges due to physiological shear stress constraints, and the coupled nature of fluid flow and mass transport to the biosensor surface. Moreover, flow significantly dilutes cell-secreted factor concentrations, hindering or even preventing detection. Predicting mass transport of secreted factors and the biosensor response is therefore critical to inform the design of microfluidic cell culture systems integrating biosensing capabilities.

Computational modeling is a useful tool for studying physical phenomena, assessing feasibility and guiding design choices in microfluidics with integrated quantitation assays.^10,11 The behaviour of surface-based biosensors has been well-characterized using theoretical and computational modeling, which provide excellent insight into convection-diffusion-reaction design considerations.^12–14 This analysis, however, is not sufficient to address the constraints of incorporating physiological levels of shear stress in the device.

Computational modeling of biosensing considerations is also generally limited to the scenario of a planar biosensor embedded within a single flow channel. This single layer design requires that the biosensor be directly adjacent to the cell layer, posing practical challenges to isolating cell culture and biosensor regions within the device. An alternative configuration is a bilayer device in which cells are cultured in the upper channel on a porous membrane support (e.g., Boyden chamber) while a surface-based biosensor is located in the lower channel underneath the membrane (Fig. 1A); the model, however, is generally applicable to other common multilayer configurations such as a hydrogel layer between two channels.^15,16 The bilayer membrane configuration is widely used in microfluidic cell culture applications,^17–20 and is particularly well-suited for microfluidic vascular models, where it has been used to study cancer cell²¹ and monocyte²² adhesion to an endothelial monolayer, endothelial permeability,^23,24 the vascular/valvular microenvironment²⁵ and the blood-brain barrier.^26,27 Bilayer culture devices have also been coupled with mass spectrometry for on-line monitoring of drug permeability.²⁸ For integrated biosensing applications, the bilayer configuration is advantageous as it positions the cells and biosensor closely together while isolating the biosensor from the elevated flow rates in the upper channel, allowing cell-secreted factors to accumulate in the lower channel.

Here we present a computational model of fluid flow, mass transport and biosensor kinetics to develop design guidelines that enable both on-chip detection of cell-secreted factors and application of physiological flow-induced shear stress to cultured cells in a bilayer microfluidic device. These guidelines detail the selection of critical design parameters of the bilayer device that will ensure optimal device functionality, which is contingent on meeting several important criteria. First, shear stress in the upper channel must be within physiological limits. Secondly, secreted biomolecule transport rates across the membrane must be maximized in order to efficiently supply analyte to the biosensor. Finally, the biosensor signal itself must also equilibrate rapidly in order to perform measurements within a practical timescale. Meeting these criteria enables the rapid and sensitive detection of secreted biomolecules from cells under a wide range of designer-specified shear stresses, adding novel and valuable functionality to microfluidic vascular models, further fulfilling their potential as lab-on-a-chip systems.

A schematic of the computationally-modeled bilayer microfluidic device is shown in Fig. 1B. Based on a model previously developed by Young et al.,²³ the cell monolayer and porous membrane support were modeled as a uniform porous medium. The microchannels, porous membrane and cell monolayer were layered in the y-direction. A small height-to-width aspect ratio was imposed on the upper and lower microchannel geometry and therefore, coupled with the low Reynolds numbers in the device, no variance of the model in the z-direction was assumed, resulting in a two dimensional model in the xy-plane. For simplicity, the biosensor geometry was modeled as a rectangular strip on the bottom surface of the lower channel. COMSOL Multiphysics 4.2 finite element analysis software was used to model and numerically simulate fluid flow, mass transport and reaction kinetics in the microfluidic device. A more detailed description of the implementation may be found in the ESI.

The governing equations and boundary conditions used to computationally model the bilayer microfluidic device are described below. A complete list of constant and simulation parameter values may be found in Table I. In order to achieve the desired device functionality, namely the integration of both biosensing and shear stress stimulation in a single device, proper selection of critical design parameter values is necessary. Table II contains a list of the design parameters considered in the guidelines developed in this paper.

Flow in the device was divided into three regions: laminar flow in the upper channel generated by an external source (e.g., pump); porous media flow through the uniform porous medium representing the cell monolayer and porous membrane; and laminar flow in the lower channel generated by a positive pressure difference between the upper and lower channel, and where no external flow was applied. The fluid, i.e., cell culture media, was assumed to have the density and dynamic viscosity of water.

Free flow in the upper and lower channels is governed by the Navier-Stokes equations for incompressible flow:

where the density ρ (kg m^-3) and dynamic viscosity μ (Pa·s) of the fluid were that of water, u is the velocity field of the fluid within the microchannels and p is the pressure. The upper channel flow boundary conditions consisted of specifying pressure values at the inlet and outlet boundaries. The lower channel backpressure was assumed to be greater than the upper channel backpressure, leading to a hydrostatic pressure head P_head between the upper and lower channel outlets. The lower channel flow boundary condition consisted of specifying P_head at the two lower channel outlet boundaries (i.e., both the left and right side boundaries are outlets), which was selected to be equivalent to 98.1 Pa, the backpressure from an upper and lower channel outlet height difference of 10 mm. A no-slip condition was assumed at all microchannel walls and the biosensor surface.

Porous media flow in the uniform porous layer representing the cell monolayer and porous membrane support is governed by the Brinkman equations:

where ε_p is the porosity, u are vectors of the velocity field within the porous medium and κ is the Darcy permeability of the porous medium (m²).

No universal model of cell secretion currently exists. Therefore, secretion was modeled as a constant local concentration on the upper boundary of the uniform porous layer in order to study the evolution of the concentration profile within the device and the biosensor response. Biomolecules were assumed to be secreted apically from the cell monolayer, pass across the cell monolayer first, followed by the porous membrane before entering the lower channel.

The concentration of secreted analyte was held constant on the surface of the permeable membrane:

where C_S is the apical concentration of analyte over the cell monolayer, L_m is the length of the porous membrane support, h_b is the lower channel height and δ_m is the combined thickness of the cell layer and porous membrane support.

Transport of secreted analyte was assumed to be equivalent to the transport of a diluted species since cellular secretion generally results in low secreted species concentrations. The concentration profile within the device is given by the following equation:

where c_a is the concentration of analyte (M), and D_S is the diffusion coefficient of a typical secreted protein (10^-10 m²s^-129). The secreted biomolecules were assumed to be smaller in size than endothelial cell tight and leaky junctions (∼2-20 nm³⁰). Cell-biomolecule interactions such as charge, steric and surface receptor binding effects were neglected, and thus solute reflection at the upper cell monolayer surface was also neglected. Oncotic pressure-driven transport was assumed to be negligible due to the protein concentration in the medium being nearly identical on both sides of the porous membrane support and cell monolayer.

No flux passes through the walls of the microfluidic device, leading to the following boundary condition on the device walls where n is the outward unit normal vector:

The biosensor on the bottom surface of the lower channel was modeled as a second order reaction representing surface receptor-biomolecule binding interactions given by the following rate equation:^12,31

where b_s is the surface concentration of bound receptors (mol m^-2), b_max is the total surface concentration of receptors, k_on is the association rate constant (M^-1s^-1) and k_off is the dissociation rate constant (s^-1).

The binding of analyte to the biosensor necessitates a boundary condition describing the flux of biomolecules to the biosensor surface:

Fig. 2A shows the upper channel shear stress at transmembrane pressures from 0 to 1,800 Pa for different upper channel height-to-length ratios h_t/L_t. The lower the ratio, the wider the range of transmembrane pressures that may be applied while remaining within physiological shear stress limits. For a given h_t/L_t ratio, the presence of backpressure P_head = 98.1 Pa in the lower channel results in a non-zero upper channel shear stress and flow rate being required to produce ΔP_T = 0, and eliminate backflow from the lower to upper channel. The threshold shear stress τ₀ was defined as the minimum shear stress for which ΔP_T > 0 (derivation found in ESI). Fig. 2B represents the plane of possible upper channel shear stresses that may be generated over a range of h_t/L_t ratios. The threshold shear stress (Eq. S3 of the supplementary material) separates this plane into a region where conditions for biosensing are met since ΔP_T > 0, and a region where biosensing is not possible since ΔP_T < 0 and no analyte transport to the lower channel occurs. By referring to Fig. 2B and selecting an appropriate upper channel h_t/L_t ratio, a wide range of physiological shear stresses may be generated from low 1-5 dyn cm^-2 post-capillary venule shear stress³⁷ to higher 40-60 dyn cm^-2 arteriole/capillary shear stress^2,37 while also ensuring ΔP_T > 0 and enabling biosensing within the device.

These findings lead to the following design guideline that ensures analyte transport to the biosensor and physiological shear stresses: the minimum desired upper channel shear stress, i.e., threshold shear stress τ₀, determines the maximum device h_t/L_t ratio. Coupled with the constraint w_t/h_t = 10, all three upper channel microchannel dimensions may therefore be selected. For example, for a minimum desired upper channel shear stress of τ₀ = 5 dyn cm^-2, the maximum device h_t/L_t is ∼0.005 (Fig. 2B). For a typical L_t = 20,000 μm, this h_t/L_t ratio yields h_t = 100 μm and w_t = 1,000 μm. Furthermore, Fig. 2B may be used to establish the ΔP_T for different levels of physiological shear stress for h_t/L_t between 10^-4 to 5×10^-3. For example, an arterial stress of 35 dyn cm^-2 in a channel with h_t/L_t = 0.005 requires ΔP_T ∼ 600 Pa (Fig. 2B). Upper channel shear stresses within arteriole shear stress values or lower were used in all subsequent simulations and therefore meet the requirement of physiological shear stress.

The properties of the cell monolayer and porous membrane support have a significant effect on apically-secreted biomolecule transport from the upper to lower channel, as the biomolecules must travel through the cell monolayer. In particular, the permeability of this combined layer to fluid and biomolecules is a critical consideration in the context of analyte transport to the biosensor surface. The hydraulic permeability for a bare porous membrane support is ∼6×10^-16 m²²³ whereas that of a fully confluent endothelial cell monolayer^32,33 with tight cell-cell junctions is ∼10^-18-10^-20 m².

Fig. 3A illustrates the numerically-simulated laminar flow profile in the upper channel and vertical flow through the cell monolayer-porous membrane support layer in the microfluidic device. Analyte transported from upper to lower channel is indicated by flux lines as well as the device concentration profile where secreted biomolecules are observed to accumulate in the lower channel. Fig. 3B shows simulated time constants for a range of Darcy permeabilities κ and transmembrane pressures ΔP_T for a membrane length L_m of 10 mm. The time constant is the time required for the analyte concentration to reach 63% of the final saturating concentration, C_S = 1 nM, in the region of interest (ROI). The ROI is defined as the surface of the lower channel as indicated in the inset of Fig. 3B. The biosensor reaction was not modeled in these simulations to examine analyte transport independent from biosensor effects on the local analyte concentration.

For a given ΔP_T, there was a significant monotonic increase in time constant value as the permeability decreased (Fig. 3B). Across the range of ΔP_T simulated, greater ΔP_T led to lower time constants. However, at Darcy permeabilities lower than ∼10^-19 m², increasing ΔP_T became a less effective strategy of reducing the time constant even if increased to 1,600 Pa, which corresponds to the upper end of physiologically-relevant shear stresses.

Fig. 3C illustrates the effect of the porous membrane length L_m on the time constants at low and high transmembrane pressures. Larger porous membrane lengths resulted in lower time constants due to the additional surface area available for analyte transport to the lower channel and the larger ΔP_T gradient lengthwise across the membrane. This effect is more apparent at low ΔP_T where the ΔP_T gradient is negligible for smaller L_m, leading to higher time constants. For example, at ΔP_T = 0 Pa, L_m = 1 mm leads to larger time constants across all permeabilities. However, as ΔP_T increases, the ΔP_T across all L_ms become sufficiently large to transport analyte to the lower channel ROI at comparable rates, leading to similar time constants. At ΔP_T = 1,600 Pa, the time constants for different membrane lengths converge as permeability increases.

These analyses are equally applicable to other multichannel configurations used in microfluidic cell culture, including configurations in which a hydrogel layer (with or without seeded cells) is used instead of a porous membrane support to separate channels.^15,16 Furthermore, the model assumes secretion in the apical direction, but could also be modified to study basal secretion where biomolecule transport is not across the entire cell monolayer thickness, but instead mostly across the porous membrane support.

The simulated range of Darcy permeabilities was determined from hydraulic conductivity values, which are used in well-known models of transendothelial transport such as the Kedem-Katchalsky equations³⁸ to describe convective transport of solute across the endothelium. The simulated Darcy permeabilities, however, do not inherently account for further reductions in solute permeability arising from endothelium size-selectivity and semi-permeability, a limitation of the model presented here. Proteins like albumin (69 kDa) can nevertheless still cross the endothelium and are similar in size to the tight junctions between endothelial cells³⁰ (∼2 nm). Indeed, biomolecules secreted in the vascular environment³ are often smaller than albumin, including endothelin-1 (2.5 kDa), nitric oxide (30 Da) and transforming growth factor-β (25 kDa). Moreover, Young et al. previously experimentally validated an analytical model of transendothelial transport based on Darcy’s law for endothelial permeability measurements²³ where the cell monolayer and porous membrane support were also considered a uniform porous medium.

We studied the relative effects of convective and diffusive transport rates through the porous membrane by considering the pore Péclet number Pe_p, which combines the critical mass transport parameters for porous media flow into a single dimensionless quantity:

A pore Péclet number Pe_p << 1 corresponds to diffusion-dominated transport and Pe_p >> 1 corresponds to convection-dominated transport. A realistic range of possible Pe_p values was determined to be ∼10^-3-10² based on the parameter space explored in Fig. 3. Fig. 4 shows the time constants for the aforementioned range of Pe_p values and different membrane lengths L_m between 0.5-10 mm. Transport is diffusion-limited for Pe_p < 10^-2, and further decreases in Pe_p did not significantly increase the time constant. For convection-dominated transport (Pe_p > 10), further increases in Pe_p did not significantly decrease the time constants. Thus, the greatest reduction in time constant values for a given membrane length occurred in the transition region, with Pe_p values of 10^-2 to 10¹. In order to design for rapid analyte transport, Pe_p > 10 in the convection-dominated region is ideal. However, given the constraints imposed by physiological upper channel shear stress and cell monolayer permeability, it is more realistic to design for Pe_p in the transition region. For example, using Eq. 10, Pe_p = 1, which corresponds to a time constant of ∼700 s, requires ΔP_T ∼1,100 Pa for κ = 10^-18 m², which leads to an upper channel shear stress of 36 dyn cm^-2 for h_t/L_t = 3×10^-3. In terms of selecting a membrane length L_m, for Pe_p > 1, time constants are relatively independent of L_m while for Pe_p < 1, larger L_m lead to lower time constants. Thus, longer L_m improve analyte transport times when lower Pe_p are required.

We next considered the effect of a surface-based biosensor on the local analyte concentration in the lower channel on the complete convection-diffusion-biosensor system. The porous membrane length, surface concentration of receptors and dissociation rate constant were fixed at L_m = 10 mm, b_max = 10^-7 mol m^-2 and k_off = 10^-3 s^-1, respectively. The cell monolayer concentration C_S was set as the equilibrium dissociation constant K_D value where, upon biosensor signal equilibration, 50% of the surface receptors are bound. K_D is an important measure of binding affinity where a lower K_D represents a higher affinity. The biosensor was positioned in the lower channel under the membrane edge nearest the upper channel inlet.

In Fig. 5A, the effect of different Pe_p on the % signal equilibration over time is shown (K_D = 100 nM, k_off = 10^-3 s^-1). The % signal equilibration represents the percentage of surface receptors bound at a given time relative to the total number of receptors bound at equilibrium. For larger Pe_p = 0.5 and 1, representing increased convective transport to the biosensor surface, the % signal equilibration rises to completion more rapidly than when Pe_p < 0.1, again highlighting the important role of analyte transport to the biosensor surface.

One mechanism by which rapid response times may be achieved even at lower Pe_p is shown in Fig. 5B. The % signal equilibration increases rapidly even at a low Pe_p of 10^-2 for c_a/K_D = 10 relative to c_a/K_D = 1, which are ratios of the local concentration relative to the equilibrium constant. Thus, the concentration of secreted analyte in the lower channel, if larger than the K_D value of the receptor, may significantly lower the biosensor signal equilibration time.

The affinity-based biosensor’s response depends on numerous parameters including the biosensor surface-bound receptor properties, analyte transport to the surface and local analyte concentration. We studied the effect of convective-diffusive analyte transport and the surface receptor binding rates on the biosensor response, as represented respectively by the pore Péclet number (Eq. 10) and the Biot number:

A Biot number << 1 indicates a reaction-limited system where the reaction timescale greatly exceeds the diffusive transport timescale. Conversely, a Biot number >> 1 indicates a diffusive transport-limited system, which typically occurs with surface receptors that exhibit a high affinity towards its ligand.

Fig. 6 shows the percentage of signal equilibration across a Biot number range spanning six orders of magnitude (L_m = 10 mm, κ = 10^-18 m²). A timescale of 1 h was selected to compare the kinetics of the system across different conditions. A higher percentage of equilibration corresponds to a lower response time, a particularly important consideration due to the effect of the cell monolayer in reducing convective transport of the analyte to the biosensor surface. A complete list of parameter values used to compute Bi is detailed in Table S1 of the supplementary material. For each Bi, analyte transport rates to the sensor surface were also varied by examining four different Pe_p values between 10^-3 to 10⁰. Furthermore, for each Biot number, simulations were performed for k_off = 10^-3 s^-1 (Fig. 6A) and k_off = 10^-4 s^-1 (Fig. 6B).

With k_off = 10^-3 s^-1, the % equilibration monotonically decreased as Bi increased for a given Pe_p value, with significant decreases occurring at Bi ∼ 10^-1-10⁰. At high Bi, the % equilibration was very low implying long biosensor equilibration times. Furthermore, the % equilibration increased as the Pe_p value increased for a given Bi, indicating that larger analyte transport rates resulted in faster biosensor signal equilibration. For example, the biosensor signal equilibration is at 96% completion for Pe_p = 10⁰ at Bi < 10⁰, but only 40-50% for Pe_p = 10^-3. At Bi > 10⁰, the % equilibration decreases rapidly as Bi increases where, at high Bi, the % equilibration after 1 h is very low even at Pe_p = 10⁰, the highest value simulated. The same trends were also observed when k_off = 10^-4 s^-1 except that the % equilibrations for a given Bi were lower. For example, at Bi = 10⁰ and Pe_p = 10⁰, signal equilibration was at 24% completion in Fig. 6B compared to 96% for the equivalent parameters in Fig. 6A, a ∼70% difference. k_off thus has a significant effect on biosensor equilibration times.

In designing the device-integrated biosensor, a balance is required between biosensor affinity, surface receptor concentration and equilibration time. Constraints are imposed on the surface receptors and analyte transport within the membrane microfluidic device, described by k_off, Bi, and Pe_p. Given that a surface receptor with k_off = 10^-4 s^-1 yields much slower response times than a receptor with k_off = 10^-3 s^-1, one can infer that k_off values ≤ 10^-5 s^-1 will yield prohibitively lengthy response times even at high Pe_p. It should be stressed that these equilibration times are an inherent property of receptors with low k_off values, and are not unique to this system, as the reaction-limited equilibration time scales proportionally to k_off.¹² Unlike a single channel device configuration, however, it is not possible to simply increase the flow rate over the biosensor surface in order to increase convective transport and thereby minimize equilibration times. Moreover, the flow rate is ultimately constrained by the need to apply physiological shear stress. Thus, surface-bound receptors with k_off >10^-4 s^-1 are recommended. Based on the significant drop in % signal equilibration after 1 h for Bi > ∼10⁰ in Fig. 6, Bi within 10^-1 to 10⁰ or less are recommended to reduce biosensor response times in the microfluidic cell culture model. For a lower channel height h_b = 0.1 mm, this allows for combinations of k_on and b_max between 10⁴-10⁶ M^-1 s^-1 and 10^-7-10^-9 mol/m², respectively, for both k_off = 10^-3 and 10^-4 s^-1, which compare favourably to receptor association constant and surface density values found in the literature.^34–36 Receptors with K_D = k_off/k_on values between 1-100 nM and 0.1-10 nM may therefore be used for k_off = 10^-3 and 10^-4 s^-1, respectively. Antibodies possess, on average, K_D values in the nanomolar range and are considered to exhibit a high level of affinity toward the analyte.^39,40 Aptamers for important cell-secreted factors in the vascular environment such as transforming growth factor-beta 1 and tumour necrosis factor-alpha have also been shown to possess K_D values in the nanomolar range (1.07 nM⁴¹ and 7.0 nM,⁴² respectively). Thus, it is possible to use receptors that exhibit high affinity towards their ligands while also achieving rapid biosensor equilibration times. Further reductions in the Biot number may be achieved by reducing h_b. However, unlike k_on and b_max, it is not possible to vary h_b by multiple orders of magnitude and thus only limited reductions in Bi are possible. Similarly, the diffusion coefficient D_s for secreted biomolecules remains on the same order of magnitude as the simulated value of 10^-10 m²/s. In conjunction with the guidelines previously described, it is therefore possible to design for a specific physiological shear stress range while also maximizing cell-secreted factor transport to the lower channel and minimizing the biosensor response time.

We developed a computational model of a microfluidic vascular model integrating flow-induced shear stress and the monitoring of cell-secreted biomolecules using an integrated biosensor. The model was applied to a bilayer configuration that isolates the biosensor from significant analyte dilution, thereby maximizing the concentration of analyte detected, but is readily extended to other multilayer configurations. Using the computational model, we developed guidelines that address both physiological and mass transport constraints as well as enable device operation to be tailored to designer-specified physiological shear stresses while ensuring rapid transport of secreted biomolecules to the biosensor. Special consideration was also given to minimizing the biosensor response time within the mass transport limits of the device configuration. These guidelines may ultimately be used to incorporate physiological shear flow and integrated biosensing into microfluidic vascular models, adding valuable functionality for use in fundamental cell biology and drug screening applications.

See supplementary material for a detailed description of the computational model mesh, derivation of the critical shear stress and biosensor reaction kinetic simulation parameters.

The authors have no conflicts of interest to declare. This work was partially funded by the Natural Sciences and Engineering Research Council of Canada and the Queen Elizabeth II Graduate Scholarship in Science and Technology Program.