Mass Transport Effects in Suspended Waveguide Biosensors Integrated in Microfluidic Channels

Label-free optical biosensors based on integrated photonic devices have demonstrated sensitive and selective detection of biological analytes. Integrating these sensor platforms into microfluidic devices reduces the required sample volume and enables rapid delivery of sample to the sensor surface, thereby improving response times. Conventionally, these devices are embedded in or adjacent to the substrate; therefore, the effective sensing area lies within the slow-flow region at the floor of the channel, reducing the efficiency of sample delivery. Recently, a suspended waveguide sensor was developed in which the device is elevated off of the substrate and the sensing region does not rest on the substrate. This geometry places the sensing region in the middle of the parabolic velocity profile, reduces the distance that a particle must travel by diffusion to be detected, and allows binding to both surfaces of the sensor. We use a finite element model to simulate advection, diffusion, and specific binding of interleukin 6, a signaling protein, to this waveguide-based biosensor at a range of elevations within a microfluidic channel. We compare the transient performance of these suspended waveguide sensors with that of traditional planar devices, studying both the detection threshold response time and the time to reach equilibrium. We also develop a theoretical framework for predicting the behavior of these suspended sensors. These simulation and theoretical results provide a roadmap for improving sensor performance and minimizing the amount of sample required to make measurements.


Theory and Model
Our model involves three simultaneous and coupled processes: fluid flow, mass transport, and a first-order surface reaction. We must solve the time-dependent partial differential equations that govern each physical process. Even for the simplified 2-dimensional geometry described in the main paper, this would be a daunting problem. However, the scope of our study allows us to simplify the calculations required tremendously. Since we are dealing with bulk analyte concentrations in the sub-micromolar regime, the effect of mass transport on the fluid velocity is negligible, and we can split up the simulations into two stages: (1) solve for fluid flow and (2) use this result in solving the coupled transient mass transport and surface reaction equations. This approach leverages the ability of COMSOL Multiphysics to perform iterative simulations. In addition, we limit the range of flow velocities studied to fall within the laminar flow regime, which allows us to use a steady-state solution to Part 1.

Fluid Flow
The steady-state Navier-Stokes equations for an incompressible Newtonian fluid with no volume forces are: where ρ is the fluid density, μ is the fluid viscosity, is the fluid velocity vector, p is pressure, and I is the identity matrix. Equations (S1) and (S2) are expressions of conservation of fluid momentum and mass respectively.

S2
In our model of a 2-dimensional cross-section of the flow cell, we neglect wall effects that might alter the velocity field in the width direction. This is an acceptable approximation when the effective sensing area of the waveguide is situated in the middle of the channel and away from the sidewalls, as it is for the devices of interest, or when the channel width is very large relative to the channel height. We solve Equations (S1) and (S2) simultaneously, subject to the following boundary conditions: (BC-1) At all solid surfaces (i.e., the top and bottom channel walls and the waveguide surface), we assume no-slip: (BC-2) At the channel inlet, we assume unidirectional pressure-driven flow with a parabolic velocity profile given by: where ∆ is the axial pressure gradient, H is the channel height, is the unit vector in the direction of flow, and y is the coordinate in the direction of the cell height, with y = 0 at the cell floor. Integrating across the channel height, we obtain: which can be inserted into (S4) to obtain the boundary condition equation used in our model, an inlet velocity field parameterized by : 6 (S6) (BC-3) At the channel outlet, we assume atmospheric pressure: Solving Equations (S1) and (S2) using the boundary conditions described here gives the steady-state velocity profile and pressure distribution in the channel. These results are shown in Figure S1

Mass Transport and Surface Reaction
The transport of analyte in bulk solution is described by the convection-diffusion equation: where [A] is the analyte concentration, D is the diffusivity of the species, and is the velocity field solved for numerically as in Section 1.1. The binding of analyte to the functionalized sensor surface is assumed to follow a first-order Langmuir model describing the reaction between a freely-diffusing antigen A and an unoccupied binding site B, forming a bound complex C s . This reaction may be expressed as: and its forward (rate f ) and reverse (rate r ) reaction rates are given by: Here, k f and k r are the forward and reverse kinetic rate constants, respectively, and [ ] denotes species concentration. The concentration of unoccupied binding sites [B] is subject to conservation of species according to the expression: where [B] m is the total concentration of binding sites. Combining Equations (S10) through (S12), we obtain the following mass-balance expression for the surface concentration of bound antigen-antibody complex: where A | is the concentration of analyte at the sensor surface. Equations (S8) and (S13) must be solved simultaneously to determine the time-dependent concentration fields. The following initial condition and boundary conditions apply: (BC-2) At the sensor surface, the following reaction/flux boundary condition couples Equations (S8) and (S13). The unit outward normal vector to the sensor surface, , is used to express the flux of analyte leaving solution as: The net mass flux across all other (non-binding) surfaces is zero: There is no diffusive mass flux out of the flow cell (analyte is removed via advection only):

Laminar Flow Assumption
The Reynolds number: is a dimensionless number that arises naturally from the dimensionless form of Equation (S1). It gives an indication of the relative importance of inertial and viscous forces in a fluid system. Here, v 0 is a characteristic velocity, l 0 is a characteristic length scale, and ν = μ/ρ is the kinematic viscosity (dynamic viscosity divided by density) of the fluid. In our case, we define separate Reynolds numbers for the channel and the sensor: and: Re (S19) where is the mean inlet fluid velocity, H is the channel height, and l = w + d is the length of the sensor in the direction of flow (see Figure 2(a) in the main text). Indicative values for the range of flowrates studied are listed in Table S1 below. For fully-developed flow in a channel, laminar flow occurs for Re c < 2,000, and we are well within this domain. We approximate our sensor as a cylinder in order to validate our assumption of laminar flow in the channel and use the literature analysis for flow around such an obstacle. At Re s « 1, flow is orderly, and at Re s ~ 10, a pair of stable vortices appears behind the sensor, but flow as a whole remains laminar [1]. The presence/absence of these vortices should not significantly affect mass transfer to the sensor, since they exist downstream of it.

Model Convergence
The finite element mesh used for this study was composed of triangular elements with sizes constrained by: (S20) and: where and are the maximum allowed mesh element sizes in the flow cell and on the surface of the waveguide, respectively, H is the channel height, d is the diameter of the waveguide, and n v and n s are tunable mesh size parameters. By independently varying n v and n s , we control the quality of the mesh generated. A finer mesh produces a more precise solution, but can require vastly greater computation power and time to solve. The final mesh used in our study is shown in Figure S2 below. The mesh is characterized by n v = 40 and n s = 40, and contains 57,916 elements.

Determining Optimum n v
The dominant processes occurring in the flow cell, far from the waveguide surface, are fluid flow and convective mass transport. Thus, a good criterion for evaluating the effectiveness of the mesh here is the resolution of the laminar flow profile generated. We evaluated the flow velocities along a cutline near the entrance to the flow cell ( Figure S3(a)) while varying n v , and compared these values to the smooth parabolic profile specified by the inlet boundary condition Equation (S6); n s was held constant at 40. As Figure S3(b) shows, the profiles quickly approach the correct shape as n v (and, consequently, the number of nodes along the cutline) increases. We quantify this convergence by calculating the root-mean-square error in each approximation of the flow profile ( Figure S3(c)). Based on this data, we chose n v = 40, which corresponds to ~0.1% RMS deviation from the predicted flow profile.

Determining Optimum n s
We determined the effectiveness of the mesh at the waveguide surface by monitoring convergence to the known equilibrium concentration of bound analyte. At equilibrium, the amount of analyte bound per unit length of the sensor is given by: where: and is the perimeter of the cross-section of the waveguide (~68.1356 μm). Figure S4(a) displays binding data obtained at different values of n s (with n v = 40, as determined in the preceding analysis).
We can see that as we increase mesh resolution, the ratio of the (calculated) contour integral on the left to the (exact) expression on the right converges to 1. The percent deviation: is plotted versus n s in Figure S4(b). Based on this data, we chose n s = 40, which corresponds to ~0.001% deviation from the accurate equilibrium concentration of bound analyte. Figure S4. Verification of n s convergence. (a) Binding data at a range of n s values. Each curve asymptotically approaches an equilibrium amount of bound analyte as t increases.
As n s increases, these equilibrium amounts converge to the analytical value given by the right side of Equation (S22). (b) Percent deviation of the bound analyte concentration from the analytical equilibrium value as a function of n s and time. Data is presented at different times to demonstrate that the degree of convergence can be arbitrarily improved by increasing n s at large enough times.

Convergence to [C s ] eq
We quantify convergence to the theoretical equilibrium concentration by using the method described in Section 2.2. Specifically, we evaluate Equation (S24) at the last time step of each simulation performed. The results are displayed in Figure S5. . The sensor-surface mesh element size has been optimized to resolve binding to the complex suspended waveguide geometry. Therefore, the equivalent mesh is able to easily model binding to the simple planar geometry, resulting in very accurate solutions for these cases.
In most cases, the deviation is ~0.001%; the n s -limited value determined in the convergence analysis of Section 2.2. The highest deviation is ~0.025%, and Figure S4(b) suggests that we could reduce this by simply extending those simulations to later times. Since our metric for equilibration is 95% of [C s ] eq , however, it is unnecessary to do so, as the degree of asymptotic convergence demonstrated is more than sufficient (0.025% « 5%).

Start of Binding Time t 1
We define the start of binding as the time t 1 (and corresponding average surface concentration [C s ] 1 ) when, on average, a single analyte molecule is bound per micrometer length along the sensor. The number of bound analyte molecules per unit sensor length is given by: where N A is Avogadro's constant, and the integral is over the perimeter of the waveguide cross-section.
With N C = 10 6 molecules/m, we get:

Effect of Detection Threshold on Sensor Elevation/Geometry Study Results
As discussed in the main text, we assumed an approximate detection threshold for the suspended waveguide sensor studied here based on previous reports involving similar devices. Figure S7 shows that the qualitative response behavior of the various sensor elevations and geometries remains consistent over a wide range of noise-limited detection thresholds. In all cases, the suspended sensors outperform the embedded ones, and performance increases as the waveguides are located nearer to mid-channel. The quantitative differences in detection times are greatest for low flow velocities and high detection limits.

Slow-Limit Velocity Expression
The average velocity of pressure-driven flow through a channel of height h is given by: For our hypothetical two-channel system, neglecting the thickness of the sensor itself, we have: A mass balance over the whole system gives: (S29) In the slow limit, ∆ is identical in both hypothetical channels (as is μ, of course). Substituting Equation (S27) into (S29): ∆ 12 (S30) Therefore: S11 (S31) which can be rearranged to give Equation (5) using the identity Equation (S28).

Fast-Limit Velocity Expression
In the fast limit, the total flowrate through the region above or below the sensor is found simply by integrating the inlet parabolic velocity profile Equation (S6) over that part of the channel: 6 (S32) Evaluating the integral and rearranging gives Equation (6) directly.