Delft University of Technology Quantitative chemical biosensing by bacterial chemotaxis in microfluidic chips

pinatubonensis towards the herbicide 2,4-dichlorophenoxyacetate. Our results show that bacterial on a scale of minutes and may be used for developing faster bioreporter assays.


Introduction
Whole-cell living bacterial bioreporters are useful and versatile tools to detect specific chemicals in environmental (Kohler et al., 2000;van der Meer and Belkin, 2010), food (Baumann and van der Meer, 2007;Daszczuk et al., 2014) or medical samples (Turner et al., 2007;Courbet et al., 2015). Bioreporters consist of engineered microorganisms or cell lines equipped with a synthetic gene circuit to produce an easily measurable signal in response to a single or group of related chemical targets (Daunert et al., 2000;van der Meer et al., 2004). Target detection by bioreporters is achieved by exploiting the huge natural diversity of sensory proteins and metabolic pathways, or mutants thereof (van der Meer et al., 2004;Park et al., 2013). The use of bioreporters for detection of chemical assaults provides an alternative to classical analytical methods with the advantage that they assess the bioavailable fraction of the toxicant. Bioreporter assays are relatively fast, cheap to produce and flexible, adaptable to in situ conditions as well as to portable analytical devices (van der Meer and Belkin, 2010;Park et al., 2013;Lim et al., 2015). Conventional whole-cell bioreporters are based on the induction of the expression of a reporter gene (e.g., fluorescent protein, beta-galactosidase or luciferase) when cells are exposed to a specific ligand (Daunert et al., 2000;Harms et al., 2006). The induction process typically requires between 30 min and several hours; the total time needed for de novo transcription upon detection of the signal [several minutes (Gotta et al., 1991;Vogel and Jensen, 1994)], translation [minutes (Proshkin et al., 2010)], maturation of proteins [minutes to hours (Iizuka et al., 2011)] and accumulation of sufficient reporter protein to provide a signal (hour-scale). For handheld devices or for inclusion within robots, assays based on regular reporter gene expression may take too long. It would, thus, be interesting to explore alternative physiological responses, which may show faster reaction to the presence of specific compound, and which do not require de novo transcription-translation-maturation. An interesting candidate is bacterial chemotaxis, where cells react on a time scale of hundreds of milliseconds (Sourjik and Wingreen, 2012) to chemical gradients.
Chemotaxis is the behaviour by motile bacteria to sense their environment and swim in the direction of or away from sources of chemicals. Chemotaxis has been intensively studied in Escherichia coli (Baker et al., 2006;Hazelbauer, 2012;Sourjik and Wingreen, 2012;Tindall et al., 2012), but various components of the system (e.g., chemoreceptors, flagellar operation) are conserved among motile bacteria (Porter et al., 2011). E. coli motility is characterized by alternating periods of straight runs and tumbling, during which cells randomly reorient their swimming direction. In a uniform environment, E. coli cells swim stochastically with alternating 'runs' ($ 1 s) and 'tumbles' ($ 0.1 s), in order to explore the maximal accessible space (Sourjik and Wingreen, 2012). In presence of an attractant, the cells bias their swimming direction towards the highest concentration of attractant, with longer runs and fewer tumbles (Sourjik and Wingreen, 2012). The periods of runs and tumbles are determined by the rotation direction of the flagella: when flagella turn counterclockwise, they form a bundle and the cell goes straight ('runs'). In contrast, when a flagellum turns clockwise, it leaves and destabilizes the bundle, and the cell changes its direction ('tumbles'). The rotation direction of flagella is controlled by a molecular pathway triggered by the binding of the attractant (or repellent) to a specific chemoreceptor at the cell surface, which via protein-protein interaction and phosphorylations is transmitted to the flagellar motor [for review, see (Wadhams and Armitage, 2004)]. Interestingly, some bacteria show chemotaxis to compounds of environmental interest, which they use as carbon and energy sources (Parales and Harwood, 2002;Parales et al., 2015). An example of such a bacterium is Cupriavidus pinatubonensis (formerly C. necator or Ralstonia eutropha) JMP134, which is chemotactic towards the herbicide 2,4-dichlorophenoxyacetic acid (2,4-D), that can be metabolized by enzymes encoded on the plasmid pJP4 (Fisher et al., 1978;Don and Pemberton, 1985;Hawkins and Harwood, 2002).
In order to potentially exploit natural chemotactic responses for biosensing of chemical compounds, a quantitative readout of chemotaxis would be imperative, ideally leading to some simple relation between the intensity of the chemotaxis readout and the chemical target concentration. Chemotaxis has been traditionally studied using different techniques such as swimming plates, capillary assays (Adler, 1966;Mesibov and Adler, 1972) or agarose plug assays (Armstrong et al., 1967), which mostly give a qualitative measurement of the chemotaxis response and take several hours or days for results. More recently, microfluidic chips have been deployed to study and quantify chemotaxis (Ahmed et al., 2010a,b;Wessel et al., 2013;Lim et al., 2015;Rothbauer et al., 2015), which could potentially enable faster and quantitative readouts of chemotactic biosensing. The crucial aspect for inducing cellular chemotaxis in a microfluidic chip is the creation of a microscale gradient of the target chemical, which can be formed either by flow or by diffusion (Kalinin et al., 2009;Ahmed et al., 2010a,b). Flow-based gradient generator chips deploy the laminar flow in parallel channels with either attractant or buffer to create a concentration gradient perpendicular to the flow as a consequence of molecular diffusion. This type of chip creates rapid and stable gradients in which chemotaxis is observed, but the cells undergo a shear stress induced by the flow (Mao et al., 2003;Lanning et al., 2008;Englert et al., 2009). In the second type of microfluidic chips, the attractant is allowed to diffuse from a source to create a concentration gradient in which cells can react, but the source is separated from the cell compartment by a physical barrier or hydrogel (Stocker et al., 2008;Kalinin et al., 2009;Ahmed et al., 2010a,b;Si et al., 2012).
Although chemotaxis has been extensively quantified and modelled, it has rarely been proposed as means for biosensing of target molecules, except for a recent study using chemotaxis fluorescence resonance energy transfer (FRET) (Bi et al., 2016). The main aim of this study was, thus, to test whether chemotactic movement of bacterial cells can be reliably quantified and calibrated as a function of attractant concentration. We quantify chemotaxis using a microfluidic chip consisting of three parallel channels, conceptually similar as those described in Kalinin et al. (2009) and Ahmed et al (2010a,b). In our chip, the parallel channels are linked by pores so shallow that cells cannot pass (Fig. 1A), using a procedure described recently (Buffi et al., 2016). Solutions of attractant and buffer are flowing in each of the side channels to create via diffusion a stable gradient across a middle observation channel, where the cells are introduced (Fig. 1B). Cell accumulation on the side of the observation channel closest to the source is recorded and quantified by epifluorescence microscopy. We tested the concept by measuring chemotactic response of E. coli towards its best known attractants (serine, aspartate and methylaspartate) at different concentrations. Results were compared to numerically modelled steady-state transport of the attractant and biomass in the microfluidic chips, and by kinetic modelling of cell-agent-biased swimming in attractant gradients. We further used a strain of C. pinatubonensis JMP134 that was reported to be chemotactic towards 2,4-dichlorophenoxyacetic acid to investigate whether chemicals of environmental concern might be detected using a chemotaxis-microfluidic biosensor chip concept.

Fabrication of a stable gradient-forming microfluidic chip
A polydimethylsiloxane (PDMS) chip was designed in which a stable solute concentration gradient can be established by flowing buffer with and without attractant ('source' and 'sink', respectively) in two side-channels ( Fig. 1A and B). The side-channels are connected to the observation channel by series of perpendicular and shallow pores ($600 nm height and 5 mm wide, Fig. 1A and B), allowing attractant diffusion from the source to the sink channel, which creates a gradient perpendicular to the observation channel where cells can be introduced at low flow rate (Fig. 1B). Modelled pressure and flow distribution across the chip in steady-state conditions (see Experimental methods) shows that the flow in the observation channel is mostly insulated from the side channels, apart from a small inflow through the pores (Fig. 1C, see streamlines).
The generation and stability of the chemical gradient was verified by addition of rhodamine (1 or 10 mM) to the source channel and measuring the fluorescence profile across the observation channel near the entry port. The diffusion coefficient of rhodamine in water [D 5 4.27 3 10 210 m 2 s 21 (Culbertson et al., 2002)] is in the same order of magnitude as that of small attractant molecules (e.g., serine: 8.8 3 10 210 m 2 s 21 and aspartate: 8.0 3 10 210 m 2 s 21 ) (Longsworth, 1955), and can serve as a proxy for the formation of attractant gradients in the experiments with bacterial cells.
Without flow in the observation channel, rhodamine was detected near the sink channel after 5 min and persisted for more than 1 h under constant flow (0.25 ml min 21 ) in the side channels (Supporting Information Fig. S1A and B, red lines). The introduction of (fluorescently labelled) bacterial cells under slow flow ($3 nl min 21 ) in the observation channel slightly changed but did not perturb the stability of the gradient (Supporting Information Fig. S1A, B, blue lines). The modelled chemical gradient, based on the steady-state flow field as presented in Fig. 1C, resembled quite closely the observed rhodamine gradient (Supporting Information Fig. S1C), although we noticed that some rhodamine absorbs to the material of the chip, slightly increasing background fluorescence compared to the model. Notice, however, how the modelled concentrations at the source border of the observation channel do not completely attain the dosed source concentration in the side flow channel, due to a concentration drop across the filter pores (Supporting Information Fig. S1D). We concluded that introducing a chemical in the source channel would rapidly ($20 min) lead to the generation of a stable gradient into which the bacterial cells could be introduced with minimal disturbance. We, thus, expected that the chemotactic cells would react to the chemical gradient by accumulating near the source. Chemotaxis of E. coli MG1655 toward serine To quantify chemotaxis as a function of different source concentrations in the developed PDMS chip, we measured the distribution of E. coli MG1655 cells (Table 1) tagged with a constitutively expressed enhanced green fluorescent protein (MG1655-gfp) across the channel transect and in a 100 mm-wide zone (around position 2200 mm in the chip, Fig. 1C) at different times after establishing the serine gradient. A DfliC nonmotile mutant of E. coli tagged with mcherry (DfliC-mcherry) was used as internal control for nonchemotactic cell distribution as a result of the flow inside the observation channel (Table 1). MG1655-gfp and DfliC-mcherry cells were added as a mixture to the reservoir of the observation channel of the chip and introduced at low flow (3 nl min 21 ) after serine gradients had established (20 min after starting flow in the side channels).
As expected, MG1655-gfp but not DfliC-mcherry cells accumulated near the source channel when serine was added, whereas in absence of serine, MG1655-gfp cells distributed homogenously across the channel ( Fig. 2A, B and D and Supporting Information Fig. S2). Both the relative accumulation and the time needed for accumulating MG1655-gfp cells depended on the serine concentration ( Fig. 2A and C). The relative cell accumulation was quantified by calculating a chemotaxis index (CI), which represents the proportion of fluorescence from cells in the 100 mm zone closest to the source compared to the total fluorescence of cells across a channel transect of 100 mm 3 600 mm (Supporting Information Fig. S3). The CI was slightly, but significantly higher for a source with 1 mM serine compared to buffer alone (p 5 9.5 3 10 217 , pair-wise t-test, one tail, equal variance). In contrast, the CI increased already after 5 min at a source concentration of 10 mM (Fig.  2C). Interestingly, the CI increased slower with a 0.1 mM serine source but was higher after 1 h compared to 10 mM, at which point 50% of the MG1655-gfp cells had accumulated in the 100 mm segment closest to the source channel. At 1 mM serine in the source, the final CI reached was the same as for 10 mM, but the cells took even longer to accumulate (almost no response before 25 min, Fig. 2C).

Chemotaxis simulation under steady-state and kinetic conditions
Mathematical simulation of steady-state cell biomass distributions in the observation channel, using the chip geometry and flow conditions, and including both chemotaxis as well as serine metabolism, largely confirmed the observed results (Fig. 3). Cells without modelled chemotactic ability in presence of 10 mM serine in the source channel formed a band of slightly higher biomass concentration in the middle of the channel without any accumulation near the source (Fig. 3C, 10 mM ser, DfliC), similar to the experimentally observed DfliC mutant behaviour (Fig. 3B, t 5 50 min, DfliC). This centering in the observation channel is the result of the incoming flow from the side channels (Fig. 1C), pushing nonswimming cells away from the filter pores. In contrast, modelled cells given chemotactic ability accumulated near the source channel, which was dependent on the serine concentration (Fig. 3B, t 5 50 min, 10, 100 mM and 1 mM ser; Fig. 3C, steady state 10, 100 mM and 1 mM ser). One can see that the serine concentration gradient is not the same at every position along the length of the observation channel as a result of the presence of the connecting pores and convective transport along the channel (Fig. 3A). As a consequence, the accumulation of cells also changes along the channel, but is most pronounced near the beginning ( Fig. 3B and C).
Measured and modelled steady-state biomass accumulation in the 100 mm wide zone across a channel transect at 2200 mm from the inlet (as indicated in Fig. 3A), were in close agreement (Fig. 4A). Cells indeed responded slightly at 1 mM serine in the source channel, with a maximum chemotaxis response at around 100 mM (Fig. 4A). The modelled chemotaxis index response curve was slightly sensitive to variations in two key chemotaxis model parameters, the serine dissociation constant (K C , optimum 25-30 mM) and the chemotactic sensitivity coefficient (v 0 : optimum value $8 310 24 cm 2 s 21 , Fig. 4A).
In order to better capture and explain the observed time response differences of E. coli cells to increasing serine concentrations ( Fig. 2C), we simulated numerically the chemotactic behaviour of individual cell agents over time as a result of biased random walk in a grid (similar to the inner microfluidic channel), exposed to the calculated steady-state serine gradient (continuum COMSOL simulations, e.g., Fig. 3A). Simulations were based on a Metropolis Markov chain (see Experimental methods). Indeed, kinetic simulations using similar chemotactic parameterization of the cells as in the steady-state modelling, confirmed both the time response as well as the final magnitude of accumulation as a function of serine source concentration ( Fig. 5A-C, compare to Fig. 2A and C). As can be seen from simulation results in Fig. 5B at 0.1 mM serine, cells closest to the source will tend to accumulate faster (red agents in Fig. 5B), but as a result of biased random walk, a small fraction of cells will remain distributed across the channel even after longer periods.

Comparison of chemotaxis of MG1655 towards aspartate and its nonmetabolizable analogue methylaspartate
Because serine is not only attractant for E. coli but can also be metabolized by the cells, its concentration gradient may change over time. Modelled serine uptake on top of chemotactic behaviour suggested that this would slightly reduce the concentration gradient but should have very little effect on the steady-state biomass distribution across the channel (Supporting Information Fig. S1E). In order to test experimentally whether attractant metabolism may have an effect on the chemotactic response time and accumulation of cells near the source, we compared E. coli distributions in the PDMS chip to sources of aspartate (Asp) or its nonmetabolizable analogue, N-methyl-D-aspartate (MeAsp) at 1 mM, 10 mM, 0.1 mM, 1 mM and 10 mM concentrations (Fig. 6).
Compared to serine, the CIs of E. coli to Asp and MeAsp in general were lower (i.e., less cells accumulating near the source, Fig. 6). E. coli MG1655-gfp responded slightly to 1 mM Asp but more rapidly to 10 mM and 0.1 mM, after which the reaction time increased again with 1 and 10 mM (Fig. 6A). The maximum CI increased from 0.19 (at 1 mM) to 0.21 (at 10 mM) and 0.32 (at 0.1 mM, 1 and 10 mM; Fig. 6A). In contrast, exposure of cells to 1 mM MeAsp did not elicit a reaction, whereas the CI at 10 mM MeAsp (0.18) was slightly lower (p 5 0.00078, pair-wise t-test, time points 30-60 min, one tail, equal variance) than at 10 mM Asp  6B). At source concentrations of 0.1 and 1 mM, however, the response of the cells to MeAsp was faster than at the equivalent Asp concentration, whereas at 10 mM MeAsp the response dropped to the level of 10 mM MeAsp (Fig. 6B). This suggests that cells are slightly more sensitive to Asp than to MeAsp, but that attractant metabolism can play a role at higher concentrations (0.1 mM-10 mM), leading to a delay of chemotaxis movement.
As expected the DfliC-mcherry mutant of E. coli did not react to Asp or MeAsp at any of the tested concentrations ( Fig. 6A and B).

Chemotaxis response of Cupriavidus pinatubonensis JMP134 towards 2,4-dichlorophenoxyacetic acid
In order to determine whether chemotactic cell movement can be used as a biosensor readout for pollutants, we tested C. pinatubonensis JMP134 (Table 1), which has been shown to be attracted by the herbicide 2,4-D (Hawkins and Harwood, 2002). Similar to E. coli, C. pinatubonensis JMP134-gfp cells accumulated over time near the source inlet in the PDMS chip and as a function of the 2,4-D concentration (Fig. 7A). However, the concentrations of 2,4-D needed to elicit observable JMP134-gfp cell accumulation were much higher than in case of E. coli for serine, Asp or MeAsp. C. pinatubonensis JMP134-gfp cells became more responsive over time upon storage in the cell reservoir of the microfluidic chip (Fig. 7A, no vs. 1h adaptation). In the first exposure, the CI of C. pinatubonensis JMP134 increased linearly over time and with increasing concentrations of 1, 5 and 10 mM, being significantly higher after 20 min compared to 0 and 0.1 mM 2,4-D (Fig. 7A, no adaptation). In the second exposure, the cells reacted within 5-10 min and the CIs at 1, 5 and 10 mM reached more or less stable values of 0.175, 0.180 and 0.190, respectively (Fig. 7A, 1 h). CI values were reasonably well explained using the steady-state continuous model (implemented in COMSOL) with a K C of 1.5 mM and v 0 of 6 3 10 25 cm 2 s 21 (Fig. 4B).
In order to better understand possible chemotactic differences, we analyzed and compared individual C. pinatubonensis and E. coli cell trajectories calculated from 13 s movie series in zones near the source and the sink channels immediately after introducing the cells in the observation channels (Fig. 7B, 0 min) and 10 min afterwards (Fig. 7B, 10 min). Cells were categorized as 'swimmers' if the angle of their average trajectory length deviated more than 108 from the flow direction and if the ratio of the effective travel distance divided by the total trajectory length was smaller than 0.9; else they were considered 'nonswimmers'. The proportion of swimmers of C. pinatubonensis JMP134-gfp increased up to twofold within 10 min in zones close to the attractant compared to the sink, for 1, 5 and 10 mM 2,4-D but not for 0.1 and 0 mM (Fig. 7B, left panel). In contrast, the number of nonswimmers was similar in both zones over time (ratio $1, Fig. 7B, right panel). Also E. coli showed a sharp increase within 10 min in the number of swimmers close to the source compared to the sink, up to a level of sevenfold for serine and MeAsp, and twofold for Asp (all at 0.1 mM source concentration), whereas the number of nonswimmers remained equal in both zones (Fig. 7B). This, thus, suggested that even when C. pinatubonensis cells were proportionally swimming as much as those of E. coli (ratio 5 2 at 10 mM 2,4-D and 0.1 mM Asp, Fig. 7B), they accumulated much less near the source (CI at 10 mM 2,4-D 5 0.18, Fig. 4B; vs. 0.30 at 0.1 mM Asp for E. coli, Fig. 4C).

Discussion
In this study, we assessed whether bacterial cell accumulation in response to a chemical gradient as a result of chemotactic attraction can be properly quantified for its usage as a biosensor. We deployed a continuous stream of cells in a microfluidic chip that enabled to produce stable chemical gradients under flow, and quantified cell accumulation over time from fluorescence measurements of labelled cells. In order to calibrate the biosensor response, we used E. coli MG1655 chemotaxis towards serine, Asp or MeAsp at different attractant concentrations (between 1 mM and 1 mM). The strongest accumulation was observed after 30 min at 0.1 mM serine (in the source channel), CIs using the Metropolis kinetic individual cell-agent simulation and steady-state continuous modelling (chi_0 and K C are the chemotaxis coefficients and have units as in Fig. 3). B. Experimental and steady-state continuous modelled CI for C. pinatubonensis and 2,4-D. C. Cell accumulation rates for E. coli and the three attractants in different concentrations, calculated from the slope of the chemotaxis kinetic profiles. [Colour figure can be viewed at wileyonlinelibrary.com] whereas the fastest response was obtained after 5 min with 10 mM serine (Fig. 2C). The response of the cells to Asp and MeAsp was slightly less sensitive than to serine with observed accumulation between 10 mM and 10 mM source concentration (Figs 2C and 6). In principle, both the steady-state cell accumulation (i.e., after 40-50 min, Fig. 4A) and the cell accumulation rate (Fig. 4C) as a function of source concentration may be deployed as a calibration curve to interpolate the response to unknown samples. Steady-state cell accumulation from chemotaxis results in a bell-shaped function with an optimum (for E. coli and serine) at around 100 mM serine and a decline at higher concentrations (Fig. 4A), which is in contrast to typical saturation-type output of gene reporter-based biosensors.
The concept of our microfluidic chip was similar as previous designs (Kalinin et al., 2009;Ahmed et al., 2010a,b), but altered by including an array of small and shallow channels, which allow chemical diffusion from the source to a sink channel, while separating the bacterial cells in an observation channel (Fig. 1). The use of shallow filters prevents passage of the cells into the source channel and abolishes the need for hydrogels or membranes as physical barrier (Kalinin et al., 2009;Ahmed et al., 2010a,b), which facilitates the fabrication of the chips. The bacterial cells enter the observation channel at low flow, minimizing shear stress. A slight drawback of the filter structures is that they allow a small part of the flow in the side channels to enter the observation channel, which pushes (nonswimming) cells away from the walls (Fig. 1C). Swimming chemotactic bacteria easily overcome the side-channel inflow and accumulate near the source channel when attractant is present (Fig.  3B). However, a further consequence of this design is a nonhomogeneous attractant gradient and nonuniform cell accumulation along the observation channel, with the highest cell accumulation near the channel's first cross-sections (Fig. 3). Since there is a constant slow influx of new cells, accumulation is reversible when the source concentration alters. Multiple consecutive cell accumulation tests can also be conducted on the same chip by temporarily increasing the flow in the observation channel to wash away the accumulated cells and then resume normal flow, for as long as the cells remain chemotactically active (in their motility buffer in the chip reservoir, Fig. 6C).
The observation that both the magnitude and the kinetic response of the cells varied with attractant concentration in the source, and do not have the same optimum, is intriguing (Fig. 2C). We confirmed by both continuous biomass and kinetic individual cell-agent modelling that steady-state cell accumulation near the source can be explained by two chemotactic parameters: the receptor-ligand dissociation constant (K C ) and a chemotaxis specificity coefficient (v 0 ) as previously proposed (Chen et al., 1998) (see Experimental methods). The steady-state continuous biomass model slightly underestimated the cell accumulation at intermediate (10 mM) and high serine concentrations (1 mM), whereas the kinetic model predicted visible accumulation at low (1 mM) concentrations. The cell accumulation at 1 mM source may not have been detectable experimentally as a consequence of the small medium inflow through the pores (Fig. 5A). Kinetic simulations required an adaptation term a, which is dependent on the concentration gradient (see Experimental methods), suggesting that the delay in response at higher attractant concentrations (0.1 and 1 mM serine) is due to cells not immediately detecting the chemical gradient because of the high local attractant concentration. Attractant metabolism by the cells is further influencing the concentration gradient (Supporting Information Fig. S1E), and particularly at higher concentrations, as our comparison of the responses of E. coli MG1655 cells towards Asp and its nonmetabolizable analogue MeAsp attest (Fig. 6).
In comparison to the E. coli response to serine, Asp or MeAsp, the C. pinetubonensis response to 2,4-D was much lower (Fig. 5B). Fitting C. pinetubonensis CIs by the steady-state continuous model required a two orders of magnitude less sensitive chemotaxis (v 0 of 6 3 10 25 instead of 8 3 10 24 cm 2 s 21 for E. coli, and a K C receptorligand dissociation constant of 1.5 mM instead of 30 mM, Fig. 4B, Supporting Information Table S1). This suggests that the type of chemotaxis by C. pinetubonensis JMP134 to 2,4-D is different than that of E. coli to serine, Asp or MeAsp. Perhaps C. pinetubonensis movement towards 2,4-D is the result of energy-taxis rather than chemotaxis with a genuine receptor for 2,4-D, or the result of an intermediate of 2,4-D metabolism being recognized intracellularly (Ni et al., 2013). The K C -value of C. pinetubonensis to 2,4-D found here (1.5 mM) is in the same range as other reported K C -values for P. putida F1 and toluene ($1.0 mM) (Wang et al., 2015), which may be representative for the same type of chemotaxis. On the other hand, measurements of individual cell movements suggested that C. pinetubonensis cells are on average as motile as E. coli.
Chemotaxis has attracted wide interest and cell movement has been frequently studied in (a variety of) flow-or diffusion-based microfluidic chips, although rarely for the purpose of potential biosensing itself. Diffusion-based microfluidic chips allow specific accumulation of bacteria in a specific zone close to the source of attractant (Stocker et al., 2008;Si et al., 2012;Wang et al., 2012;Crooks et al., 2015;Wang et al., 2015), and, as we demonstrate, cell accumulation is a reasonable quantitative readout of chemotaxis. Our results are in good agreement with measurements performed by others with E. coli RP437 in a different type of microfluidic chip, but without kinetic observations (Kalinin et al., 2009). Our used model values for both K C and v 0 parameters of E. coli and C. pinatobunensis adhere to previous experimental data (Wang et al., 2015). As an alternative quantification method, one could compare distributions of a motile to a nonmotile strain, as a few groups have illustrated (Mao et al., 2003;Lanning et al., 2008;Englert et al., 2009). Also the abundance of individual cell trajectories extracted from long exposure images or movies of fluorescent bacteria has been used as a parameter to characterize chemotaxis response (Fig.  7B) (Englert et al., 2009;Kalinin et al., 2009), although this requires more complex data treatment. Finally, instead of measuring the cellular accumulation, one could infer chemotactic response from FRET between, for example, labelled CheY-YFP and CheZ-CFP in E. coli cell suspensions in stop-flow cuvettes (Sourjik and Berg, 2002;Sourjik and Berg, 2004).
Collectively, our results demonstrate that cell accumulation in microfluidic chips can be quantified as a biosensor. For environmental purposes it is not so exciting to use E. coli to measure serine concentrations, but we demonstrated with C. pinetubonensis and others with P. putida F1 (Wang et al., 2015) that it should be possible to deploy other strains that naturally swim towards toxic compounds. While the response of C. pinatubonensis to 2,4-D was not sensitive enough (1-10 mM range) for a competitive 2,4-D biosensor, our setup enabled to verify the sensitivity of the response. For chemotaxisbased biosensors, one would, thus, ideally need strains with K C and v 0 values in the same range as those of E. coli for serine, in which case cells will accumulate visibly at mM source concentrations. Chemotactic response and readout in the microfluidic chip was relatively rapid (5-10 min, 10 mM serine, Fig. 3C), which is faster than the readout of typical de novo gene induction bioreporters (30 min or longer). We also showed that multiple subsequent source samples can be detected by washing away the accumulated cells (Fig. 6C), although the design of this particular microfluidic chip is not really suitable for multisample analysis. Sample exchange in the current chip is not practical because of the large channel dead volumes, but the microfluidic structure could be modified to include a second PDMS layer with valves that facilitate sample manipulation and reduce dead volumes (Buffi et al., 2016). Microscopic detection of accumulating cells as in its current embodiment could be improved by light scattering or fluorescence detectors within or close to the chip, as shown by Truffer et al. (2014). We conclude that quantitative detection of chemotactic attraction through cell accumulation is a promising alternative to gene reporters for faster whole cell biosensors.

Bacterial strains and growth conditions
A specific motile strain of E. coli MG1655, obtained from the E. coli Genetic stock center (Yale) (CGSC#8237), was transformed with plasmid pME6012-P tac -gfp (which enables constitutive gfp expression from the tac promoter, Table 1). This strain, MG1655-gfp, was grown in M9 minimal medium at 378C with 180 r.p.m. shaking supplemented with 4 g l 21 of glucose, 1 g l 21 of Bacto TM casamino acids (BD difco), Hutner's trace metals (Gerhardt and Murray, 1981) and 1 mM of MgSO 4 and 25 mg ml 21 of tetracycline (medium hereafter called M9-Glc-Tc25) (Salman and Libchaber, 2007). A nonflagellated E. coli MG1655 DfliC mutant [obtained from the Keio collection (Baba et al., 2006)] transformed with plasmid pME6012-P tac -mcherry was used as a control for the effect of flow in the channel (Table 1). Strain DfliC-mcherry was cultured in M9-Glc-Tc25 and additional kanamycin at 50 mg ml 21 .

C. pinatubonensis JMP134 tagging with gfp
A derivative of C. pinatubonensis JMP134 was produced carrying a single chromosomal copy of a P tac -gfp transcriptional fusion delivered by mini-Tn5 transposition. The minitransposon was delivered into C. pinatubonensis JMP134 in a triparental mating (Ditta et al., 1980). C. pinatubonensis transconjugant colonies were checked for green fluorescence and motility and were regrown in MM with 5 mM 2,4-D as sole carbon and energy source. Three individual and independent colonies of C. pinatubonensis JMP134-gfp (with possible different insertion positions of the mini-Tn5) were then stored at 2808C. Only one of the clones (strain 5198, Table 1) was used for the chemotaxis assays but all the clones showed similar chemotaxis towards 2,4-D in agarose plug assays.

Chip design and fabrication
A microfluidic chip was designed using CleWin4 software and fabricated in PDMS (polydimethylsiloxane). The chip had three parallel channels linked by shallow pores of 5 mm wide, 100 mm long and 650 nm high ( Fig. 1A and B). Buffer (sink) and attractant (source) are flowed in separate side channels, each 1 mm wide and 15 mm high. Constant flow and attractant diffusion through the pores lead to the creation of a gradient perpendicular to the observation channel (600 mm wide and 15 mm high) where the cells are introduced (Fig. 1A and B).
A silicon mold of the designed chip was produced in the cleanroom by standard soft photolithography and etching processes (Buffi et al., 2016). A first step of photolithography and etching was performed to create an imprint for the pores, because these have a different height. In a second photolithography step, the imprints for the flow channels were created (Supporting Information Fig. S4). The chip was then fabricated as a negative imprint on the silicon mold using PDMS. Syglard 184 silicone elastomer base and curing agent (Dow Corning) were mixed in a 10:1 ratio, degassed under vacuum and poured onto the silicon mold. The polymerization was performed at 808C overnight. The chips (a block of approximately 2 3 2 cm) were cut and peeled from the silicon mold, after which flow channel inlets were punched with a 1.5 mm diameter Harris Uni-Core puncher (TED PELLA, Inc). Finally, the PDMS chip was bonded to a standard glass microscopy slide (1 mm thickness, RS France, Milian) using oxygen plasma treatment for 6 s (FEMTO plasma cleaner from Diener Electronic, settings: 0.6 mbar, 100 W, Supporting Information Fig. S4)

Gradient formation
Formation of the chemical gradient in the PDMS chip was tested by using 1 or 10 mM of Rhodamine B in water as source and water alone as sink fluid, introduced at a flow rate of 0.25 ml min 21 . Fluorescence was recorded by epifluorescence microscopy in a 100 mm-wide zone around position 2200 mm of the observation channel ( Fig. 3A and B) every 2 min during 20 min (exposure times 110 ms). After 20 min, an MG1655-gfp cell suspension ($10 8 cells ml 21 ) was connected at a flow rate of 3 nl min 21 . Rhodamine fluorescence was again recorded immediately, 5, 10, 15 and 20 min after connecting the cell suspension (exposure times 200 ms). Fluorescence was recorded using a DFC 350 FX R2 Leica camera mounted on an inverted DMI 4000 Leica microscope using a N PLAN 103 objective. Fluorescence profiles were calculated from digital images using the ImageJ open source software (http:// imagej.nih.gov/ij/).

Preparation of cells for chemotaxis assays
Overnight cultures of MG1655-gfp and DfliC-mcherry in M9-Glc with appropriate antibiotics were diluted 100-fold in the same medium and grown 3 h until reaching early exponential phase (culture turbidity at 600 nm (OD 600 ) of 0.5). After checking the motility of the cells under phase-contrast microscopy, 1 ml of each culture was centrifuged at 2400 g for 5 min, after which the supernatant was decanted and the cell pellet resuspended in 1 ml of motility buffer (motility buffer is 10 mM potassium phosphate, 0.1 mM EDTA, 10 mM lactate, 1 mM methionine, pH 7.0) (Berg and Brown, 1974). Centrifugation and resuspension were repeated twice, after which both MG1655-gfp and DfliC-mcherry pellets were, respectively, resuspended in 250 ml and 500 ml of motility buffer. Both cell suspensions were mixed at a 1:1 ratio (250 ml: 250 ml) for introduction into the microfluidic chip. Cultures suspended in motility buffer were used for chemotaxis experiments within 1 h. For chemotaxis assays with C. pinatubonensis JMP134-gfp, an overnight culture grown on MM with 10 mM succinate plus 0.5 mM 2,4-D and 50 mg ml 21 kanamycin was diluted 16-fold in MM with 5 mM succinate and 0.1 mM 2,4-D. The culture was incubated at 308C until its turbidity reached an OD 600 between 0.5 and 0.6 (exponential phase). This cell suspension was then directly used for the chemotaxis assay in order to preserve maximal cellular motility (centrifugation or filtering caused too much loss of swimming cells).

On-chip chemotaxis assays
Two Hamilton 2.5 ml glass syringes (model 1002 RN with an RN(22/51/3) large needle, Hamilton) were connected to $40 cm long tygon tubings (internal ø 0.51 mm, wall 0.91 mm, Tygon ST, Ismatec) and a custom-made metal adapter (length $5 cm, bend to 908, made from steel capillary internal ø 0.25 mm, external ø 1/16 inches, Metrohm). One syringe was filled with motility buffer containing attractant (e.g., 10 mM serine), whereas the other was filled with motility buffer (for E. coli). In case of chemotaxis assays with C. pinatubonensis we used MM instead of motility buffer. A 100 ml Hamilton glass syringe (model 1710 RN with a RN22/51/3 small needle, Hamilton) connected to 40 cm Tygon tubing was filled with the cell suspension. The metal adapters with connected tubing and syringes were gently pushed into the inlets of the side channels on the PDMS chip, the syringes were mounted on a double syringe pump (NE-4000, NewEra, Pump System Inc.) and flow was started at 0.25 ml min 21 . After 20 min, 7 ml of cell suspension was pipetted in the inlet for the middle (observation) channel and the metal adapter for the tubing and 100 ml syringe with the further cell suspension was gently inserted to minimize flow disturbance. The flow in the side channels was maintained at 0.25 ml min 21 , whereas the observation channel was operated at 3 nl min 21 using a separate syringe pump (NE-1000, NewEra, Pump System Inc.). The glass slide with PDMS chip was mounted on an inverted Leica DMI 4000 microscope, and focus was maintained on the first millimeter of the observation cell channel (between positions 1500-2500 mm in Fig. 3A) using an N PLAN 103 objective. Phase-contrast (10 ms) and fluorescence (gfp, 110 ms and mcherry, 200 ms) was digitally imaged every 5-10 min during 1 h using a black-and-white DFC 350 FX R2 camera (Leica). A select number of experiments was complemented by video series of images taken every 5-10 min during 13 s in order to record cell trajectories. Images were exported as 16-bit TIF and analyzed for (fluorescence) intensity in both channels using ImageJ. A zone of 100 3 600 pixels ($101 3 610 mm) was drawn at a distance of 400 microns from the beginning of the filter (approximately the transect at 2200 mm in Fig. 1C, Supporting Information Fig. S3). Fluorescence profiles were extracted using ImageJ by averaging the grey values across a 100-pixels-width segment spanning the complete channel transect (600 mm, $600 pixels). Grey values of ten 100 pixel-lines were then again averaged (obtaining 60 values, each of a 100 3 10-pixel zone).
The fluorescence values were finally normalized by the total fluorescence in the zone. For the CI (chemotaxis index) the ratio of fluorescence was calculated between the 100 pixels ($100 microns) segment closest to the source compared to the total transect fluorescence. The CAR (cell accumulation rate) was then quantified as the slope of CI over time (from at least three time points).
Individual cell trajectories were extracted from three zones of cropped 200 3 200 pixels near the sink and the source channels, by using the mosaic 2D tracking plugin in ImageJ (Sbalzarini and Koumoutsakos, 2005). Coordinates of the trajectories were exported to MATLAB (version R2015_b, MathWorks, Inc.) for further analysis. Total trajectory lengths, effective travel distances (start to end coordinates) and average trajectory angles compared to flow (08) were calculated for each individual cell. The distributions of those three parameters were then used to separate cells into 'swimmer' and 'nonswimmer' categories, with nonswimmers characterized by a ratio of trajectory length divided by effective travel distance of 1 (610%) and a trajectory angle of 086108. Next, the ratios of swimmers to nonswimmers were calculated per attractant and buffer (sink) zone, averaged across three pixel areas per zone and statistically compared with a Welch two sample t-test.

Steady-state model description
Geometry. The microfluidic device had three parallel channels (domains X 1 , source; X 2 , sink and X 3 , cells), connected by two filter zones (domains X 4 and X 5 ; Supporting Information Fig. S5A). For a reasonable balance between computational burden and model accuracy, we chose to represent the device in a two-dimensional space (i.e., top-view). However, a shallow-channel approximation was implemented to correct for the large friction forces induced by the parallel walls on the fluid flow. All the geometry dimensions and physical model parameters are listed in Supporting Information Table S1.
Fluid flow. A correct representation of fluid flow (hydrodynamics) through the microfluidic device is essential for the further description of solute (attractant) transport and microbial chemotaxis. Due to the very low Reynolds number (Re < 0.02), the inertial term in the momentum balance was neglected. The creeping flow is, therefore, described by the Stokes equations for stationary, laminar and incompressible flow (Bruus, 2008): where the variables are p (liquid pressure) and u (velocity vector), with a constant dynamic viscosity l of water at 208C. In microfluidic devices, the hydraulic resistance induced by the top and bottom parallel walls of the flow channels is very important, leading to a large contribution in the pressure drop. To correct the two-dimensional momentum balance for this shallow-channel effect, an additional body force was introduced in Eq. (1) where h is the channel height: 14 lm for the main flow channels (X 1 , X 2 and X 3 ) and 0.25 lm for the filters (X 4 and X 5 ). The source and sink channel inlets C i,1 and C i,2 were defined as laminar flows with the same mass flow rate F in,l . Another laminar flow with flow rate F in,m was set at the observation channel inlet C i,3 . The outflows C o,1 , C o,2 and C o,3 were all assigned a zero gauge pressure, being open to the air. All the other external boundaries were no-slip walls (u 5 0), while flow continuity existed between all flow domains.
Solute transport. The spatial distribution of chemoattractant (serine) or fluorescent tracer molecule (rhodamine B) concentrations c S in all domains resulted from the material balances including convection and diffusion terms: where D S is a constant diffusion coefficient of the solute. If the chemoattractant substance was taken up by the microbial cells, an additional reaction term R S entered the solute balance (4) only on the domain X 3 (in the other domains R S 5 0): The chemoattractant reaction kinetics was described by a saturation rate expression, with v max the maximum specific uptake rate, K m the Michelis-Menten coefficient and c X the local concentration of cells.
In the source channel inflow C i,1 a constant concentration c S 5 c S,i was applied, while in the other two inflows C i,2 Chemotaxis-microfluidic biosensor 253 and C i,3 there was no solute (c S 5 0). All outflows C o,1 , C o,2 and C o,3 were purely convective (no-diffusion boundary type). All the other external boundaries were insulated (noflux conditions), while solute flux continuity was assumed between all domains (e.g., Fig. 1C).
Cell transport. The spatial distribution of microbial cells in domain X 3 was represented by a concentration c X . The concentration of cells results from a material balance with an additional chemotactic term besides the convection and isotropic motility, according to the Alt's cell balance modifications of the Patlak-Keller-Segel model of chemotaxis (Keller and Segel, 1971;Chen et al., 1998): D X is a constant motility coefficient for the cells swimming against the gradient of cell concentration rc X (similar to isotropic molecular diffusion). The chemotactic sensitivity function models the oriented motility of cells along the gradient of substrate rc S , proportional with a chemotactic coefficient D ch and the local concentration of cells c X . The chemotactic coefficient was then assumed to be inversely dependent on the attractant concentration c S according to (Chen et al., 1998): with v 0 the chemotactic sensitivity coefficient and m a scaling factor [typical values for E. coli are v 0 5 10 25 -10 24 cm 2 s 21 , and m53 (Chen et al., 1998)]. K C is the receptorligand dissociation constant [typical value for the Tsr receptor and its ligand serine is 30 mM (Kalinin et al., 2009)]. We included an upper limit for the chemotactic coefficient, motivated by a maximum reachable density of cells c X,max so that D ch 5 D ch;0 12c X =c X ;max À Á . A constant concentration c X,i was fed through inflow C i,3 . The outflow C o,3 was purely convective, thus, the normal components of the diffusion and chemotaxis fluxes were both set to zero. All the other boundaries of domain X 3 (walls and filter/channel edges) were insulated so that the total flux of cells through these was zero.
Model solution. The steady state model equations were solved by a finite element method, implemented in COM-SOL Multiphysics (v5.2, COMSOL Inc., Burlington, MA). To improve the numerical convergence, the equations were solved in a sequential approach, first fluid flow, then solute transport and finally the cell transport. In the case of chemoattractant uptake, the equations for solute and cell transport with chemotaxis and reaction were solved simultaneously. A triangular finite element mesh of max. size 50 lm was created for domains X 1 , X 2 and X 3 , while a rectangular mesh with elements of 1 lm by 10 lm was used for each filter region X 4 and X 5 . In the cell channel near first filters the mesh was refined to max. 20 lm and min. 1 lm sizes in order to capture the very strong cell concentration gradients (Supporting Information Fig. S5B).

Kinetic model description
Metropolis algorithm. In order to describe the dynamics of bacterial chemotaxis and cell accumulation over time in the microfluidic chip, we constructed an individual cellbased model based on the so-called Metropolis algorithm (Glazier et al., 2007;Tu, 2013).
Simulation in absence of chemoattractant. In absence of chemoattractant the cells have no preferred direction of motion, and their path can be described by a random walk. At the start of simulation, n cells (n 5 5000) are placed at random positions on a grid corresponding to a width of 600 mm and a length of 500 mm of the microfluidic channel (simulating positions between 2000 and 2500 mm in Fig.  3A). At each time step, for every cell i at point p i 5 (x i ,y i ) a new position p next with d corresponding to the distance between p i and p next i , traveled within one step of the simulation and u a random angle uniformly chosen between 0 and 3608. This defines a continuous space Markov chain. In the following, the transition probability associated with a move from location p to p' will be denoted by q 0 (p,p'). Cells tending to move beyond the y-boundaries are stopped there (i.e., stick to the channel walls). No boundaries are put in place for the x-direction, as cell accumulation over time is only counted within a 100 mm window between coordinates 2200 and 2300 mm that is encompassed by the larger simulation window (2000-2500 mm).
Simulation in presence of chemoattractant. In presence of chemoattractant the cells will bias their movement direction according to a probability l b to reach position p 5 (x,y) on the grid, given by the Gibbs-Boltzmann probability measure: where H is a general function (described below), b a free parameter fitted from the experimental data shown in Fig.  2 and Z b ð Þ the partition function of the system, defined by: and with K being the set of all possible states for p. Originally, the problem was to compute average values of physical observations under the Gibbs-Boltzmann probability measure. However, in most frameworks, the partition function cannot be computed theoretically or numerically, so that it is impossible to sample this distribution by direct simulation. Therefore, we designed a Markov chain with computable transition probabilities having the Gibbs-Boltzmann distribution as a steady-state. In our framework, we consider the exploratory transition kernel q 0 p; p 0 ð Þ, the probability to move from site p to p 0 in absence of attractant. The decision in the simulation to go from site p to p' is obtained from the variation of free energy DH5H p 0 ð Þ 2H p ð Þ. H p ð Þ contains the main relevant information on site p, such as the local chemoattractant concentration. For the Metropolis chain the probability q b p; p 0 ð Þto go from position p to position p 0 is described as: where DH ð Þ pos 5 max 0; DH f g. The function H represents the bacterial sensing of the chemical gradient, and is dependent on the relative change of attractant at any position and the chemotactic sensitivity and receptor saturation of the bacterial cells, similar to Eqs. (7) and (8). In a first representation H*, the function is: with c S (p) representing the ligand concentration at position p5(x,y), K C being the receptor-ligand dissociation constant and n 0 being a unitless factor to account for chemotactic sensitivity (similar in scale as the v 0 used in the continuum model). Used parameter values are listed in Supporting Information Table S2. Attractant gradients were taken from the continuum-based COMSOL models at position 2200 mm for each of the used serine concentrations (1 and 10 mM, 0.1 and 1 mM). For simplicity, we applied a c S,0 which is the average across the gradient at the given source concentration, whereas c S (p) is the actual concentration at any position p across the simulated channel width. Equation (13) is further corrected by a factor a for the difference in measured and true chemotactic sensitivity, according to the Stevens' power law (Stevens, 1957). Finally, we multiply H* by a factor a that describes the adaptation of the cells over time to the chemoattractant, and which takes the form of a logistic function: ð Þ with k being the curve steepness, s 0 the value of the sigmoid's midpoint along the s-axis and s the simulation step axis (discrete values). A total number of 1200 simulation steps was carried out, which corresponds in the experiments to 60 min. Putting all those terms together leads to the following equation for the function H that was used in Eq. (10).

Model solution
For each simulation step the next position p next i for every cell is computed using the following sequence: If instead the new position has a higher free energy, the chain will take the new unfavourable position with p b p target i ; p i <1, and keep the old one otherwise. The rational for keeping a (small) probability of chosing a bad position (i.e., with higher free energy) is that it might help in escaping local free energy minima. Kinetic simulations were scripted and run in Julia (Supporting Information Methods Metropolis script) (Bezanson et al., 2017). The cell distribution at fixed simulation steps within a 100-mm width zone between coordinates 2200 and 2300 mm (as in Fig. 3A) was extracted from the total grid (600 3 500 mm 2 ) and is presented in Fig. 5A. The chemotaxis index (i.e., the proportional accumulation of cells in the first 100-mm across the channel) was calculated from fivefold repeated simulations for every serine concentration and for absence of serine.

Supporting information
Additional Supporting Information may be found in the online version of this article at the publisher's web-site:  Fig. 1C. (E) Modelled effect of serine metabolism of E. coli on the established gradient at 10 mM serine source concentration and position 2200 mm in the observation channel. Note that for ease of comparison to the chip design, the distance across the channel is represented from 300 to 2300 mm. Fig. S2. Nonchemotactic cell distribution of E. coli DfliCmcherry toward serine. Images showing the distribution of DfliC-mcherry at the 2100-2200 mm location on the chip over time (0 -40 min) as a function of the indicated serine concentration. Top is sink channel, bottom is source channel. Fig. S3. Chemotaxis index measurement setup. Chemotaxis response was quantified in a zone of 600 3 100 microns at a distance of 400 microns from the beginning of the filters. Fluorescence intensity profiles were extracted from the fluorescence images using ImageJ and normalized by the total fluorescence in the zone of measurement. The chemotaxis index was calculated as the proportion of fluorescence in the 100 mm segment closest to the source of attractant compared to the total fluorescence across the channel. Fig. S4. Microfluidic chip fabrication procedure. The fabrication procedure starts with a silicon wafer (1). A photolithography process produces a layer of resist at the filter position that protects this zone during the etching step (2). The etching results in the formation of the negative of the 650 nm high channels of the filters (3). A second step of photolithography produces the mold of the channels with a resist layer of 14 microns high (4). This inverted mold is used multiple times to produce the PDMS chips, by pouring PDMS on it and let polymerize (5). Once polymerized, the PDMS is peeled off the inverted mold and, after punching holes for the inlets, is bonded to the glass slide by a plasma treatment (6). Fig. S5. (A) Model geometry, dimensions, domains and boundaries. X 1 : Source domain (fed with chemoattractant solution), X 2 : Sink domain (fed with water), X 3 : Cells domain (fed with a suspension of cells), X 4 and X 5 : Filter domains (separate the cells from source and sink channels). C i,1 , C i,2 , C i,3 : Inflows, C o,1 , C o,2 , C o,3 : Outflows. The geometry dimensions are listed in Table S1. (B) Finite element mesh detail in the neighborhood of the filter region. Table S1. Parameters of the continuum steady-state model. Table S2. Parameters, variables and functions used in the Metropolis model.