A novel approach to chemotaxis: active particles guided by internal clocks
Luis Gómez Nava, Robert Großmann, Marius Hintsche, Carsten Beta, Fernando Peruani
AA novel approach to chemotaxis: active particles guided by internal clocks
Luis G´omez Nava, ∗ Robert Großmann,
1, 2, ∗ Marius Hintsche, Carsten Beta, and Fernando Peruani † Universit´e Cˆote d’Azur, Laboratoire J. A. Dieudonn´e,UMR 7351 CNRS, Parc Valrose, F-06108 Nice Cedex 02, France Institut f¨ur Physik und Astronomie, Universit¨at Potsdam,Karl-Liebknecht-Strasse 24/25, Haus 28, 14476 Potsdam, Germany (Dated: June 12, 2020)Motivated by the observation of non-exponential run-time distributions of bacterial swimmers, wepropose a minimal phenomenological model for taxis of active particles whose motion is controlled byan internal clock. The ticking of the clock depends on an external concentration field, e.g. a chemicalsubstance. We demonstrate that these particles can detect concentration gradients and respond tothem by moving up- or down-gradient depending on the clock design, albeit measurements of thesefields are purely local in space and instantaneous in time. Altogether, our results open a new routein the study of directional navigation, by showing that the use of a clock to control motility actionsrepresents a generic and versatile toolbox to engineer behavioral responses to external cues, such aslight, chemical, or temperature gradients.
In the canonical picture of bacteria with run-and-tumble motility, such as
Escherichia coli or Salmonella [1,2], bacteria display exponentially distributed run-times –despite recent findings that suggest the possibility ofnoise-induced heavy-tailed distributions [3, 4] – and per-form chemotaxis by regulating the associated tumblingfrequency. By measuring the chemical concentrationthrough clustered arrays of membrane receptors [5] andsubsequent signal processing via a complex biochemicalcascade [6] that leads to an effective memory [6–9], thebacterium is able to detect and respond to an externalconcentration gradient. It extends the duration of runs,i.e. decreases the number of tumbles when heading inthe direction of increasing attractant concentration, thusperforming a biased random walk towards the attractantsource. Recently, it was discovered that several bacte-rial species differ from this classical picture. Notably,in the soil bacterium
Pseudomonas putida ( P. putida )that displays a run-and-reverse motility pattern withmultiple run modes [10, 11], the distribution of run-times is non-exponential and exhibits a refractory pe-riod, cf. Fig. 1a and Refs. [12, 13]; an observation thatsuggests that run-times are not controlled by a Poisso-nian mechanism [14]. By adapting the reversal statisticsin response to a chemical gradient,
P. putida is able toperform chemotaxis as indicated in Fig. 1b. However,the non-exponential run-time distribution suggests thatchemotaxis in
P. putida , and potentially also in otherbacterial species, may involve a taxis mechanism thatis fundamentally different from the one reported for
E.coli . Similar non-exponential bell-shaped run-time distri-butions were reported for
Myxococcus xanthus [15] and
Paenibacillus dendritiformis [16], which display run-and-reverse motility similar to
P. putida . Non-exponentialrun-times were observed also for the marine bacterium
V. alginolyticus [17] and, notably, also for the rotation ∗ L.G.N. and R.G. contributed equally to this work. † [email protected] time of the flagellar motors of Escherichia coli [18]. Fur-thermore, it has been reported that surface explorationof
Escherichia coli [19] is not consistent with the canon-ical run-and-tumble picture of bacterial motility, but in-volves multiple motility modes.Motivated by the non-Poissonian run-time statisticsof
P. putida – though not pretending to be a realistic,biochemical description of chemotaxis in
P. putida – wepropose a minimal, phenomenological, generic model fortaxis of active particles that inherently produces non-exponential run-time distributions. Particles move atconstant speed and experience velocity reversals. Thekey element of the model is that the occurrence of reversalevents is controlled by an internal clock. The clock mim-ics the fact that the directional response of a microor-ganism to an external signal or field, such as a chemicalconcentration or temperature gradient, requires to sensethis signal, internalize and process it, presumably involv-ing cascades of biochemical events [18, 20], and to executea behavioral response (e.g. a reversal), after which the mi-croorganism continues sensing, processing and respond-ing to the signal. We assume that this complex cycle canbe reduced to a series of stochastic checkpoints or steps,where some or all of them depend on the external signal.These steps are represented by the “ticks” of a clock. Thetransitions between two consecutive ticks are modeled asPoissonian processes with concentration-dependent tran-sition rates. Ergo, the clock dynamics is, by definition,Markovian: it does not incorporate or presuppose mem-ory in any manner. On the other hand, the distribution ofthe times between two consecutive behavioral responsesis non-exponential as observed in
P. putida experiments,cf. Fig. 1a. Experimental details on cell culturing, thechemotaxis chamber and imaging can be found in theSupplemental Material (SM) and Refs. [12, 21–24].In this study, we demonstrate that the design of theclock controls the long-time motility of the particles:some of the clock designs result in signal-insensitive par-ticles and a variety of other designs lead to particles dis-playing actual taxis. The taxis responses include either a r X i v : . [ phy s i c s . b i o - ph ] J un up-gradient or down-gradient biased motion. In short,we show how a clock can be used to guide active par-ticles subjected to external stimuli, and to obtain non-exponential run-time distributions. Model. – We consider active particles that move at con-stant speed v in a one-dimensional system of size L withreflecting boundary conditions. Particles are exposed toa temporally constant external field c ( x ). The equationof motion of one of these particles is given by dx ( t ) dt = v ( t ) = v s ( t ) , (1)where x ( t ) denotes the position of the particle and s ( t ) ∈{− , } indicates its direction of motion. The variable s ( t ) undergoes stochastic transitions such that s ( t ) →− s ( t ). These reversal events are controlled by a stochas-tic M -tick clock. A reversal occurs every time the clockcompletes a full cycle, i.e. in the transition from M to1, as illustrated in Fig. 2. Thus, the state of a par-ticle is characterized by its position x ( t ), its orienta-tion s ( t ) as well as the internal state (or tick) of theclock n ∈ { , , . . . , M } . Clock dynamics. – The transition between tick n to n + 1 is modeled by a Poisson process with a concentra-tion dependent transition rate γ n ( x ) = f n [ c ( x )], wherethe function f n [ c ] denotes the dependence of the n -thrate on the external field c ( x ). For simplicity, we willhenceforth denote the rates as functions of x , keeping inmind, however, that the dependence arises via the gra-dient in c ( x ). If all transition rates are independent ofthe external field c ( x ) and equal, i.e., γ n = β with a pos-itive constant β , the spatially homogeneous case stud-ied in [25] is recovered. It is important to notice thatthe model can be easily formulated in two dimensionsas explained in the SM, however, we focus on the one- FIG. 1. Run-time distribution and chemotactic responseof
P. putida . Left: probability distribution function (PDF) φ of run times. Right: cumulative distribution function (CDF)of run-times discriminating between up-gradient (Φ u ) anddown-gradient (Φ d ) runs. The purple-dashed and black-solidcurves in the left panel correspond to fittings with an ex-ponential and a γ -distribution, respectively. A χ -test indi-cates that an exponential distribution is rejected at signif-icance level 0 .
05. The insets show that the same qualita-tive behavior is obtained in simulations with a clock modelas introduced in the main text. Parameters: M = 2, γ = 5 s − , γ ( x ) = [0 . x/L ) + 0 . − s − , L = 750 µ m and v = 30 µ m / s. M γ M ( x ) γ ( x ) γ ( x ) γ ( x ) γ ( x ) γ ( x ) tx ( t ) reversal FIG. 2. Illustration of the internal clock dynamics (left panel)and the resulting particle motion (right panel). The transitionrates γ n ( x ) = f n [ c ( x )] may be concentration dependent. Aparticle reverses its direction of motion, i.e., s ( t ) → − s ( t ),after a full clock cycle has been completed as indicated by ared square on the clock. A segment of a space-time trajectoryis represented by a gray line in the right panel. Additionally,the internal dynamics is overlaid at those instances in time inwhich clock ticks occur. dimensional scenario for simplicity and without loss ofgenerality here. Let P n ( t ) denote the probability to findthe clock in state n at time t . The temporal evolution of P n ( t ) can be expressed via the Master equation [26] dP ( t ) dt = − γ (cid:0) x ( t ) (cid:1) P ( t ) + γ M (cid:0) x ( t ) (cid:1) P M ( t ) (2a) dP n ( t ) dt = − γ n (cid:0) x ( t ) (cid:1) P n ( t ) + γ n − (cid:0) x ( t ) (cid:1) P n − ( t ) (2b)with n ∈ { , , . . . , M } . Notice that this is a closed chainof states. Run-time distribution. – To compute the run-time dis-tribution, we solve for the first passage time of a directedwalk from state 1 to state M within the clock. Let usconsider a particle that moves in direction s at time t and is in state 1. We follow its dynamics in space given by x ( t ) = x + s v ( t − t ) and stop it when it leaves state M .We can obtain the probabilities P n ( t ) by recursion: P ( t ) = e − (cid:82) tt dt (cid:48) γ ( x ( t (cid:48) ) ) , (3a) P n ( t ) = (cid:90) tt dt (cid:48)(cid:48) γ n − (cid:0) x ( t (cid:48)(cid:48) ) (cid:1) P n − ( t (cid:48)(cid:48) ) e − (cid:82) tt (cid:48)(cid:48) dt (cid:48) γ n ( x ( t (cid:48) ) ) . (3b)The run-time distribution is given by φ ( t | x , s ) = γ M (cid:0) x ( t ) (cid:1) P M ( t ); an example of φ for a M = 2 clock isshown in Fig. 1. Spatio-temporal dynamics. – Since we are dealing witha genuine Markov process, the exact, full spatio-temporaldynamics of the problem can be described in terms of aMaster equation for the probability density P ± n ( x, t ) tofind a particle in position x at time t oriented along thedirection s = ± n : (cid:104) ∂ t ± v ∂ x + γ n ( x ) (cid:105) P ± n ( x, t ) = (cid:40) γ M ( x ) P ∓ M ( x, t ) , n = 1 ,γ n − ( x ) P ± n − ( x, t ) , n ≥ . (4)This is a system of 2 M coupled linear partial differentialequations. Notice that the dynamics of P ± ( x, t ) is specialsince it contains the dynamics of reversals. Eq. (4) canbe expressed in a concise matrix form for the vector P ( x, t ) = (cid:0) P +1 ( x, t ) , P − ( x, t ) , . . . , P + M ( x, t ) , P − M ( x, t ) (cid:1) T as follows: ∂ t P ( x, t ) = C · ∂ x P ( x, t ) + L ( x ) · P ( x, t ) , (5)where the matrix C , which encodes the motility ofparticles, is diagonal with the coefficients C ij =( − j v δ ij and L ( x ) contains the transition rates ofthe clock dynamics as well as the reversals. Thetime-dependent solution may be written as P ( x, t ) =exp [( t − t ) M ( x )] P ( x, t ), where we introduced the op-erator M ( x ) = C ∂ x + L ( x ) on the right hand sideand P ( x, t ) abbreviates the initial condition [27]. Clock design & taxis response. – If there is anobservable taxis response, the density ρ ( x, t ) = (cid:80) Mn =1 [ P + n ( x, t ) + P − n ( x, t )] to find a particle at position x should be spatially modulated in the long time limit, i.e., ρ ( x ) = lim t →∞ ρ ( x, t ) should not be constant. We exploitthe fact that the time-independent solution of the Mas-ter equation (4) is unique [28]. In this way, we can easilyverify whether a non-constant ρ ( x ) should be expectedfor a given clock design.We classify clock designs into two categories, homoge-neous and inhomogeneous clocks, addressed separately inthe following. (1) Homogeneous clocks: We call a clock homogeneous if all transition rates γ n ( x ) are identical functions of theconcentration: γ n ( x ) = f [ c ( x )] for n = 1 , , . . . , M withan arbitrary function f [ c ]. We start by considering aspecial case of great relevance: a clock model with onlyone tick ( M = 1). In this particular case, the Masterequation reduces to (cid:16) ∂ t ± v ∂ x (cid:17) P ± ( x, t ) = γ ( x ) (cid:104) P ∓ ( x, t ) − P ± ( x, t ) (cid:105) . (6)Even though the transition rate γ depends on c ( x ), thestationary solution is spatially independent: P ± ( x ) =(2 L ) − . In short, a particle moving at constant speedis unable to detect a chemical gradient if the reversalprocess is controlled by a Poisson process ( M = 1) witha concentration dependent transition rate. Surprisingly,a homogeneous stationary solution, P ± n ( x ) = (2 M L ) − , isalso obtained for all homogeneous clocks with M >
1, de-spite the run-time distributions are γ -shaped. We stressthat homogeneous clocks will, however, lead to a mea-surable run-time bias at the single particle level, i.e. run-times up-gradient differ compared to their down-gradientcounterparts, cf. the inset in Fig. 3a and Eq. (3). Never-theless, there is no accumulation at concentration max-ima nor minima. In short, homogeneous clocks with mul-tiple ticks do not induce a long-time taxis response asillustrated in Fig. 3a,d. We thus conclude that the mea-surement of stationary concentration profiles is an indis-pensable piece of information, whereas the observationof a run-time bias in the individual trajectories provides only insufficient insight into the long-time chemotaxisperformance. (2) Inhomogeneous clocks: If at least two transi-tion rates are different from one another, the station-ary state is not spatially homogeneous and nontrivialstationary profiles of ρ ( x ) can develop (see Fig. 3b,c,e,and f). For arbitrary spatial dependencies of the tran-sition rates γ n ( x ), we cannot find the exact stationarysolution of Eq. (4) analytically, but various approxima-tion methods can be applied.In Ref. [29], a drift-diffusion approximation of theFokker-Planck type [26] ∂ t ρ ( x, t ) (cid:39) − ∂ x (cid:104) f ( x ) ρ ( x, t ) (cid:105) + ∂ x (cid:104) D ( x ) ρ ( x, t ) (cid:105) (7)for the long-time dynamics of the density ρ ( x, t ) was pro-posed in a related context, derived under the assump-tion that the mean distance traversed by a particle inbetween two reversals is shorter than the characteristiclength scales at which the chemical gradient varies. Theapproximation scheme is based on a slow mode reduc-tion of the full Master equation [Eqs. (4)]: as the particledensity ρ ( x, t ) is a conserved quantity, its dynamics isslow, thus enabling the systematic adiabatic eliminationof fast degrees of freedom. In this way, drift f ( x ) anddiffusion D ( x ) can be analytically obtained for any clockmotif [29]. This approximation allows the prediction ofthe density profile ρ ( x ) in the long-time limit based onthe stationary solution of Eq. (7).In this work, we present a perturbative approach thatenables us to predict the stationary solution for the vec-tor P ( x ) including all of its components beyond thescalar, stationary density ρ ( x ). In Figs. 3 and 4, wecompare particle-based simulations with the perturba-tive solution whose derivation is sketched below (see SMfor additional technical details). The central idea of theperturbation theory is that the transition rates are onlyweakly modulated in space allowing us to split the ma-trix L ( x ) into a constant and a spatially dependent part, L ( x )= L + ε L ∆ ( x ), where the second matrix is defined insuch a way that (cid:82) L dx L ∆ ( x ) = 0. Now, we assume that ε is a small parameter allowing to expand the stationarysolution as a power series: P ( x ) = (cid:80) ∞ n =0 P µ ( x ) ε µ . In-serting this ansatz into Eq. (5) and collecting orders in ε yields a systematic way to construct the stationary so-lution. This approach allows us to demonstrate ( i ) thatactive particles controlled by inhomogeneous clocks dis-play a taxis response and ( ii ) how the clock design de-termines the type of response. Figures 3b and 3c showthe emerging density profiles for particles with M = 2-clocks subjected to the same external field, where thetransition rates are γ = α ( x ) and γ = β in Fig. 3b,while γ = β and γ = α ( x ) in Fig. 3c. The importantobservation here is that ρ ( x ) increases monotonically forthe former clock and decreases for the latter one as onemoves in up-gradient direction with respect to c ( x ); no-tably, the average tumbling frequency increases in bothcases in the up-gradient direction. In Fig. 3e and 3f, a -4 -2 FIG. 3. Stationary distributions ρ ( x ) for two homogeneous, ( a ) and ( d ), and four inhomogeneous clocks, ( b ),( c ), ( e ) and ( f ),with M = 2 and M = 5 internal states. Points represent particle-based simulations and lines show the density profiles aspredicted by the perturbation theory. The external field c ( x ), common for all cases shown, is displayed as an inset in panel (d).If all rates are equal (homogeneous clocks), a uniform density profile develops – particles are thus unable to respond to thechemical gradient [panels ( a ) and ( d )]. However, there can be a run-time bias, calculated from Eq. (3) as indicated in the insetof panel ( a ). A nonuniform density profile is observed if the symmetry of the internal dynamics is broken. The probability tofind a particle is proportional to the modulation of the first transition rate α ( x ) [panels ( b ) and ( e )]. By inverting the designof the clock [panels ( c ) and ( f )], the gradient in the resulting density profile switches sign with respect to the former case. Forall cases, the run-time distribution φ ( t ) (shown in each panel as an inset using the same bin width of ∆ t = 0 .
7) is bell-shapedwhen there is more than M = 1 internal state. Note that the shape of the run-time distribution as well as the steepness ofthe density gradient depend on the design of the clock. Parameters in arbitrary units: α ( x ) = c ( x ), c ( x ) = 0 .
75 + 0 . x , β = 1, v = 0 . L = 1. For particle-based simulations, N = 10 particles were simultaneously tracked using a stochastic Eulerscheme [26] with ∆ t = 10 − and density distributions were averaged over time neglecting initial transients. similar scenario is shown for an M = 5-clock for com-parison. While the qualitative trends remain the same,the gradient in the nonuniform density profiles becomessteeper in the M = 5-case.In Fig. 4, a detailed comparison of individual-basedsimulations and the perturbation theory is shown for twoclocks with two ticks ( M = 2), confirming that the pre-sented perturbative approach can indeed predict not onlythe overall density profile ρ ( x ) but also the individualprobability densities P ± i ( x ) to find a particle in state i at position x moving up- or down-gradient, respectively,in the stationary state. Efficiency of taxis response. – To quantify the efficiencyof the taxis response of particles whose reorientation iscontrolled by clocks, we measure the emerging particlecurrent that arises once a chemical gradient is instanta-neously switched on at t = 0, given an initially flat den-sity distribution in a homogeneous environment for t ≤ M for clockswhere γ ( x ) = α ( x ), while all γ n = β with n > λ − spent by the clock between ticks 2 and M bysetting γ n = β = λ ( M − λ is a constant. Inthis way, clocks with larger M correspond to more accu-rate clocks in the sense that the standard deviation of thetime between ticks 2 and M divided by its mean valuescales as 1 / √ M −
1. Applying Eq. (7) to an initially ho-mogeneous density distribution ( ρ ( x, t = 0) = ρ ), wederive the initial current j M ( x, t = 0) = ρ (cid:104) f ( x ) − ∂ x D ( x ) (cid:105) , (8a)= ρ ( M − · v · β + ( M − α ( x ) α ( x ) (cid:2) β + ( M − α (cid:3) · dα ( x ) dx , (8b)where the first line is the general expression and the sec-ond line follows for the particular clock model under con-sideration, derived by application of the drift-diffusionapproximation proposed in [29]. We obtain j = 0 forclocks with only one tick. Fig. 5 shows that the ra-tio ε j = j M /j ∞ increases above zero only for M >
FIG. 4. Comparison of the prediction of the perturbationtheory and individual-based simulations at the level of theprobabilities P i ( x ) = P + i ( x ) + P − i ( x ) to find a particle instate i at position x in the stationary state for the clock mo-tifs shown in Fig. 3b,c. The insets indicate the splitting ofthese probabilities into up- and down-gradient motion, i.e. theprobability densities P ± i . Lines correspond to predictions ofthe perturbation theory, points show individual-based simu-lations. Parameters are identical to Fig. 3. creases with the clock accuracy, saturating for large val-ues of M . Note that the asymmetry of initial densitycurrents is compensated by a nonuniform density distri-bution of particles in the steady state. Concluding remarks. – We studied the ability of activeparticles to perform taxis controlled by internal clocks.Our results reveal that clock designs with homogeneoustransition rates γ n ( x ), i.e., γ n ( x ) = f [ c ( x )] for all n withan arbitrary function f [ c ], cannot generate a taxis re-sponse. It is important to realize that – despite the ab-sence of a taxis response – the frequency of reversals maydepend on the particle position, increasing or decreas-ing up-gradient depending on the clock design. We haveshown that the clock has to fulfill the following require-ments in order to generate a taxis response: ( i ) a numberof ticks M ≥
2, ( ii ) at least one transition rate must de-pend on the external field, i.e., depends on x , and ( iii )inhomogeneous transition rates, i.e., some or all transi-tion rates should be pairwise distinct. Regarding the ef-ficiency of the taxis response, we found that it increaseswith M , i.e. with the clock accuracy, though saturatingfor large values of M .We stress that the obtained taxis responses are not dueto the fact that we focused on velocity-reversing particles. All results hold true if we replace velocity reversals bytumbling events: if – instead of reversing with probabil-ity one – the particle reverses its direction of motion withprobability p and keeps it otherwise, all results reportedhere hold qualitatively true (cf. SM). By putting veloc-ity reversals and tumbling events on an equal footing,it becomes evident that the crucial difference betweenthe chemotactic mechanism of run-and-tumble bacteria,such as E. coli , and the one proposed here is the funda-mentally distinct dynamics of the motility-control mech-anisms, and not the type of turning maneuver (a tumbleor a reversal) which is activated by the control mecha-nism. The two most compelling differences between thesetwo mechanisms can be summarized as follows:(a) In a uniform external field, i.e. in the absence of agradient, run-and-tumble bacteria display an exponentialdistribution of run-times, and thus the temporal sequenceof tumbling events is well-characterized by a single rate.In the here-proposed model, which is phenomenologicallyconsistent with observations of run-and-reverse bacteriasuch as
P. putida , the distribution of run-times is neverexponential for
M >
1, but γ -shaped, cf. Fig. 1. This im-plies that the temporal sequence of reversals (or tumbleevents) cannot be parametrized by just a single rate, anobservation indicating that the dynamics that controlsthe triggering of such events is more complex than a sim-ple Poisson process [18, 20]; we show that it may well berepresented by a clock model.(b) The canonical chemotaxis strategy of run-and-tumble bacteria relies on a memory: when experiencinga temporal increase in chemoattractant concentration,they decrease their tumble frequency. In the proposedclock-controlled taxis mechanism, in contrast to the clas-sical chemotaxis picture, biased up-gradient motion canoccur even when the tumbling frequency increases as the . . . M t a x i s e f fi c i e n c y ε j = j M / j ∞ λ = βM − M α ( x ) λβ FIG. 5. Taxis efficiency, characterized by the relative cur-rent ε j = j M /j ∞ , see text and Eq. (8), as a function of thenumber of ticks M , which measures the accuracy of the clock.The mean frequency λ for the occurrence of the last M − M by choosing β = ( M − λ .The efficiency of sensing an external field increases with theaccuracy of the clock. There is no chemotactic responsefor M = 1. Parameters in arbitrary units (cf. Fig. 3): po-sition of the current measurement ¯ x = 0 . α ( x ) is a linearfunction with α ( x ) = 0 .
75 + 0 . x , λ = 1, v = 0 . L = 1. bacterium moves up-gradient. Moreover, depending onthe clock design, we have observed ( i ) absence of taxis,( ii ) up-gradient motion or ( iii ) down-gradient motion,even though the tumbling frequency was increasing withincreasing concentration in all three cases, as illustratedin Fig. 3. This highlights the relevance of the clock ar-chitecture and indicates that the taxis direction is notdictated by the concentration dependency of the tum-bling frequency, but rather by the clock design. In otherwords, we conclude from our results that it is generallynot sufficient to measure the run-time bias only in orderto infer the chemotaxis strategy of a microogranism butthe measurement of the stationary concentration profilein a chemical gradient provides necessary, complemen-tary information. We recall in this context that homoge-neous clocks, where all transition rates are equal, imply arun-time bias at the single-particle level but will not yielda nonuniform concentration profile, i.e. particles will notaccumulate at concentration maxima nor minima. Notealso that the measurement of these profiles is experimen-tally managable, e.g. via microfluidic maze structures asreported recently [30].Altogether, the proposed clock-controlled taxis mech-anism is a powerful conceptual toolbox that allows us toengineer a large variety of taxis responses, which may alsoplay a key role in the design of simple robots [31] or forcontrolled assembly of active colloids [32, 33]. While theproposed phenomenological model is qualitatively consis-tent with experimental observations of P. putida , it is im- portant to stress that we currently have no further mech-anistic evidence that
P. putida or other bacteria operateby such a clock. We insist on the phenomenological na-ture of the proposed clock model for sensing a signal, itsinternalization and processing, presumably involving cas-cades of biochemical events depicted by stochastic check-points. Identifying the intracellular mechanisms regu-lating the chemotactic responses of
P. putida as well asother bacteria is a long-term experimental challenge be-yond the scope of our current work. We thus hope thatthis study will open the door to a new series of experimen-tal and theoretical works that advance our understandingof the directional navigation of bacteria.
ACKNOWLEDGMENTS
L.G.N., R.G. and F.P. acknowledge financial sup-port from Agence Nationale de la Recherche via GrantNo. ANR-15-CE30-0002-01. L.G.N. was addition-ally supported by CONACYT PhD scholarship 383881and R.G. by the People Programme (Marie Curie Ac-tions) of the European Union’s Seventh Framework Pro-gramme (FP7/2007-2013) under REA grant agreement n.PCOFUND-GA-2013-609102, through the PRESTIGEprogramme coodinated by Campus France. M.H. andC.B. thank the research training group GRK 1558 fundedby Deutsche Forschungsgemeinschaft for financial sup-port. [1] H. C. Berg and D. A. Brown, Nature , 500 (1972).[2] H. C. Berg,
E. coli in Motion (Springer Science & Busi-ness Media, 2008).[3] E. Korobkova, T. Emone, J. M. Vilar, T. S. Shimizu, andP. Cluzel, Nature , 574 (2004).[4] Y. Tu and G. Grinstein, Phys. Rev. Lett. , 208101(2005).[5] M. J. Tindall, E. A. Gaffney, P. K. Maini, and J. P.Armitage, WIREs Syst. Biol. Med. , 247 (2012).[6] A. Celani and M. Vergassola, Proc. Natl. Acad. Sci. USA , 1391 (2010).[7] M. J. Schnitzer, Phys. Rev. E , 2553 (1993).[8] M. Cates, Rep. Prog. Phys. , 042601 (2012).[9] M. Flores, T. S. Shimizu, P. R. ten Wolde, andF. Tostevin, Phys. Rev. Lett. , 148101 (2012).[10] M. Hintsche, V. Waljor, R. Großmann, M. J. K¨uhn,K. M. Thormann, F. Peruani, and C. Beta, Sci. Rep. , 16771 (2017).[11] Z. Alirezaeizanjani, R. Großmann, V. Pfeifer,M. Hintsche, and C. Beta, Sci. Adv. , eaaz6153(2020).[12] M. Theves, J. Taktikos, V. Zaburdaev, H. Stark, andC. Beta, Biophys. J. , 1915 (2013).[13] M. Theves, J. Taktikos, V. Zaburdaev, H. Stark, andC. Beta, Europhys. Lett. , 28007 (2015).[14] Recent experiments with P. putida indicate that the na-ture of the run-time distribution, i.e. an exponential vs.a non-exponential shape, is strongly influenced by the growth conditions of the bacteria: the data presented hereare obtained from bacteria, which are grown on benzoateas a carbon source, whereas in [11], cells cultured in Tryp-tone Broth (TB) and imaged in casamino acid gradientsshowed an exponential run-time distribution.[15] Y. Wu, A. D. Kaiser, Y. Jiang, and M. S. Alber, Proc.Natl. Acad. Sci. USA , 1222 (2009).[16] A. Be’er, S. K. Strain, R. A. Hern´andez, E. Ben-Jacob,and E.-L. Florin, J. Bacteriol. , 2709 (2013).[17] L. Xie, T. Altindal, S. Chattopadhyay, and X.-L. Wu,Proc. Natl. Acad. Sci. USA , 2246 (2011).[18] E. A. Korobkova, T. Emonet, H. Park, and P. Cluzel,Phys. Rev. Lett. , 058105 (2006).[19] E. P. Ipi˜na, S. Otte, R. Pontier-Bres, D. Czerucka, andF. Peruani, Nat. Phys. , 610 (2019).[20] F. Wang, H. Shi, R. He, R. Wang, R. Zhang, and J. Yuan,Nat. Phys. , 710 (2017).[21] J. Sambrook and D. Russell, Molecular Cloning: A Labo-ratory Manual (Cold Spring Harbor: Cold Spring HarborLaboratory Press, 2001).[22] C. S. Harwood, R. E. Parales, and M. Dispensa, Appl.Environ. Microb. , 1501 (1990).[23] C. S. Harwood, M. Rivelli, and L. N. Ornston, J. Bacte-riol. , 622 (1984).[24] O. Pohl, M. Hintsche, Z. Alirezaeizanjani, M. Seyrich,C. Beta, and H. Stark, PLoS Comput. Biol. , e1005329(2017).[25] R. Großmann, F. Peruani, and M. B¨ar, New J. Phys. , Stochastic Methods: A Handbook for theNatural and Social Sciences , Springer Series in Synerget-ics (Springer, 2010).[27] H. Risken,
The Fokker-Planck Equation: Methods of So-lution and Applications (Springer Verlag, 1996).[28] N. Van Kampen,
Stochastic Processes in Physics andChemistry , North-Holland Personal Library (Elsevier Sci-ence, 2011).[29] L. G. Nava, R. Großmann, and F. Peruani, Phys. Rev. E , 042604 (2018).[30] M. M. Salek, F. Carrara, V. Fernandez, J. S. Guasto, andR. Stocker, Nat. Commun. , 1877 (2019).[31] M. Mijalkov, A. McDaniel, J. Wehr, and G. Volpe, Phys.Rev. X , 011008 (2016).[32] T. B¨auerle, A. Fischer, T. Speck, and C. Bechinger, Nat.Commun. , 3232 (2018).[33] H. Karani, G. E. Pradillo, and P. M. Vlahovska, Phys.Rev. Lett. , 208002 (2019). UPPLEMENTAL MATERIALA novel approach to chemotaxis: active particles guided by internal clocks
Luis G´omez Nava, Robert Großmann, Marius Hintsche, Carsten Beta, and Fernando Peruani
Appendix A: Experiments with P. putida1. Cell culture
P. putida
KT2440 was streaked on 1.5% agar (AppliChem, Germany) containing LB medium (AppliChem, Ger-many) and grown at 30 ◦ C. A single-colony isolate was used to inoculate 10 ml of M9 medium [1] supplemented with5 m M sodium benzoate as a carbon source [2, 3]. The culture was grown for 36 h on a rotary shaker (300 min − ,30 ◦ C) to stationary phase. It then was diluted 1:100 into 25 ml of fresh M9 medium and grown 12 h to an OD of approximately 0.3 in exponential phase. Bacteria were washed three times by centrifugation at 1000g for 10 minand carefully resuspended in 10 ml motility buffer (1 × − M potassium phosphate, 6 . × − M NaCl, 1 × − M EDTA and 0.5%w/v glucose; pH 7.0). Cells were diluted further to an OD of 0.05 before filling them into thechemotaxis device.
2. Microfluidics and imaging
We used a µ -Slide Chemotaxis 3D (ibid., Martinsried, Germany) to generate stable linear gradients of benzoate.Filling was performed as detailed in [4] with concentrations in the source and sink reservoir of 0 . M and 5 m M sodiumbenzoate, respectively. Imaging was done using an IX71 inverted microscope with a 20 × UPLFLN-PH objective (bothOlympus, Germany) in phase contrast mode with an attached Orca Flash 4.0 CMOS camera (Hamamatsu Photonics,Japan). Video data was acquired for two minutes at 20 Hz, approximately one hour after filling the channel. Thegradient linearity and stationarity was checked previously [4]. The viewport and focal plane were set in the center ofthe channels observation region, about 35 µ m from top and bottom. Image processing and tracking was performed aspreviously described in [4].
3. Filtering and smoothing
To filter out damaged, non-motile cells and those swimming on very wobbly trajectories, we discarded tracks basedon a number of criteria. Tracks below a mean speed of 5 µ m s − , longer than 20 s, above the 80th percentile of mediancurvature and with a total displacement below 5 µ m are discarded. To further smooth these tracks, we applied a 3point, second order Savitzky-Golay filter. Because tracks often begin or end with a reversal we cut of the first andlast 0 .
4. Run-and-tumble recognition
Characterization of reversal events was performed with the algorithm of Theves et al. [5]. This algorithm de-termines reversals by evaluating local extrema in the speed and turn rate time series and counting a reversal forsufficiently deep peaks and troughs.
Appendix B: Active particles with clock-controlled tumbling in one and two dimensions
As indicated in the main text, the appearance of a stationary profile ρ ( x ) does not depend on the responsetriggered by the clock, i.e., whether this is a velocity reversal or a tumbling event. Here, we show first – usingindividual based simulations – that the same qualitative results, described in the main text, are recovered if reversalsare replaced by tumbling events. We also explain how the model can be generalized to two dimensions and provide nu-merical evidence that indicates that the clock mechanism is also effective to produce a taxis response in two dimensions. FIG. B.1. Numerical simulations of particles in a one-dimensional box with reflecting boundary conditions using a lineargradient as done in Fig. 3 in the main text, but with a tumbling mechanism instead of reversals: the probability to reverseis p = 1 /
2, and q = 1 / M to tick 1. We use α ( x ) = c ( x ), c ( x ) = x/L , c l = 0, c u = 1, β = 1, v = 0 .
01 and L = 1.
1. Clock-controlled tumbling in one dimension
We implement tumbling in one dimension as follows: When a transition between tick M to 1 occurs, we flip thedirection of motion of the particle with a probability of p = 1 / q = 1 − p .Fig. B.1 displays the stationary profiles obtained in simulations with this one-dimensional tumbling mechanism. Thesimilarity with the reversal dynamics is evident (cf. Fig. 3 of the main text). The results show that nontrivial stationaryprofiles appear independently of the mechanism how the direction of motion changes (reversal vs. tumbling), therebyproviding a proof of concept that clock-controlled tumbling allows for guided transport and chemotaxis.
2. Clock-controlled tumbling in two dimensions
The appearance of the profile ρ ( x ) as a result of the symmetry breaking of the internal clock of particles isnot particular for the one-dimensional model studied in the main text. To illustrate this fact, we implemented atwo-dimensional version of our model using individual-based simulations.A nontrivial profile also emerges in two dimensions whenever an inhomogeneous clock is used. In two-dimensions,the motility is modeled by the differential equation ˙ x ( t ) = v b e ( t ) , (B1)where x ( t ) is the position of the particle, b e ( t ) = (cos ϕ ( t ) , sin ϕ ( t )) is a unit vector pointing in the direction of motionand v is a constant. The direction of motion ϕ ( t ) is changed each time the transition from clock tick M to 1 occurs.A new direction of motion ϕ is picked up at random from the uniform probability distribution p ( ϕ ) = 1 / (2 π ) on theinterval [0 , π ].The obtained results are shown in Fig. B.2 for the transition rate α ( x, y ) = c ( x, y ) = x/L , and β = 1. Asobserved in the one-dimensional version of the model, the emergence of a taxis response also requires the use of aninhomogeneous clock in two dimensions. FIG. B.2. Stationary probability densities ρ ( x, y ) obtained from individual-based simulations using N = 10 particles thatmove at constant speed v = 0 .
01 inside a box of dimensions L × L (with L = 1) with reflecting boundary conditions. Theleft panel displays the stationary profile ρ ( x, y ) that emerges for a M = 2 clock with γ ( x, y ) = γ ( x, y ) = α ( x, y ). The rightpanel corresponds to an inhomogeneous M = 2 clock with γ ( x, y ) = α ( x, y ) and γ ( x, y ) = β , with β = 1. In both panels, α ( x, y ) = c ( x, y ) = x/L . For each panel, the run-time distribution is plotted as an inset. Appendix C: Perturbation theory1. Weakly modulated transition rates
In this section, we explain the mathematical details of the perturbation theory which can be used to approximatethe stationary solution of the Master equation corresponding to the model under consideration. Analogous to thenotation in the main text, we abbreviate the probabilistic dynamics as follows: ∂ t P ( x, t ) = C · ∂ x P ( x, t ) + L ( x ) · P ( x, t ) . (C1)The components of the vector P ( x, t ) = (cid:0) P +1 ( x, t ) , P − ( x, t ) , P +2 ( x, t ) , P − ( x, t ) , . . . , P + M ( x, t ) , P − M ( x, t ) (cid:1) T (C2)are the probability density functions to find a particle in position x at time t with an internal state n , where n ∈{ , , . . . , M } , moving from left to right (+) or, conversely, from right to left ( − ). The matrix C encodes the motilityof particles and L ( x ) denotes space-dependent transition rates constituting the clock dynamics. The former matrix C is diagonal and its components read C ij = ( − j v δ ij . The concrete form of L ( x ) depends on the specific clock modelunder consideration. As an example, we specifically state these matrices for the clock model shown in Fig. 3b in themain text: C = − v v − v
00 0 0 v , L ( x ) = − α ( x ) 0 0 β − α ( x ) β α ( x ) 0 − β α ( x ) 0 − β . (C3)The dynamics has to be supplemented by appropriate boundary conditions as well as the corresponding initial distri-bution of particles. Since we will derive the stationary solution only, the initial condition is of minor importance inthis context. In this section, we will study periodic boundary conditions only, i.e., we consider all space-dependentfunctions to be periodic, where the linear system size L determines the period. We clarify in Appendix C 2 that thisapproach can be directly applied to reflecting boundary conditions as well.The central assumption of the perturbation theory is that the position dependence of the transition rates is weak.More precisely, we split the matrix L ( x ) into a constant and a spatially dependent part according to L ( x ) = L + ε ∆( x ) L S , (C4)where the scalar function ∆( x ) encodes the spatial dependence of the rates, the constant matrix L S states which ofthe components are spatially modulated and ε is defined by the condition:0 = Z + L/ − L/ dx ∆( x ) . (C5)For the specific clock model mentioned above, the components of the matrices L and L S read L = − α β − α β α − β α − β , L S = − − . (C6)Henceforth, we consider ε to be small such that the solution of the Master equation can be expanded in a Taylorseries in this expansion parameter. In short, we propose the ansatz P ( x, t ) = ∞ X µ =0 P µ ( x, t ) ε µ . (C7)Inserting into Eq. (C1) yields the following conditional equations for P µ ( x, t ): ∂ t P µ ( x, t ) = C · ∂ x P µ ( x, t ) + L · P µ ( x, t ) + ∆( x ) · ( , µ = 0 , L S · P µ − ( x, t ) , µ ≥ . (C8)In the following, we sketch the derivation of the solution of these equations to the lowest, nontrivial order.In the lowest order perturbation theory, all transition rates are constant and, hence, the stationary state P isspatially homogeneous. It is thus easily found from the condition L · P = 0 which allows us to determine P up toa normalization constant that is given by the particle number conservation – the sum over all components of P mustequal the uniform probability distribution function ρ = 1 /L .The stationary solution of the first, nontrivial order is determined by the following differential equation: (cid:18) − C ddx − L (cid:19) · P ( x ) = ∆( x ) L S · P . (C9)There is a mathematical subtlety since this equation does not possess a unique solution as the operator on the lefthand side is not invertible. To see this, let P ( x ) = ¯ P ( x ) be particular solution. Consequently, the function P ( x ) =¯ P ( x ) + c P is a solution as well for all real coefficients c ∈ < . This ambiguity is resolved by reconsidering thenormalization condition as detailed below. First, we integrate Eq. (C9) over the entire system yielding1 L Z + L/ − L/ dx (cid:18) − C ddx − L (cid:19) · P ( x ) = −L · p = 0 , (C10)where p µ = L − R + L/ − L/ dx P µ ( x ) was introduced. Therefore, we obtain the intermediate result p ∝ p , since L · P = 0and, thus, L · p = 0. However, we have already normalized the lowest order solution p correctly – the sum overall components of p equals 1 /L – and, consequently, we conclude that p must equal zero. In short, the first ordercorrection P ( x ) has to be determined such that the integral over this function vanishes,0 = Z + L/ − L/ dx P ( x ) , (C11)thereby ensuring the solvability of Eq. (C9) unambiguously. The explicit solution P ( x ) = 1 L ∞ X n = −∞ ,n =0 ∆ n e − iq n x G ( q n ) ·L S · P (C12a)= 2 L ( ∞ X n =1 < h ab G ( q n )∆ n i cos( q n x ) + = h ab G ( q n )∆ n i sin( q n x ) ) ·L S · P (C12b)is found by Fourier transformation of Eq. (C9) with respect to the spatial variable x . It is easily verified by insertinginto this equation. In the expression above, we abbreviated the wavenumbers q n = 2 πn/L with n ∈ Z \{ } , theFourier coefficients ∆ n of the spatial modulation ∆( x ) defined by∆ n = Z + L/ − L/ dx e iq n x ∆( x ) (C13) xα r ( x ) L ˜ L/ α l α u α xα p ( x ) - ˜ L Lα l α u α FIG. C.1. Illustration of the analogy of reflecting boundary conditions (left) and a corresponding system with periodic boundaryconditions which has the double length with respect to the former. In the stationary state, the probability distribution P ( x )on the interval x ∈ (0 , ˜ L ) will be similar if α p ( x ) is constructed from α r ( x ) by mirror reflection, cf. Eq. (C15). as well as the propagator G ( q ) = ( iq C − L ) − , defined for q = 0.In short, the condition L · P = 0 in combination with the expression for the first order, nontrivial correc-tion P ( x, t ), cf. Eq. (C12), determine the stationary solution to first order P ( x ) = P + ε P ( x ) + O (cid:0) ε (cid:1) , (C14)which can be rewritten for a specific clock model by specification of ∆( x ) as well as the matrices L as well as L S asexplained in Appendix C 3.
2. Reflecting vs. periodic boundary conditions
So far, we considered a perturbation theory for stationary density profiles which is applicable to systems withperiodic boundary conditions. However, one can utilize the expressions derived above also for systems with reflectingboundary conditions as discussed in the following, cf. Fig. C.1. For the sake of clarity, we consider a linear system [sys-tem size ˜ L and x ∈ (0 , ˜ L )] where one of the transition rates, denoted by α r ( x ), depends explicitly on space. Nowthe question is: what happens at the boundaries of the system? If a particle, moving from right to left, touches theboundary in x = 0, it will be reflected. In other words, a particle, which was moving into the direction of decreasingtransition rates α r ( x ), will move towards increasing transition rates after the reflection. Now, this situation mayanalogously be observed in a periodic system which is as twice as long [ x ∈ ( − ˜ L, ˜ L )], where the space dependenttransition rate is replaced by α p ( x ) = ( α r (+ x ) , x ≥ ,α r ( − x ) , x < . (C15)In the stationary state, the number of particles passing through the point x = 0 from left to right is equal to the numberof particles passing this point moving in the opposite direction due to the symmetry of the problem. Therefore, theparticle flux in x = 0 vanishes. This line of arguments holds similarly for the other boundary at x = ˜ L . Accordingly,the dynamics of the two systems, one with reflecting and the extended system with periodic boundary conditions,is completely equivalent on the interval x ∈ (0 , ˜ L ). In short, we can map every system with reflecting boundaryconditions onto an extended system with periodic boundary conditions and a space dependent transition rate α p ( x )that is constructed by mirroring α r ( x ) as expressed by Eq. (C15).
3. An example
To fix ideas, we study the particular solution for the clock model shown in Fig. 3b in the main text for a lineargradient as it is represented in Fig. C.1. The matrices L and L S were already given explicitly in Appendix C 1. Thelowest order solution, to be determined from L · P = 0, reads P = 12 ˜ L ( α + β ) ββα α . (C16) . . . .
270 0 . . . . P + ( x ) , P − ( x ) , P + ( x ) , P − ( x ) position x P ± , ( x ) α ( x ) α ( x )0 . .
030 0 . x ρ ( x ) c l c u . x c ( x ) ( a ) . . . . position x P +1 ( x ) P − ( x ) P +2 ( x ) P − ( x ) α ( x ) β . .
030 0 . x ρ ( x ) ( b ) . . . . position x P +1 ( x ) P − ( x ) P +2 ( x ) P − ( x ) βα ( x ) 0 . .
030 0 . x ρ ( x ) ( c ) FIG. C.2. Stationary distributions for the same motifs presented in Fig. 3 of the main text with M = 2. Points representparticle-based simulations and lines show the density profiles as predicted by the perturbation theory. In the large panels, thestationary probability density P ± n ( x ) to find a particle at position x with orientation s = ± n is plotted.Insets represent the stationary probability density ρ ( x ) obtained from the sum over all components P ± n ( x ). The externalfield c ( x ), common for (a), (b) and (c), is shown as an inset in panel ( a ). Parameters in arbitrary units: α ( x ) = α b (1 + c ( x )), α b = 9 . c ( x ) = 0 . x/L , c l = 0, c u = 0 . β = 10, v = 1, M = 2, L = 1. For particle-based simulations, N = 10 particleswere simultaneously tracked using a stochastic Euler scheme with ∆ t = 10 − and density distributions were averaged over timeneglecting initial transients. The parameter α denotes the mean transition rate α = ( α l + α r ) /
2, i.e., α r ( x ) = α + ε ∆( x ). Further, we parametrizethe spatial modulation ∆( x ) according to ∆( x ) = − x ˜ L . (C17)Now, we can make immediate use of the equations derived in Appendix C 1, however, care must be taken that thesystem size was reparametrized: ˜ L = L/
2. Accordingly, the wavenumbers q n are determined by q n = 2 πn/L = πn/ ˜ L .To utilize Eq. (C12), we need to calculate the Fourier coefficients ∆ n corresponding to the spatial modulation ∆( x ),which read ∆ n = Z + L/ − L/ dx e iq n x ∆( x ) = − Ln π h ab − ( − n i (C18)for n = 0 and ∆ = 0, consistently with the condition (C5). Due to the symmetry ∆( − x ) = ∆( x ) of the spatialdependence of the transition rate, these Fourier coefficients are real numbers. What is left to be done is the calculationof the propagator G ( q ) = ( iq C − L ) − – obtained from the inversion of a (4 × P ( x ) = 1 L ( α + β ) · ββα α + ε · L ∞ X n =1 β ∆ n α + β +( q n v ) cos( q n x ) β − α β − α β + α β + α + q n v sin( q n x ) − − + O (cid:0) ε (cid:1) . (C19)Alternatively, one may rewrite the stationary probability density to first order in perturbation theory by insertion ofthe Fourier modes ∆ n as well as the original system size ˜ L : P ( x ) ’
12 ˜ L ( α + β ) · ββα α − ε · π ∞ X n =0 β/ (2 n + 1) α + β +( q n +1 v ) cos( q n +1 x ) β − α β − α β + α β + α + q n +1 v sin( q n +1 x ) − − . (C20)The stationary probability distribution for the position of a particle, irrespective of its internal state, is eventuallygiven by the sum over all components of the vector P ( x ). The perturbation theory results are compared to individualbased simulations in the main text. [1] J. Sambrook and D. Russell, Molecular Cloning: A Laboratory Manual (Cold Spring Harbor: Cold Spring Harbor LaboratoryPress, 2001).[2] C. S. Harwood, M. Rivelli, and L. N. Ornston, J. Bacteriol. , 622 (1984).[3] C. S. Harwood, R. E. Parales, and M. Dispensa, Appl. Environ. Microb. , 1501 (1990).[4] O. Pohl, M. Hintsche, Z. Alirezaeizanjani, M. Seyrich, C. Beta, and H. Stark, PLoS Comput. Biol. , e1005329 (2017).[5] M. Theves, J. Taktikos, V. Zaburdaev, H. Stark, and C. Beta, Biophys. J.105