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

Fluid flow and flow-induced shear stress are critical components of the vascular microenvironment commonly studied using microfluidic cell culture models. Microfluidic vascular models mimicking the physiological microenvironment also offer great potential for incorporating on-chip biomolecular detection. In spite of this potential, however, there are few examples of such functionality. Detection of biomolecules released by cells under flow-induced shear stress is a significant challenge due to severe sample dilution caused by the fluid flow used to generate the shear stress, frequently to the extent where the analyte is no longer detectable. In this work, we developed a computational model of a vascular microfluidic cell culture model that integrates physiological shear flow and on-chip monitoring of cell-secreted factors. Applicable to multilayer device configurations, the computational model was applied to a bilayer configuration, which has been used in numerous cell culture applications including vascular models. Guidelines were established that allow cells to be subjected to a wide range of physiological shear stress while ensuring optimal rapid transport of analyte to the biosensor surface and minimized biosensor response times. These guidelines therefore enable the development of microfluidic vascular models that integrate cell-secreted factor detection while addressing flow constraints imposed by physiological shear stress. Ultimately, this work will result in the addition of valuable functionality to microfluidic cell culture models that further fulfill their potential as labs-on-chips. 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 will ensure optimal device which is contingent on


INTRODUCTION
The forces exerted by flowing blood on vascular cells in vivo play critical roles in regulating vascular cell biology 1 . 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 4 .
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 5 .
However, in situ microfluidic biosensing of cell-secreted biomolecules is typically performed under static fluid conditions 6 . Assay integration is an oft-cited advantage of microfluidic cell culture devices [7][8][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][13][14] . This analysis, however, is not sufficient to address the constraints of incorporating physiological levels of shear stress in the device.

Page 4 of 28
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][18][19][20] , and is particularly well-suited for microfluidic vascular models, where it has been used to study cancer cell 21 and monocyte 22 adhesion to an endothelial monolayer, endothelial permeability 23,24 , the vascular/valvular microenvironment 25 and the blood-brain barrier 26,27 . Bilayer culture devices have also been coupled with mass spectrometry for on-line monitoring of drug permeability 28 . 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. Surface-based sensor of length LS is located on bottom surface of lower channel between x = A and x = B.
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.

Computational Model Description
A schematic of the computationally-modeled bilayer microfluidic device is shown in Fig.1B.
Based on a model previously developed by Young et al. 23 , the cell monolayer and porous Page 6 of 28 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.

Governing Equations & Boundary Conditions
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 1. 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 2 contains a list of the design parameters considered in the guidelines developed in this paper. 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, is 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, are vectors of the velocity field within the porous medium and κ is the Darcy permeability of the porous medium (m 2 ).

Cell layer secretion
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 CS is the apical concentration of analyte over the cell monolayer, Lm is the length of the porous membrane support, hb is the lower channel height and δ is the combined thickness of the cell layer and porous membrane support.

Transport of secreted analyte
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 ca is the concentration of analyte (M), and DS is the diffusion coefficient of a typical secreted protein (10 -10 m 2 s -1 29 ). The secreted biomolecules were assumed to be smaller in size than endothelial cell tight and leaky junctions (~2-20 nm 30 ). 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:

Biosensor binding kinetics
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 bs is the surface concentration of bound receptors (mol m -2 ), bmax is the total surface concentration of receptors, kon is the association rate constant (M -1 s -1 ) and koff is the dissociation rate constant (s -1 ).
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.

Design Consideration: Rapid analyte transport
We studied the relative effects of convective and diffusive transport rates through the porous membrane by considering the pore Péclet number Pep, which combines the critical mass transport parameters for porous media flow into a single dimensionless quantity:  For larger Pep = 0.5 and 1, representing increased convective transport to the biosensor surface, the % signal equilibration rises to completion more rapidly than when Pep < 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 Pep is shown in Fig.5B. The % signal equilibration increases rapidly even at a low Pep of 10 -2 for ca/KD = 10 relative to ca/KD = 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 KD value of the receptor, may significantly lower the biosensor signal equilibration time.

Design consideration: Rapid biosensor equilibration times
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.  Table S1. For each Bi, analyte transport rates to the sensor surface were also varied by examining four different Pep values between 10 -3 to 10 0 . Furthermore, for each Biot number, simulations were performed for koff = 10 -3 s -1 (Fig.6A) and koff = 10 -4 s -1 (Fig.6B).
With koff = 10 -3 s -1 , the % equilibration monotonically decreased as Bi increased for a given Pep value, with significant decreases occurring at Bi ~ 10 -1 -10 0 . At high Bi, the % equilibration was very low implying long biosensor equilibration times. Furthermore, the % equilibration increased as the Pep 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 Pep = 10 0 at Bi < 10 0 , but only 40-50% for Pep = 10 -3 . At Bi > 10 0 , the % equilibration decreases rapidly as Bi increases where, at high Bi, the % equilibration after 1 h is very low even at Pep = 10 0 , the highest value simulated. The same trends were also observed when koff = 10 -4 s -1 except that the % equilibrations for a given Bi were lower. For example, at Bi = 10 0 and Pep = 10 0 , signal equilibration was at 24% completion in Fig.6B compared to 96% for the equivalent parameters in Fig.6A, a ~70% difference. koff 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 koff, Bi, and Pep. Given that a surface receptor with koff = 10 -4 s -1 yields much slower response times than a receptor with koff = 10 -3 s -1 , one can infer that koff values ≤ 10 -5 s -1 will yield prohibitively lengthy response times even at high Pep. It should be stressed that these equilibration times are an inherent property of receptors with low koff values, and are not unique to this system, as the reaction-limited equilibration time scales proportionally to koff 12 . 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 koff >10 -4 s -1 are recommended.
Based on the significant drop in % signal equilibration after 1 h for Bi > ~10 0 in Fig.6

CONCLUSION
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

A) B)
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.

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