Logarithmic sensing in Bacillus subtilis aerotaxis

Aerotaxis, the directed migration along oxygen gradients, allows many microorganisms to locate favorable oxygen concentrations. Despite oxygen’s fundamental role for life, even key aspects of aerotaxis remain poorly understood. In Bacillus subtilis, for example, there is conflicting evidence of whether migration occurs to the maximal oxygen concentration available or to an optimal intermediate one, and how aerotaxis can be maintained over a broad range of conditions. Using precisely controlled oxygen gradients in a microfluidic device, spanning the full spectrum of conditions from quasi-anoxic to oxic (60 n mol/l–1 m mol/l), we resolved B. subtilis’ ‘oxygen preference conundrum’ by demonstrating consistent migration towards maximum oxygen concentrations (‘monotonic aerotaxis’). Surprisingly, the strength of aerotaxis was largely unchanged over three decades in oxygen concentration (131 n mol/l–196 μ mol/l). We discovered that in this range B. subtilis responds to the logarithm of the oxygen concentration gradient, a rescaling strategy called ‘log-sensing’ that affords organisms high sensitivity over a wide range of conditions. In these experiments, high-throughput single-cell imaging yielded the best signal-to-noise ratio of any microbial taxis study to date, enabling the robust identification of the first mathematical model for aerotaxis among a broad class of alternative models. The model passed the stringent test of predicting the transient aerotactic response despite being developed on steady-state data, and quantitatively captures both monotonic aerotaxis and log-sensing. Taken together, these results shed new light on the oxygen-seeking capabilities of B. subtilis and provide a blueprint for the quantitative investigation of the many other forms of microbial taxis.

in COMSOL Multiphysics 4.4 (see Materials and Methods). Oxygen diffusion dynamics in the test channel were simulated for two gradients: 0%-20% and 0%-10% oxygen (dashed lines in Fig.   S1). We then set out to quantify how the spatial profile of oxygen varied as a function of time for both gradients.
To measure oxygen concentrations in the test channel we flowed in the test channel a 167 ppm solution of ruthenium tris(2,2'-dypiridyl) dichloride hexahydrate (RTDP) in 66% ethanol in water, at a flow rate of 200 nL/min. RTDP is a fluorescent dye sensitive to oxygen: the larger the oxygen concentration, the smaller the intensity of the fluorescence that RTDP emits. Consistently with previous studies 16 we used the Stern-Volmer equation I 0 /I = 1 + K q [O 2 ] to convert the fluorescence intensity I in an oxygen concentration [O 2 ]. First we need to estimate I 0 , the fluorescence intensity in absence of oxygen (100% nitrogen) and the quenching constant K q . To do so we flowed pure nitrogen (0% oxygen) in the source and sink channels, waited 10 minutes to make sure the gas concentration in the channel was equilibrated to uniform, and then acquired a fluorescence image of the channel. Background estimation and correction was carried out as in 16 ; to this aim we extracted background fluorescence in the test channel fitting a second order polynomial across the x axis (i.e. the direction of the gradient) to intensities of areas 100 µm in the left and right PDMS walls -as there is no dye in the PDMS, and PDMS is not autofluorescent at the RTDP emission wavelength, we reasoned that any fluorescence in these areas can be classified as background.
This procedure yielded an estimate of background fluorescence in the test channel -obtained using the fitted polynomial-that we used for correction by subtraction to all the intensity profiles we acquired 16 . As commonly noted, the quality of the micrographs decreased quickly in the vicinity of PDMS walls; as this made a reliable measurements of signals very close to the boundaries of the test channel challenging, we decided to analyse oxygen concentrations between 10 and 450 µm. In the same manner we measured a second reference intensity, I air , by flowing air (20.8% oxygen) in both the source and sink channels. This allowed us to calculate the quenching constant To assess the accuracy of our mathematical model in predicting the spatiotemporal profile of oxygen, we generated (in two separate experiments) the two gradients simulated with our model, namely 0%-10% and 0%-20% oxygen. For each case, we quantified the background fluorescence, flowed in the source and sink the gas mixes appropriate to generate the desired gradient (e.g. nitrogen in the sink and 20% oxygen in the source for 0%-20%) and acquired fluorescence images every 10 seconds for 5 minutes. We then converted the fluorescence intensity values into oxygen concentrations with the procedure described above. The results of this approach are presented in We note that, in the device used for our experiments, if we denote by O source and O sink the concentration (expressed in %) of oxygen flown in the source and sink channels, respectively, cells are exposed to >90% of the gradient from O sink to O source , and <10% of the gradient occurs within the lateral PDMS boundaries separating the source and sink channels from the test channel. This can be easily observed in Table 1. When O source =100% and O sink =0%, the boundary conditions in the test channel are C(0 µm)= 6.04E-5 M, i.e. 4.7% of 1.3E-3 M (oxygen saturation in water in the lab), and C(460 µm)= 1.24E-3 M, i.e. 95.4% of 1.3E-3 M. This corresponds to a total drop in oxygen concentration within the test channel of ∼90.7%, to be compared to a 100% drop between the source and sink channels. This also means that ∼ 9.3% of the gradient is retained in the PDMS walls and is not available to the cells.
Bacterial diffusivity D B In order to measure the diffusivity of B. subtilis we tracked and analyzed bacterial trajectories in uniform concentrations of oxygen ranging from 0% to 100% (Fig. S2).The (2D projection) Mean Squared Displacement (MSD) of cell, subjected rotational can be written as: where V is cell's swimming speed, t is time and τ R is the characteristic time-scale associated to rotational diffusion. We measure V directly (Fig. S2C) from bacterial trajectories and obtain τ t and, therefore τ R via fitting 33 ,32 . In agreement with what has been reported in literature 34 we measure a tumbling time τ t = 1 2 τ R 0.71s at low oxygen concentrations (O 2 <1%) and higher tumbling times τ t 1.18s for O 2 > 1% (Fig. S2B). Consistently with previous reports our data also suggest the swimming speed increases with the concentration of oxygen (see Fig. S2C) up to ∼1% O 2 . We can use these observations to derive the translational diffusion coefficient: Mechanistic derivation of advection-diffusion equation In this section, we will show how an advection-diffusion equation for densities, of the type that we fit to data, might be reasonable. As little is known about the mechanistic basis of B. subtilis aerotaxis 35 our approach is as follows.
We will first review an accepted and experimentally validated model of E. coli, and show how it leads to an advection-diffusion equation of the desired form. We will then see how this mechanism would be modified by incorporating knowledge about the differences between E. coli and B.
subtilis chemotaxis, and we will show that the same advection-diffusion equation results in spite of this difference (albeit with very different parameters). As aerotaxis and chemotaxis in B. subtilis employs the same receptor mechanism [11], we will postulate that this same model applies to aerotaxis.
We organize this section by first discussing a general approach to advection-diffusion approximations, before specializing to the E. coli and B. subtilis models.
Preliminaries Let p(x, y, ν, t) be a density function describing a population of "particles" or agents (for example, bacteria), modeled in a 2N + m dimensional phase space, where at time t, x = (x 1 , . . . , x N ) ∈ R N (N = 1, 2, 3; we soon specialize to N = 1) denotes the position of the agent, y = (y 1 , . . . , y m ) ∈ Y ⊂ R m ≥0 denotes the internal states of the agent (we will soon specialize to m = 1), and ν ∈ V ⊂ R N denotes its velocity. Also, S(x) = (S 1 , . . . , S M ) ∈ R M denotes the concentration of signals in the environment which are sensed by each agent at space location x (we will soon specialize to M = 1). The external signal S is assumed to be constant in time (steady state assumption on chemoattractant), but is allowed to depend on space coordinates.
We assume that the following system of ordinary differential equations describes the evolution of the intracellular state, in the presence of the extracellular signal S(x) at the current location of the agent: where f : S with respect to space (local gradient of chemoattractant). In most models, f depends explicitly only on y and S, but we allow this additional generality in the theory.
We assume also given an instantaneous reorientation ("tumbling") rate λ = λ(y, S(x), S (x)) (often, λ depends only on certain combinations of y and S(x), represented by the "activity" of re-ceptors), the evolution of p is governed by the following transport (or "Fokker-Planck" or "forward Kolmogorov") equation 36 (omitting arguments of functions p and f , for readability): where the nonnegative kernel T (y, ν, ν ) is the probability that the agent changes the velocity from ν to ν if a change of direction occurs. Also V T (y, ν, ν ) dν = 1.
The main goal here is to derive an approximate macroscopic model for chemotaxis using the microscopic model (4), i.e., we want to find an equation to approximately describe the evolution of the marginal density: by adapting methods from Grunbaum [24] and Othmer [25]. We will assume that the external signal is isotropic in two state directions, so that in effect we can study one-dimensional motion.
A general equation in one dimension From now on, we study the movement of agents in one dimension have constant speed, so that the velocities are ν ∈ {ν, −ν}, whereν is a positive number, which we'll think of as a parameter in the equations. We will write f + (y,ν, S, S ) instead of f (y,ν, S, S ) and f − (y,ν, S, S ) instead of f (y, −ν, S, S ), and omit the bars fromν from now on. Similarly, for p, we let p ± (x, y, t) denote the density of particles that at time t, are located at position x, with the internal state y, and with the constant speed ν, and moving to the right (+) or left (−) respectively.
The internal state evolves according to the following ODE system: where f ± : R ≥0 × R × R × R → R are continuously differentiable functions in each argument that describe the evolution of internal state of agents which move to the right (+) and left (−) respectively.
Note that we are allowing f to depend on the direction of movement as well as ν and S , the derivative of S with respect to space.
In our examples, f + = f − only depends on y and S, but we can consider the more general dependence in these preliminary derivations.
We describe the tumbling rate by introducing: λ(y, S, S ) = g(y, S, S ), where g is a continuous function.
Then, according to Equation (4), p ± (x, y, t) satisfy the following coupled first-order partial differ-ential equations: See [25] for existence and uniqueness of solutions of (8)-(9) We assume given a forward-invariant set I ⊂ R ≥0 , i.e., if y(0) ∈ I, then y(t) ∈ I, for all t ≥ 0, with the property that p ± 0 (x, y) are supported on I, i.e., p ± 0 (x, y) = 0, when y / ∈ I. (In each of the examples to be considered below, such a set I will be constructed, by appealing to Lemma 1 in Section below). In other words, p ± (x, y, t) = 0, ∀x, y / ∈ I, t ≥ 0.
The objective is to derive an approximate equation for the macroscopic density function using the microscopic model (8)-(9), by adapting a technique from [25]. To this end we introduce a flux variable j as well as moments associated to n and j: Note that by Equation (10) all the moments are well defined.
Next, we assume f + = f 0 + νf 1 , and f − = f 0 − νf 1 , where the Taylor expansions of f 0 and f 1 , with respect to the internal state y, are given as follows: for some A i 's and B i 's that are functions of S, S , and ν 2 . (We formally assume that these expansions exist.) Also we consider the following Taylor expansion for g(y, S, S ): where the a i 's are functions of S, and S .
In addition, we assume A 0 = 0, because this is satisfied in our examples. Then by multiplying (8) and (9) by 1, ν, and/or y, adding or subtracting, and integrating with respect to y on R ≥0 , and applying the fundamental theorem of calculus and integration by parts, we obtain the following equations for macroscopic density and flux and their first moments: Note that by Equation (10), p ± = 0 outside the interval I, therefore, for any i = 0, 1, . . .

Parabolic scaling
In this section, we introduce a parabolic scaling to derive an approximate chemotaxis equation from the moment equations (16)-(19). Let L, T , ν 0 , y 0 , and N 0 be scale factors for the length, time, velocity, internal state, and particle density respectively, and define the following dimensionless parameters (we use hats to denote the dimensionless forms of the parameters): The parabolic scales of space and time are given by: for any arbitrary . Therefore, the dimensionless form of moment equations (16)-(19), for = T ν 0 L , become: Next, we write Equations (25)-(28) in a matrix form, as follows: whereŵ = n,ĵ,n 1 ,ĵ 1 T and the matrices P , Q, and R defined as follows: Assuming the regular perturbation expansion for w, and comparing the terms of equal order in in (29), we get: From the last equality of Equation (31), we can derive the following equation forĵ 1 1 : By substitutingĵ 1 1 into the second equality of Equation (31), we obtain the following equation Now we compare the terms with order 2 : Note that (1, 0, 0, 0) T is in the kernel of R and the right hand side of (33) is in the image of R.
Therefore their inner product is zero: Equation (32) together with Equation (34) give the following equation for n 0 in the dimensionless variables: Since (35) leads to the following chemotaxis equation in dimensionless variables: Changing back to the original (dimensional) variables, we obtain the following PDE: Examples

E.coli
The following simplified one-dimensional model provides a phenomenologically accurate model of the chemotactic response of E.coli bacteria to MeAsp; see for example 39 , 37 . The internal state evolves according to an ordinary differential equation: which describes the methylation state of receptors, where a is a number between 0 and 1 that quantifies the fraction of active receptors, and is written as follows: in terms of free energy differences due to methylation and ligand respectively: where K I and K A are dissociation constants for inactive and active Tar receptors, respectively.
This arises from an MWC 38 model of clusters of N receptors that rapidly switch between active and inactive states, In summary, we write: and K, K I , and K A are nonnegative constants and K I < K A .
With appropriate parameter choices 39 , 37 , this model fits very well the response of E. coli to the ligand α-methylaspartate.
E. coli tumbling rate is controlled by the concentration of cheY-P. In this simplified model, one thinks of phosphorylation state of cheY as directly proportional to activity, assuming fast phosphotransfer. Thus, one takes the jump (or "tumbling" for bacteria) rate in the form: Here a 0 denotes a steady-state kinase activity, H a motor amplification coefficient, and τ the average run time. We write It is convenient to use y = e αm as a state variable, instead of the methylation level m. So the equations can be rewritten as follows: provided that we pick Observe that F m = e α /y when expressed in terms of the new variable y. The parameters p, q, K, N , and H are all positive, and, from its definition, it is clear that q is between 0 and 1.
The objective is to derive a parabolic equation for the macroscopic density function. It is convenient to define a new internal state variable as follows: Then, a simple calculation shows that and For convenience of notation, let us define G(S) := log See 57 for a proof.
Let L, T , ν 0 , and N 0 be scale factors for the length, time, velocity, and particle density respectively, and define the following dimensionless quantities: A simple calculation shows that: All other parameters remain the same as in Equations (20)-(22), and Equation (24), for y 0 = 1 T .
Note that for = ν 0 T L , we have the following analogous result to Lemma 1, in hyperbolic scale: Definition 1 (shallow condition). If G (S) G (S) ≤K, whereK = O(1), we say S has a shallow gradient.
Lemma 2. Assume that i.e., S has a shallow gradient. Then, for any i ≥ 1, See 57 for a proof.
Remark 1. Equation (45) is equivalent to the following condition for G (S): or equivalently ν c Note that for exponential signal S(x) = e ρx , using condition (47), when ρ is small enough, we are in a shallow gradient regime. For linear signal S(x) = ax + b, using condition (47), when a is small enough, we are in a shallow gradient regime.
Using the notations of Equations (13)- (14), In order to derive an advection-diffusion approximation using Equation (37), we just need to find the first two terms of the Taylor expansion of λ(w) in (42). We do that next.
A simple calculation shows that where Q(w) is the sum of higher orders of w in the Taylor expansion. Plugging the new values of a 0 and a 1 into Equation (37), we get the following advection diffusion equation: where

Modifications for B. subtilis
It is known that the activity of B. subtilis chemotatic receptors increases in the presence of attractants. This means, in effect, that the roles of K I and K A are inverted in the formula for activity: now K I > K A . Furthermore, tumbling (due to CW rotation of flagella) is induced by lack of activity, which we may model by replacing a by the fraction of inactive receptors, 1 − a, in the simplified E. coli model considered earlier.
Thus, we now assume that the internal state evolves according to the following ODE system: where we use the following form for activity: a = 1 1 + K S + K I (S + K A ) y N and p, q, K, and N , K I , and K A are positive constants, where now K I > K A . Recall that q is between zero and one.
We assume now the following form for the tumbling rate: where A and R are positive constants. (We assume that A > q, which is the case if A = 1.) The objective is to derive a parabolic equation for the macroscopic density function.
As in the previous example, let w = p(a − q). Then, a simple calculation shows that Since dw dt is exactly the same as in Example , we get the same expressions for A i 's and B i 's, namely: In order to derive an advection-diffusion approximation using Equation (37), we just need to find the first two terms of the Taylor expansion of λ(w). We do that next.
A simple calculation shows that where Q(w) is the sum of higher orders of w in the Taylor expansion. Plugging the new values of a 0 and a 1 into Equation (37), we get the following advection diffusion equation: where that can be also rearranged in a more compact form, gives us Equation (2) as presented in the main text: with Thus, a formula of exactly the same form as for E. coli has been obtained.
Our mathematical model best captures aerotaxis in B. subtilis In order to assess how the model we propose compares to alternative solutions in literature we grouped previous advection-diffusion chemotaxis models in 3 main classes: KS, LS and RTBL models (see following section). Each of these models has a different expression of the chemotactic speed V C and they range from fully phenomenological (e.g. KS) to biophysically-informed approaches (like RTBL). The vast majority of the other advection-diffusion models used to capture chemotaxis can be derived from the ones we consider in the following.
We compared the performance of the models by plotting the prediction (Fig. S3-S9) of the best combination of parameters the optimization algorithm found over 100 iterations and its prediction error (Fig. 5, see Eq. 5 in the main text). For each model we also plot the distribution of prediction errors of the 100 solutions to the optimization problem.
Notably, for the KS model the genetic algorithm consistently identified a single solution to the parameter optimization problem (Fig. S3), hence the tight distribution in Fig. 5. Similar results in terms of prediction accuracy (and therefore SSE, see Fig. 5) can be achieved using the best solution identified for the LS model (Fig. S4). A significant improvement, instead, can be achieved using the RTBL model ( Fig. S5 and Fig. 5): the best parameter set found in this case achieves an SSE significantly smaller than in the previous cases (0.95·10 −1 compared to 1.84 ·10 −1 for the LS and 1.90 ·10 −1 for the KS models). However the model we propose displays the smallest prediction error (0.73 ·10 −1 , Fig. 5) and therefore best captures the body of experimental data we describe (Fig. 1C).

SI Materials and Methods
Growth protocol

Microfluidic fabrication, experimental operation and image analysis
In order to generate oxygen gradients, the source and sink channels were each connected to a gas-mixing unit, supplied by gas tanks (Air Gas, MA). We used 100% nitrogen as well as 0.1%, 1%, 20% and 100% oxygen/nitrogen mixtures. Each gas-mixing unit was composed of two high- in the source and sink channels and allowed them to diffuse within the device. Of note, the presence of the 220 µm thick PDMS wall separating the test channel from the sink channel implied that the minimum oxygen concentration in the test channel was higher than the concentration in the sink channel. Similarly, the maximum concentration in the test channel was lower than the concentration in the source channel. For example, a 0%-100% case (0 M in the sink channel and ≈ 8 mM, on the other end, at the interface between PDMS and the source channel) corresponds to an oxygen gradient ranging from 4.6% (60 µM) to 95% (1.24 mM, 100% oxygen in water corresponding to 1.3 mM) in the test channel (see Table 1 in the Supplementary Information). Bacteria were then injected in the test channel and glass coverslips were used to seal its inlet and outlet of the test channel to suppress any residual flow. Cells reached steady state distribution within 5 minutes after the injection (Fig. 4). We then used an automated acquisition routine to capture 30,000 phase-contrast images of the same location along the test channel (equidistant from the inlet and outlet) at 67 ms intervals over 33 min (20 objective; Andor Zyla camera with 6.5 µm/pixel (leading In order to fully characterize the model we need to identify each of its three parameters K 1 , K 2 and χ 0 -note that D B is measured experimentally (see Supplementary Information and Fig. S2). To this aim we developed a genetic-algorithm-based multi-experimental fitting procedure designed to find the combination of parameter values that minimized the sum of the squared errors between model predictions and experimental data where n = 33 is the number of experimental designs, w(x) is a vector of weights increasing linearly from 1 to 1000 (empirically found to ensure the best results in terms of prediction error were attained), B E (x) are the experimental data and B S (x) the simulated accumulation profiles via numerical integration (∆x = 10 nm) of: with test channel width W = 460 µm and f (ξ) = 1/((K1 + C(ξ))(K2 + C(ξ))) for the model in Eq. 2 in the main text. This expression of B(x) can be obtained plugging Eq. 2 in Eq. 1 in the main text, using the linearity of the oxygen gradient (i.e. ∇C independent of ξ) and posing ∂B ∂t = 0. At each iteration the genetic algorithm generated a number of random solutions, ranked them based on Eq. 55, the worst solutions, selected the best ones and applied "cross-over" and "mutation" to obtain new solutions to be evaluated at the next iteration 50 . The search for a solution stopped when a stall was detected, i.e., when the average change in SSE(χ 0 , K 1 , K 2 ) over 50 iterations was smaller than 10 −6 . The reported parameter set is the best combination identified over 100 repetitions of this procedure. We adopted the same method to identify the parameter values for all models (see Supplementary Information).

Robustness analysis of parameter estimates
Although very powerful at solving complex optimization problems, Genetic Algorithms do not provide any guarantee of convergence. As a consequence of this, a set of "optimal parameters" obtained as a result of the optimization, might actually be a local, rather than a global solution -these are solutions that optimize the objective function in a sufficiently large neighborhood of, but not the entire, space of parameters. Yet, at the end of the parameter optimization process we would ideally identify a set of values that minimizes the cost function (Eq. S62) globally rather than locally.
To assess whether the values obtained from the Genetic Algorithm could be outcompeted by other combinations of values, we decided to adopt a Naive Grid Search approach. The principle behind this method is simple: the set of values each parameter can take is discretized and all the combination of discretized parameters are evaluated using the cost function. The more fine-grained the discretization is, the more this approach resembles an exhaustive search. The main limitation of this approach is that for large numbers of parameters and/or parameter values the number of objective function evaluations quickly increases and ultimately makes the problem intractable.
As customary in these cases, we assigned to each parameter identification task (i.e., each model among the ones we considered) a budget of "function evaluations" equal to 10 5 . For each of the i parameters in that model, we identified a physically feasible set of values, and discretized it into M values, with M being the closest integer to 10 5 i . We then evaluated the cost function for each of these combinations and, for each model, the value of the minimum cost identified by the Naive Grid Search method was compared to the minimum found by the Genetic Algorithm (Fig. S10).
For both the KS and the LS models (1 and 2 parameters, respectively) we confirmed that the Naive Grid Search identified values of the optimal parameters substantially undistinguishable from the ones returned by the Genetic Algorithm. For the RTBL and the Finite Range Log-sensing regime, instead, the Naive Grid Search algorithm returned values different from the Genetic Algorithm and, in both cases, characterized by higher value the cost function -suggesting that the parameter values identified by the Naive Grid Search are not global optima. These results indicate that it is unlikely that the parameter sets identified by the Genetic Algorithm for our model represent local optima and that they are instead the global optima we sought.

Model validation on transient aerotaxis
As a stringent validation of the model, we tested its performance in predicting the population migration in a transient aerotaxis experiment. At the start of the experiment, sink and source channels both contained a flow of 21% oxygen and cells were allowed time to equilibrate to their steady state distribution, which was uniform given the uniform oxygen concentration (Fig. 4). At time zero we started flowing 0% and 0.05% oxygen in the sink and source channels, respectively, and recorded the spatial distribution of bacteria across the test channel at 100 frames/s for 4 min To produce To test the ability of this model to recapitulate our experimental results, and compare its prediction capabilities with other models, we adopted the same approach reported in the main text for the "finite-regime log-sensing" model we propose (see "Derivation and Identification of the Mathematical Model"). We ran 100 instances of a genetic algorithm meant to identify the value of χ 0 (the only free parameter in this model) that minimizes the average mismatch between model prediction and experimental results over the whole dataset. It should be noted that, based on Eq. (4) (see main text) and assuming D B does not depend on space, the steady state distribution of bacteria B(x) can be rewritten as: where we observe that χ They succeeded in this effort and proposed a formulation of V C directly proportional to the chemo-tactic sensitivity coefficient χ 0 and inversely proportional to the squared sum of K and the chemoattractant concentration C: By using population scale measurements of bacterial fluxes, not only were Lapidus and Schilller able to identify the values of χ 0 and K, they also showed the predictions of their model were in good agreement with the experimental results.
It is worth noting that, while achieving good performance in capturing the experimental results in 51 , the LS model does not support logarithmic sensing. Moreover, as our understanding of the cascade of signaling events leading to chemotaxis furthered, an increasing number of approaches focused on bridging single cell behavior and population level phenomena.
To assess the ability of this model to capture our data we followed the same approach described for the KS model. In this case, however, the parameters to be identified are both χ 0 /D B and K.
We set non-negativity constraints for this identification task following the same line of reasoning mentioned above and recorded the results of the 100 optimization procedures. The solutions to the optimization problem is plotted in Fig. S4 (χ 0 /D B = 57.89 and K = 1.39 · 10 −5 M).

RTBL model
The RTBL model, developed by Rivero and colleagues 43 , achieves a macroscopic characterization of bacterial chemotaxis using microscopic variables involved in the chemotactic response of single cells (e.g. receptor occupation and swimming speed). In order to derive their model Rivero and colleagues considered two sub-populations of bacteria (p + and p − ) exposed to a chemoattractant gradient in a 1D domain. Each bacterium can either proceed from left to right or viceversa; this will determine which subpopulation it belongs to. Tumbling makes a bacterium switch from one group to the other; just like we would expect to happen in-vivo, the probability of tumbling depends on the time derivative of the number of bound receptors. In this framework, following the steps reported in Appendix A in 55 , one can derive the expression of the chemotactic speed V C : where V is the swimming speed of bacteria. Consistently with what we previously reported, we probed the ability of the RTBL model to capture our dataset running 100 instances of our optimization procedure. In this case the parameters to be identified were three: χ 0 , K and V . For all of them we set non-negativity constraints, following the considerations we previously discussed; moreover we restricted V , the swimming speed, to not exceed 40 µm/s (we set this constraint according to experimental quantification of bacterial swimming speed we obtained while measuring D B ). Also in this case we collected statistics on the prediction error of the solutions identified during the 100 runs of the optimization procedure (Fig.   5) and we plotted the results from the simulation of the best among the 100 solutions identified by the genetic algorithm in Fig. S5 (χ 0 = 7.10 · 10 −8 m 2 /s, K = 7.01 · 10 −6 M and V = 39.4 · 10 −5 m/s). Table 1: Oxygen concentrations inside the test channel. For each oxygen mixture flown in in the sink and source the actual concentrations within the test channel, as well as the number of replicates, are reported here. In each case the bacteria were exposed to a linear gradient with minimum C(0 µm) and maximum C(460 µm).