A Mathematical Description of Bacterial Chemotaxis in Response to Two Stimuli
AA Mathematical Description ofBacterial Chemotaxis in Response to Two Stimuli
Jeungeun Park ∗ Zahra Aminzare † Abstract
Bacteria are often exposed to multiple stimuli in complex environments, and their efficientchemotactic decisions are critical to survive and grow in their native environments. Bacterialresponses to the environmental stimuli depend on the ratio of their corresponding chemorecep-tors. By incorporating the signaling machinery of individual cells, we analyze the collectivemotion of a population of
Escherichia coli bacteria in response to two stimuli, mainly serineand methyl-asparate (MeAsp), in a one-dimensional and a two-dimensional environment. Undersuitable conditions, we show that if the ratio of the main chemoreceptors of individual cells,namely Tar/Tsr is less than a specific threshold, the bacteria move to the gradient of serine,and if the ratio is greater than the threshold, the group of bacteria move toward the gradientof MeAsp. Finally, we examine our theory with Monte-Carlo agent-based simulations.
Key words.
Chemotaxis, Multi-scale dynamics, Population dynamics, Intracellular decision mak-ing, Fokker-Planck equations, Advection-diffusion equations, Monte-Calro simulations.
Mathematics Subject Classification (2020).
The preferred movement of a bacterium along the gradient of chemical substances, the so-calledchemotaxis, includes a directed movement (run) and a relatively short random turning (tumble).See e.g., [1] and [2] for
Escherichia coli ( E. coli ) and
Salmonella typhimurium chemotaxis. Eachbacterium carries an internal state which may be modeled by a system of ordinary differentialequations. In the presence of a stimulus in the environment, each cell changes its direction atrandom, with a tumbling rate which depends on the internal state, biasing moves toward morefavorable environments or away from noxious substances.In their natural environment, bacteria are often exposed to multiple chemical stimuli. To navigatetoward a favorable environment, they choose their direction of movement based on environmentalperception, individual preferences, and interaction with others. Also, each individual’s decisioncharacterizes the behavior of a group of bacteria. Thus, understanding how bacterium choosesbetween multiple stimuli is essential to study bacterial chemotaxis at a population level. ∗ Department of Mathematical Sciences, University of Cincinnati, OH, USA. [email protected] † Department of Mathematics, University of Iowa, IA, USA. [email protected] a r X i v : . [ q - b i o . Q M ] J un hemical signals are often detected via five main chemoreceptors of E. coli , namely Tar, Tsr, Tap,Trg, and Aer [3]. In [4], where responses of
E. coli to two chemoattractant signals are demonstrated,it is shown that the expression levels of the most abundant receptors, Tar and Tsr, are determinedby the bacterial density in a batch-mode culture within the growth phase; in turn, the ratio of thesereceptors differentiates their chemical preferences. Inspired by the experimental results of [4], ourgoal of this work is to incorporate the bacterial decision-making process into a mathematical modelinvestigate the corresponding collective behavior observed in [4].Although modeling the movement of bacteria in response to one stimulus has been studied verywell (see [5] for a review on multi-scaling model approaches for chemotaxis), a few studies haveaddressed the behavior of bacteria in the presence of multiple chemical stimuli.In this work, we consider a population of bacteria in a one-dimensional and a two-dimensionalspatial domain occupied by two stimuli that their temporal rates are assumed to be zero. Weemploy a Fokker-Planck type master equation (also known as balance equation [6]) to describe thebacterial chemotaxis. This (microscopic) model enables us to incorporate the internal dynamics of
E. coli representing the chemotaxis signaling pathway [7, 8]. However, since this equation is notmathematically tractable, using techniques from [9], we derive a (macroscopic) advection-diffusionequation, which is analogous to the Keller-Segel model [10], in which the temporal rate of thechemical concentrations is zero, as is assumed here. Note that unlike the microscopic equation, themacroscopic equation does not depend on the internal dynamics of bacteria explicitly, however, itscoefficients are expressed in terms of the internal state variables.Several works have been studied the behavior of bacteria in response to external signals using themicroscopic and macroscopic models. For instance [9, 11] studied
E. coli chemotaxis in a one-dimensional and n − dimensional space in response to single stimulus, respectively. These studieswere generalized in [12] to multiple time-dependent signals by applying a general type of receptorbased-response laws [13, 14]. However, these works considered a toy model for the internal dy-namics. In [15, 16], the authors studied the microscopic and macroscopic dynamics of bacterialchemotaxis in the presence of a single signal in a one-dimensional space and for any arbitrary in-ternal dynamics. In [17], the authors incorporated particularly E. coli signalling pathway [7] into aone-dimensional macroscopic equation to explain observations from [18, 19, 20] in which the ratioof Tar and Tsr affects bacterial thermotaxis and pH taxis.Our contributions towards understanding the dynamics of a population of bacteria in response totwo stimuli are as follows. First, we incorporate a relatively general class of one-dimensional inter-nal dynamics into a one- and a two-dimensional microscopic equation. Under suitable conditions,we derive a macroscopic equation from the corresponding microscopic equation. Second, we useour macroscopic model for a population of
E. coli with a mechanistically realistic, while a mathe-matically tractable, model of internal dynamics and analyze the response of
E. coli to two stimuliin a one- and a two-dimensional environment. Numerical solutions of our model agree well withMonte-Carlo agent-based simulations. Besides, our mathematical model qualitatively captures theexperimental results of [4].The remainder of the paper is organized as follows. In Section 2, we first review the internal dy-namics of
E. coli which describe how the cells can produce runs and tumbles. Then, given a generalinternal dynamics of bacteria, we introduce a (forward) Fokker-Planck equation which describesthe dynamics of a probability distribution of a population of bacteria. In Section 3 (respectively,Section 5), we first derive a one-dimensional (respectively, two-dimensional) advection-diffusion2quation which approximates the Fokker-Planck equation with a general internal dynamics. Then,we focus on a population of
E. coli with a specific internal dynamics. In Section 4, (respectively,Section 6), for different combinations of stimuli, we solve the one-dimensional (respectively, two-dimensional) advection-diffusion equation numerically and compare the solutions with Monte-Carloagent-based simulations. In Section 7, we conclude with a brief summary and discussion of futuredirections. In Appendix A, we summarize the models with parameter values that we use for theinternal dynamics in Section 2 and for the derivation of the macroscopic equation in Sections 3 and5. The appendix also provides a brief description of the Monte-Carlo agent-based simulation andan overview of our numerical simulations with input data.
E. coli bacteria
We briefly review the internal dynamics of
E. coli which transfer a signal of the environment intoa run (see [7, 21, 22] for more details). Then, we derive a probabilistic equation which describes microscopic dynamics of a population of bacteria with a given internal dynamics, [9, 15, 23]. Later,in the following section, we use the microscopic equation to derive a macroscopic equation whichapproximates the dynamics of a population of bacteria by integrating the internal dynamics of allthe bacteria.
E. coli : A brief review
E. coli bacteria use four to six helical flagella that are connected to rotary motors in their cell wallto swim. Their swimming patterns are characterized as a random walk, consisting of long runs( ∼ ∼ . S , and the tumbling rate, denoted by λ ,represent the input and the output, respectively. As explained above, binding the ligand to thereceptor inhibits the activity of CheA, denoted by a . On the other hand, the methyl group (denotedby m ) in the receptors enhances the activity of a . Therefore, a = G ( S, m ) can be described as anincreasing function of m and a decreasing function of S . CheA CheA-P
CheR + CH - CH CheB
CheB-P
CheY-P CheZ CheY receptor
CheW a m S output (tumbling rate) input (ligand) ligand Figure 1:
Left:
E. coli signaling path-way. Binding ligands to receptors, thesignal is transduced to the flagellar mo-tor via six cytoplasmic chemotaxis pro-teins. Right: An input-output rep-resentation of
E. coli signaling path-way. The internal signaling pathwayshown in left is reduced to the inter-action between the methylation level m and the kinase activity a . This inter-action, which depends on ligand con-centration S (input), controls the motorrotation by changing the tumbling rate(output). See Section 2.1 for detaileddescription. As described earlier, the kinase activity of CheA enhances the CheB-P level, and CheB-P reducesthe methylation level of the receptors. Consequently, the kinase activity a reduces the methylationlevel m , indirectly. So, the dynamics of m can be described by dm/dt = F ( a ), where F is adecreasing function of a .There are several models for F and G as in [7, 8, 35, 36]. Since our results can be applied to almostany scalar functions F and G , for ease of calculations, we choose a simple model, as discussed in(1) and (2) below.Note that the tumbling rate λ is controlled by the level of CheY-P, which is affected by the kinaseactivity. Therefore, λ can be modeled by an increasing function of a , as described in (5) below.4ollowing the experimental set up in [4], we consider two stimuli: S and S , which, respectively,stand for methyl-asparate (MeAsp) and serine, and can be sensed by chemoreceptors Tar andTsr. Furthermore, since the experiments in [4] are designed to keep the external signals S and S constant in time, we assume that S and S only depend on the spatial variable x and areindependent of time t : S = S ( x ) and S = S ( x ) . Following [17, 37, 38], we let a heterogeneous Monod-Wyman-Chageux (MWC) model [39] describethe kinase activity a : G ( S , S , m ) = 11 + η ( m ) η ( S ) η ( S ) , (1)where η ( m ) η ( S ) η ( S ) is derived from the total free energy difference between the active andinactive states. According to [7, 22, 40, 41, 42], the methylation-dependent free energy gives η ( m ) = exp( N α ( m − m )) , where N is the number of the responding receptor dimers in the cluster, and α and m denotethe free-energy per added methylation group and a reference methylation level, respectively. Theligand-depdent free-energy obtains η i ( S i ) = (cid:18) S i /K iI S i /K iA (cid:19) Nr i , where K iI and K iA are the dissociation constants of the corresponding ligand ( i = 1 for MeAsp, i = 2 for serine) to the inactive and the active receptor ( i = 1 for Tar, i = 2 for Tsr). The constantparameters r and r are the fraction of receptors Tar and Tsr in the receptor cluster, respectively.We assume that r + r = 1 and r N and r N are the number of the receptors binding to thecorresponding ligand.The average methylation level, m , of receptors evolves slowly and can be described by the followingequation [7, 22]: dmdt = F ( a ) = a − aτ a , (2)where τ a (cid:29) a is a constant which represents the adaptation level of a , i.e.,when a > a , dm/dt < m and consequently a decrease. When a < a , dm/dt > m and consequently a increase.It is more convenient to use a as a state variable instead of the methylation level m . Taking timederivative of a gives: dadt = ∂a∂m dmdt + ∂a∂S ∇ x S · d x dt + ∂a∂S ∇ x S · d x dt . (3)Using (1) for a = G ( S, m ), we obtain ∂a∂m = αN a (1 − a ) , ∂a∂S i = N a ( a − r i /K iI − /K iA (1 + S i /K iI )(1 + S i /K iA ) . For i = 1 ,
2, we assume that for any x , K iI (cid:28) S i ( x ) (cid:28) K iA ,
5s in [7, 22]. This assumption guarantees scale-invariant behavior of
E. coli in response to externalsignals, which was mathematically predicted in [43] and experimentally verified in [44]. Scale-invariance property of a system means that the system does not distinguish between an input(here, S or S ) and its scaled version (e.g., p S or p S ). For more details see [8] and [45]. Usingthis assumption, we make the following approximation1 /K iI − /K iA (1 + S i /K iI )(1 + S i /K iA ) ≈ S i . Therefore, dadt = αN a (1 − a ) a − aτ a + N a ( a − (cid:16) r ∇ x S · d x /dtS + r ∇ x S · d x /dtS (cid:17) . (4)From now on, we use (4) as a one-dimensional internal dynamics of E. coli bacteria.
Remark 1.
E. coli bacteria can also sense pH changes, and their internal dynamics during pH taxisis analogous to that during chemotaxis. For example, according to [17, 19, 20], Tar receptors areattracted to a decrease of pH, but Tsr receptors show the opposite response. Taking into accounttwo chemical stimuli with different pH levels, we can apply the heterogeneous MWC model anduse the following assumptions to derive the internal dynamics for pH: K I (cid:28) S ( x ) (cid:28) K A , and K A (cid:28) S ( x ) (cid:28) K I , which yield 1 /K I − /K A (1 + S /K I )(1 + S /K A ) ≈ S , and 1 /K I − /K A (1 + S /K I )(1 + S /K A ) ≈ − S . Remark 2.
In this work, we are interested in the total receptor kinase activity of the entirereceptor cluster. Thus, we do not consider two different methylation dynamics for two differenttype of receptors as in [17].As a result of the slow adaptation process (2), bacteria use their methylation state as a short-termmemory store to compare changes of stimuli temporarily during the run. This process helps thebacteria to run or tumble effectively toward their preferred location. Based on the experimentalresults, on the transition to the tumbling state and rotational diffusion, the tumbling rate functioncan be described as λ ( a ) = λ + 1 τ (cid:16) aa (cid:17) H , (5)where λ , H, and τ denote the rotational diffusion, the Hill coefficient of flagellar motor’s responsecurve, and the average run time, respectively, and a is as given in (2). Note that since a dependson S , we may write λ = λ ( a, S ) (see Section 2.2 below). More details about the physical meaningof these parameters can be found in [17, 22, 40]. The parameter values are shown in Table 1.6 .2 Deriving a Fokker-Planck equation describing a population of bacteria In what follows, we describe the motion of a population of bacteria by incorporating their internaldynamics.Let p ( x , a , ν , t ) be a probability density function describing a population of bacteria, modeled ina 2 N + M + 1 dimensional phase space, where time t ∈ R , x = ( x , . . . , x N ) ∈ R N (we willspecialize to N = 1 ,
2) denotes the position of a cell centroid, a = ( a , . . . , a M ) ∈ A ⊂ R M (we willspecialize to M = 1) denotes the internal dynamics of the cell, and ν = ( ν , . . . , ν N ) ∈ V ⊂ R N denotes its velocity, d x /dt = ν . The vector S ( x , t ) = ( S ( x , t ) , . . . , S K ( x , t )) ∈ R K represents theconcentration of extracellular signals in the environment (we will assume that S only depends on x as in Section 2.1 and [4]).Let the following system of ordinary differential equations describe the evolution of the intracellularstate, in the presence of the extracellular signal S : d a dt = f ( a , S ) , (6)where f : R M × R K → R M is a continuously differentiable function with respect to each component,i.e., f ∈ C ( R M × R K ).Assuming constant velocity, dν i /dt = 0, the evolution of p = p ( x , a , ν , t ) with turning rate λ = λ ( a , S ) is governed by the following forward Fokker-Planck equation describing a velocity-jumpprocess [6, 23]: ∂p∂t + ∇ x · ν p + ∇ a · f p = − λ ( a , S ) p + (cid:90) V λ ( a , S ) T ( a , ν , ν (cid:48) ) p ( x , a , ν (cid:48) , t ) d ν (cid:48) , (7)where the non-negative kernel T ( a , ν , ν (cid:48) ) is the probability that the bacteria changes the velocityfrom ν (cid:48) to ν , and (cid:90) V T ( a , ν , ν (cid:48) ) d ν (cid:48) = 1 . Equation (7) is not tractable mathematically and is hard to be validated by typical experimentaltechniques. The goal is to use the microscopic model (7), and derive a macroscopic model forchemotaxis in a one-dimensional space (in Section 3) and a two-dimensional space (in Section 5),i.e., an equation for the marginal density n ( x , t ) = (cid:90) V (cid:90) A p ( x , a, ν , t ) da d ν , with N = 1 or 2, M = 1, and K = 2; n ( x , t ) is the number of individuals which at time t arelocated at position x , whatever their internal dynamics and velocity are.Note that our theory works for any arbitrary K . However, we are interested in two extracellularsignals, so we only consider K = 2. In this section, we assume that the bacteria move in a one-dimensional space, i.e., a finite interval[0 , L ] where we assume L is sufficiently large. We let p ± ( x, a, t ) = p ( x, a, ± ν, t ) denote the density of7he bacteria, located at x ∈ [0 , L ], moving to the right and left, respectively; and let f ± = f ± νf describe their corresponding internal state. Here, ν > ν is constant. Then the Fokker-Planck equation (7) becomes ∂p + ∂t + ν ∂p + ∂x + ∂∂a (cid:2) f + ( a, S ) p + (cid:3) = 12 λ ( a, S )( p − − p + ) , (8) ∂p − ∂t − ν ∂p − ∂x + ∂∂a (cid:2) f − ( a, S ) p − (cid:3) = 12 λ ( a, S )( p + − p − ) . (9)Following [15], under a decay condition for p ± , some conditions on the internal dynamics (forexample, shallow conditions for the stimuli– see Proposition 1 below), moment closure techniques,and parabolic scaling, a general advection-diffusion equation for the marginal density n ( x, t ) = (cid:90) A ( p + ( x, a, t ) + p − ( x, a, t )) da can be derived from Equations (8)-(9) as follows ∂n∂t = ∂∂x (cid:18) ν α ∂n∂x − α B ν α ( A − α ) n (cid:19) . (10)Here, α i , A i , and B i are the Taylor constants of λ , f , and f , respectively: λ = α + α a + · · · ,f = A + A a + · · · ,f = B + B a + · · · . All the Taylor constants depend on S = S ( x ) and we assume that A = 0, A (cid:54) = 0, a (cid:54) = 0, A (cid:54) = a ,and B (cid:54) = 0. We omit the derivation of the one-dimensional advection-diffusion equation (10),since the derivation is very similar to (and easier than) the two-dimensional advection-diffusionequation (39), which is given in Section 5 below. Remark 3.
In [15], the authors assumed that the non-negative kernel T ( a, ν, ν (cid:48) ) is the probabilitythat the bacteria changes the velocity from ν (cid:48) to ν , if a change of direction occurs. Therefore, ina one-dimensional space, T ( a, ν, ν (cid:48) ) = 1, and hence the right hand side of (8)-(9) for [15] has nofactor 1 /
2. In this work, we do not assume such an assumption; therefore T ( a, ν, ν (cid:48) ) = 1 /
2. Theassumption in [15] leads to the following equation instead of (10): ∂n∂t = ∂∂x (cid:18) ν α ∂n∂x − α B ν α ( A − α ) n (cid:19) . In [9, 12], Equation (10) is derived for a toy model that captures the essential excitation andadaptation components. Here, (10) can be used for any continuous tumbling function λ and alarger class of internal dynamics f ± = f ± νf , (see the following section for more details). E. coli bacteria
In what follows, we determine the terms in the advection-diffusion equation (10) for a population of
E. coli bacteria in a spatial domain [0 , L ] equipped with two chemical gradients MeAsp, denoted by8 ( x ), and serine, denoted by S ( x ). We further assume that S and S are respectively increasingand decreasing functions on [0 , L ], i.e., MeAsp accumulates near x = L and serine accumulatesnear x = 0. As we discussed in Section 2.1, in a one-dimensional space, the internal state of E. coli evolves according to the following ODE system: dadt = f ± ( a, S , S ) = f ( a, S , S ) ± νf ( a, S , S ) , (11)where, as described in (4), f ( a, S , S ) = ατ a N a ( a − a )( a − ,f ( a, S , S ) = N a ( a − (cid:18) γ γ S (cid:48) S + 11 + γ S (cid:48) S (cid:19) . (12)Here, S (cid:48) i = dS i /dx and γ := r /r denotes the ratio Tar/Tsr. Recall that r + r = 1, so indeed r = γ γ and r = γ . All the parameters used in this section are described in Section 2.1. Proposition 1.
Assume that the density functions p ± satisfy the decay condition p ± ( x, a, t ) ≤ C ( x, t ) e − c ( x,t ) a for some functions C, c : R × [0 , ∞ ) → R > and the stimuli S and S satisfy the shallow condition (cid:12)(cid:12)(cid:12)(cid:12) γ γ S (cid:48) ( x ) S ( x ) + 11 + γ S (cid:48) ( x ) S ( x ) (cid:12)(cid:12)(cid:12)(cid:12) ≤ min { q, − q } pν , ∀ x ∈ [0 , L ] , (13) where q = a and p = ατ a represent the adapted value and the the speed of adaptation, respectively.Then, for the given internal dynamics (11) , the dynamics of a population of E. coli, n ( x, t ) , canbe approximated by the advection-diffusion ∂n∂t = ∂∂x (cid:18) D ∂n∂x − χ (cid:18) γ γ S (cid:48) S + 11 + γ S (cid:48) S (cid:19) n (cid:19) , (14) where the diffusion coefficient D and the advection constant χ are as follows: D = ν λ + rq H , χ = rN Hq H ( q − ν ( λ + rq H )( N pq ( q − − λ − rq H ) . (15) Proof.
The proof is similar to the case of one stimulus, see [15], and the case of two-dimensionalspace which is given in Section 5 below.Note that the condition (13) holds if either the adaptation rate p is large or γ , S and S are chosenso that the left hand side (LHS) of (13) is small, i.e., the shallow condition is equivalent to eithersmall changes in the environment or fast adaptation. See the examples given in Section 4 for moredetails.Now we determine the boundary conditions of (14). Following the experimental set up in [4], wewant the population of the bacteria to be conserved in time, i.e., for any t ≥ ddt (cid:90) L n ( x, t ) dx = D (cid:16) ∂n∂x ( L, t ) − ∂n∂x (0 , t ) (cid:17) − χ (cid:0) V ( L ) n ( L, t ) − V (0) n (0 , t ) (cid:1) , (16)9here V ( x ) = γ γ S (cid:48) ( x ) S ( x ) + 11 + γ S (cid:48) ( x ) S ( x ) . (17)The following zero flux boundary conditions at x = 0 and x = L guarantee (16). For any t ≥ ∂n∂x (0 , t ) = χD V (0) n (0 , t ) and ∂n∂x ( L, t ) = χD V ( L ) n ( L, t ) . (18)In the following lemma, we provide sufficient conditions which guarantee existence and uniquenessof solutions of (14) with boundary conditions (18). Lemma 1.
Let V ( x ) be continuous on [0 , L ] and n ( x ) be a smooth non-negative function. Then, (14) with boundary condition (18) and initial condition n ( x ) admits a unique solution of the form n ( x, t ) = (cid:80) ∞ n =1 X n ( x ) T n ( t ) . Moreover, n ( x, t ) is uniformly bounded in x and t. This lemma can be proved by the method of separation of variables in a standard way: We canapply Sturm-Liouville theory [46] to solve the eigenproblem in which the first eigenvalue can bealso explicitly estimated to guarantee the uniform boundedness of the solution in time. For a proofsee Appendix A.1.
The bacterial responses to MeAsp and serine depend on the ratio of their chemoreceptors Tar andTsr, i.e., γ = Tar/Tsr. The goal is to find a positive γ ∗ and show that for γ > γ ∗ the bacteria tendto move toward a gradient of increasing MeAsp (i.e., accumulate near x = L ) and for γ < γ ∗ theymove toward a gradient of increasing serine (i.e., accumulate near x = 0).Let Φ( x ) be the steady state solution of the advection-diffusion equation (14) with boundarycondition (18). If S , S , and γ are chosen such that V ( x ) satisfies the condition in Lemma 1, thenthe solution of (14) converges to Φ( x ) as t → ∞ . Indeed, in the following examples, V ( x ) satisfiesthe condition in Lemma 1.Assuming that the bacteria start from a point x ∈ (0 , L ), they move toward a gradient of increasingMeAsp (respectively, serine) and accumulate near x = L (respectively, x = 0), if the steady statesolution of the advection-diffusion equation (14) admits a maximum on the right (respectively, left)sub-interval ( x , L ] (respectively, [0 , x )). Therefore, in what follows, we find conditions that Φ( x )admits a maximum on the right sub-interval ( x , L ] or the left sub-interval [0 , x ).To compute the steady state solution of (14), we let ∂n/∂t = 0, which gives a constant flux,i.e., D ∂n/∂x − χV ( x ) n = constant . Assuming zero flux boundary conditions (18), the constantbecomes zero and a simple calculation shows that the steady state solution satisfiesΦ( x ) = Φ( c ) exp (cid:26) χD (cid:90) xc V ( y ) dy (cid:27) . (19)10e choose c such that Φ( c ) >
0. Indeed, there is such a c by (16): ddt (cid:90) L n ( x, t ) dx = 0 ⇒ (cid:90) L n ( x, t ) dx = constant > ⇒ lim t →∞ (cid:90) L n ( x, t ) dx = (cid:90) L Φ( x ) dx = constant > ⇒ there exists c such that Φ( c ) > V as a function of both x and γ , V = V ( x, γ ). Considering the fact thatΦ (cid:48) ( x ) = χD V ( x, γ )Φ( x ) and Φ( x ) >
0, Φ takes a unique maximum at x ∗ ∈ [0 , L ] if, for any γ > V does not change sign or V is a non-increasing function of x and V ( x ∗ , γ ) = 0. Now we areready to find γ ∗ in the following lemma. Lemma 2.
Assume the bacteria start at x ∈ [0 , L ] and for any γ > , ∂V /∂x ≤ . Also, assumethat S and S are respectively increasing and decreasing functions on [0 , L ] . Then there exists γ ∗ > such that V ( x , γ ∗ ) = 0 and for γ > γ ∗ the bacteria accumulate on the right side of x andfor γ < γ ∗ they accumulate on the left side.Proof. A simple calculation shows that V ( x, γ ) = 0 if and only if γ ( x ) = S (cid:48) ( x ) /S ( x ) S (cid:48) ( x ) /S ( x ) . Let γ ∗ := γ ( x ) . For γ > γ ∗ , V ( x , γ ) >
0, therefore, since ∂V /∂x ≤
0, Φ takes its maximum(either x = L or x = x ∗ < L ) on the right side of x , and hence the bacteria accumulates towardthe right side of x . Similarly, if γ < γ ∗ , V ( x , γ ) <
0, and hence Φ takes its maximum (either x = 0 or x = x ∗ >
0) on the left side of x , and hence the bacteria accumulates toward the leftside of x .We refer to γ and γ ∗ as the bifurcation parameter and bifurcation value, respectively, since at γ = γ ∗ the direction of the bacterial changes. See Figure 2 below.Note that if the bacteria are initially distributed on [0 , L ] instead of locating on a single point x ,we consider γ ∗ = γ ∗ ( L/
2) as the bifurcation value.In the following section, we consider two sets of stimuli: (i) S linear and increasing, S linear anddecreasing; (ii) S exponential and increasing, S exponential and decreasing. We also assume thatthe bacteria are located at x = L/ V ( x, γ ) is a decreasing function on[0 , L ]. Hence, the conditions of Lemma 2 hold and, therefore, γ ∗ can be determined based on theinitial location of the bacteria, i.e., x = L/ To show that the advection-diffusion equation (14) with boundary condition (18) is a good approx-imation for the microscopic description of
E. coli chemotaxis, we run a Monte-Carlo agent-basedsimulation. A detailed description of the Monte-Carlo simulation is given in Appendix A.3.The following computational setting of our Monte-Carlo agent-based simulation is motivated bythe experimental set up in [4]. 11 patial Domain.
A one-dimensional channel of length of 400 µm ( x ∈ [0 , Stimuli.
Along the two sides of the channel two opposing chemical signals, S ( x ) and S ( x ), flowand diffuse across the channel. Two opposing linear and two opposing exponential chemical signalsare considered in Sections 4.1 and 4.2, respectively. Initial Condition. At t = 0 (sec), an ensemble of 100,000 agents is located in the center of thechannel ( x = 200). Boundary Conditions.
When a cell reaches a boundary, we relocate the cell to stay inside thedomain, i.e., zero flux boundary condition is applied.
Simulation Duration.
We simulate the bacterial behavior for 200 sec, t ∈ [0 , t = 200 . To illustrate distributions of the cells, we display histograms with 100 equal-sized bins.We use an explicit finite difference method to numerically solve the advection-diffusion equation(14) with the boundary condition (18).In the following examples, we compare the solutions of the macroscopic equation (14) with boundaryconditions (18) with results of the Monte-Carlo simulation. Further, for each case, we computethe bifurcation value γ ∗ defined in Section 3.2. To measure bacterial preference, we define thechemotactic migration coefficient (CMC):CMC x ( t ) = mean( x ( t )) − . (20)In the Monte-Carlo simulation, mean( x ( t )) is the average of individual positions x i at time t acrossthe channel, i.e., mean i ( x i ( t )). For a solution n ( x, t ) of (14), mean( x ( t )) is the expectation valueof the probability density n ( x, t ), i.e., (cid:82) L xn ( x, t ) dx . The absolute value of CMC x determines thedisplacement of the bacteria in x -direction. The sign of CMC x indicates their preference to theright or left. When CMC x > x < x = 200 (respectively, left, i.e., below x = 200). To demonstrate responses of
E. coli to two opposing linear gradients MeAsp and serine, andfollowing the experimental set up in [4], we let S ( x ) = 0 . x + 130 and S ( x ) = − . x + 20 (21)represent the concentrations of MeAsp and serine at each point x ∈ [0 , γ > V ( x, γ ) = γ γ . . x + 130 + 11 + γ − . − . x + 20 (22)is decreasing on [0 , V (200 , γ ) is an increasing function of γ , and V (200 , γ ∗ ) = 0 for γ ∗ ≈ . γ > .
985 (respectively, γ < . emark 4. In this example, for any x ∈ [0 , L ], S > S , | S (cid:48) | > | S (cid:48) | and | S (cid:48) /S | > | S (cid:48) /S | .Therefore, one may expect that the bacterial always choose to move towards MeAsp ( S ). However,as we proved in Lemma 2, when the ratio Tar/Tsr is small enough ( γ < γ ∗ ), the bacteria movetoward the gradient of serine. Figure 2 displays the relation between γ and the initial positionof the bacteria, x . The dotted curve γ ∗ ( x ) = (cid:0) S (cid:48) /S S (cid:48) /S (cid:1) ( x ) = x − x satisfying V ( x , γ ∗ ) = 0represents the bifurcation values in which the bacterial direction changes. As it is shown in Figure2, γ ∗ is an increasing function in x , that is, | S (cid:48) /S | increases faster than | S (cid:48) /S | as x increases.This means that if the bacteria start from near the right end point, a stronger force (a larger γ ∗ ) isneeded to drag them toward the gradient of serine ( S ). In the following section with exponentialgradients, although S > S and S (cid:48) > S (cid:48) everywhere, the needed force γ ∗ to drag the bacteria tothe gradient of serine is always equal to 1. The reason is that | S (cid:48) /S | ≡ | S (cid:48) /S | , in that case.Figure 2: Change of signs of V in (22)as x and γ vary. For ( x , γ ) in the darkred (respectively, blue) region, V be-comes positive (respectively, negative)as shown in the color bar. The dottedcurve is a set of ( x , γ ∗ ) where V = 0.The solid point at (200 , . To examine the result of Lemma 2, we choose two values for γ , γ = 1 . > γ ∗ ≈ .
985 and γ = 0 . < γ ∗ ≈ . E. coli obtained from the Monte-Carlo agent-based simulation and numerical solution of theadvection-diffusion (14). Three snapshots at times t = 10 , ,
200 (sec) are shown. As expected,the snapshots of a solution of (14) and the snapshots of a solution of Monte-Carlo simulation moveto the right when γ > γ ∗ , as shown in Figure 3(c), and they move to the left when γ < γ ∗ , asshown in Figure 3(a).Figures 3(b, d) display the corresponding CMC x which, as expected, is positive when γ > γ ∗ andthe bacteria accumulates on the right and is negative when γ < γ ∗ and the bacteria accumulateson the left.In Figure 3, the adaptation speed rate p is 0 . γ and p are chosen such that theshallow condition (13) holds. Therefore, by Proposition 1, the advection-diffusion equation (14)approximates the Fokker-Planck equations (8)-(9). A comparison of numerical solutions of (14)and the solutions of Monte-Carlo simulations in Figure 3 confirms this result. We now repeat the discussion of Section 4.1 for two opposing exponential gradients MeAsp andserine. We let S ( x ) = 130 e . x and S ( x ) = 8 e − . x − (23)13igure 3: (a) and (c): Comparisons of Monte-Carlo simulation and numerical solutions of (14)for two linear gradients (21) at times t = 10 , ,
200 (sec) with γ = 0 . < γ ∗ and γ = 1 . > γ ∗ ,in which the snapshots move to the left and right, respectively. (b) and (d): Comparisons of thecorresponding CMC x .represent the concentrations of MeAsp and serine at x ∈ [0 , γ ∗ , which determines the direction of bacteria, we apply Lemma 2.A simple calculation shows that V of this example is equal to V ( x, γ ) = 0 . γ − γ + 1 . For any γ > V ( x, γ ) is non-increasing on [0 , V (200 , γ ) is an increasing function of γ and V (200 , γ ∗ ) = 0 for γ ∗ = 1. Therefore, by Lemma 2, for γ > γ < γ , γ = 1 . > γ ∗ = 1 and γ = 0 . < γ ∗ =1. Figures 4(a, c) display distributions of the normalized density of E. coli obtained from theMonte-Carlo agent-based simulation and numerical solution of the advection-diffusion (14). Threesnapshots at times t = 10 , ,
200 (sec) are shown. As expected, the snapshots of a solution of(14) and the snapshots of a solution of Monte-Carlo simulation move to the right when γ > γ ∗ , asshown in Figures 4(c) and they move to the left when γ < γ ∗ , as shown in Figures 4(a).Figures 4(b, d) display the corresponding CMC x which, as expected, is positive when γ > γ ∗ andthe bacteria accumulates on the right and is negative when γ < γ ∗ and the bacteria accumulateson the left.Note that for the given exponential stimuli, the values of γ and p = 0 .
05 are chosen such that theshallow condition (13) holds. As discussed in Section 4.1, Proposition 1 and Figure 4 confirm thatthe numerical solutions of (14) agree well with the solutions of Monte-Carlo simulations.14igure 4: (a) and (c): Comparisons of Monte-Carlo simulation and numerical solutions of (14) fortwo exponential gradients (23) at times t = 10 , ,
200 (sec) with γ = 0 . < γ ∗ and γ = 1 . > γ ∗ ,in which the snapshots move to the left and right, respectively. (b) and (d): Comparisons of thecorresponding CMC x . In this section, we assume that the bacteria move in a two-dimensional space. Applying momentclosure techniques and parabolic scaling [9, 11, 12, 15], we derive an equation for the density of cellsat the population level that carries the description of an internal state of individuals in responseto the extracellular signals.As introduced in Section 2.2, let p ( x , a, ν, θ, t ) be a density function that describes a population ofagents at time t and location x = ( x, y ) (cid:62) with velocity ( ν , ν ) = ( ν cos θ, ν sin θ ) and an internalstate a. For the sake of simplicity, by fixing a constant speed ν , we let p θ ( x, y, a, t ) denotes thedensity of bacteria centered at ( x, y ) (cid:62) which move to the direction (cos( θ ) , sin( θ )) (cid:62) , θ ∈ [0 , π ),with the speed ν .According to the forward Fokker-Planck equation (7), for θ ∈ [0 , π ), p θ ( x, y, a, t ) satisfies ∂p θ ∂t + ∂∂x ( ν cos( θ ) p θ ) + ∂∂y ( ν sin( θ ) p θ ) + ∂∂a ( f θ ( a, S , S ) p θ )= 12 π λ ( a, S , S ) (cid:90) π ( p η ( x, y, a, t ) − p θ ( x, y, a, t )) dη, (24)where f θ and λ describe the internal dynamics and tumbling rate, respectively.In the presence of two extracellular signals S ( x, y, t ) and S ( x, y, t ), the evolution (6) of the internalstate of the bacteria that move to the direction (cos( θ ) , sin( θ )) (cid:62) with the speed ν is governed by15he following ordinary differential equation. dadt = f θ ( a, S , S )= f ( a, S , S ) + ν cos( θ ) f ( a, S , S ) + ν sin( θ ) f ( a, S , S ) , (25)where the real-valued functions f θ , f , f , and f are continuously differentiable. We assume that f , f and f have the Taylor expansions with respect to a as follows: f = A + A a + A a + · · · ,f = B + B a + B a + · · · ,f = B + B a + B a + · · · . Also, we assume that the tumbling rate λ = λ ( a, S , S ) has the Taylor expansion λ = α + α a + α a + · · · . All the Taylor constants are functions of S and S .At a fixed time t , consider a population of bacteria with internal dynamics (25) and tumbling rate λ that are located in ( x, y ). We want to show that, under some conditions, the population ofbacteria, which can be described by n ( x, y, t ) = (cid:90) R (cid:90) π p θ ( x, y, a, t ) dθda, solves an advection-diffusion equation: ∂n∂t = 12 ∂∂x (cid:18) ν α ∂n∂x − ν α B α ( A − α ) n (cid:19) + 12 ∂∂y (cid:18) ν α ∂n∂y − ν α B α ( A − α ) n (cid:19) . Following the techniques from [9] and [15], we define the fluxes as j (1) ( x, y, t ) = (cid:90) R (cid:90) π ν cos( θ ) p θ ( x, y, a, t ) dθda,j (2) ( x, y, t ) = (cid:90) R (cid:90) π ν sin( θ ) p θ ( x, y, a, t ) dθda, and the higher moments of the density and the fluxes as n i ( x, y, t ) = (cid:90) R (cid:90) π a i p θ ( x, y, a, t ) dθda, i = 1 , , . . . ,j (1) i ( x, y, t ) = (cid:90) R (cid:90) π a i ν cos( θ ) p θ ( x, y, a, t ) dθda, i = 1 , , . . . ,j (2) i ( x, y, t ) = (cid:90) R (cid:90) π a i ν sin( θ ) p θ ( x, y, a, t ) dθda, i = 1 , , . . . . (26) Assumption 1.
For any θ ∈ [0 , π ), the density functions p θ satisfy the decay condition p θ ( x, y, a, t ) ≤ C ( x, y, t ) e − c ( x,y,t ) a (27)for some functions C, c : R × [0 , ∞ ) → R > . Assumption 2.
For any i ≥ , we assume that n i , j (1) i , and j (2) i are negligible compared to n ,j (1)0 , j (2)0 , n , j (1)1 and j (2)1 . This assumption is made for the purpose of more tractable calculations.
Assumption 3. A = 0 , A (cid:54) = 0 , a (cid:54) = 0 , A (cid:54) = a , B (cid:54) = 0 , and B (cid:54) = 0 . This assumption guarantees the existence of unique solutions for some equations (see (32)-(33)below).Multiplying (24) by 1, ν cos( θ ), ν sin( θ ), and/or a and integrating the resulting equations withrespect to a and θ over R and [0 , π ), respectively, we obtain the following six equations: ∂n∂t + ∂j (1) ∂x + ∂j (2) ∂y = 0 ,∂j (1) ∂t + ν ∂n∂x = − α j (1) − α j (1)1 − (cid:88) k ≥ α k j (1) k ,∂j (2) ∂t + ν ∂n∂y = − α j (2) − α j (2)1 − (cid:88) k ≥ α k j (2) k ,∂n ∂t + ∂j (1)1 ∂x + ∂j (2)1 ∂y = A n + A n + B j (1) + B j (1)1 + B j (2) + B j (2)1 + (cid:88) k ≥ A k n k + (cid:88) k ≥ B k j (1) k + (cid:88) k ≥ B k j (2) k ,∂j (1)1 ∂t + ν ∂n ∂x = A j (1) + ( A − α ) j (1)1 + ν B n + ν B n + (cid:88) k ≥ ( A k − α k − ) j (1) k + ν (cid:88) k ≥ B k n k ,∂j (2)1 ∂t + ν ∂n ∂y = A j (2) + ( A − α ) j (2)1 + ν B n + ν B n + (cid:88) k ≥ ( A k − α k − ) j (2) k + ν (cid:88) k ≥ B k n k . (28)Here, we used the decaying condition (27) which, for i = 0 , , , . . . , yields (cid:90) R (cid:90) π a i ν cos( θ ) sin( θ ) p θ ( x, y, a, t ) dθda = 0 , (cid:90) R (cid:90) π a i ν cos ( θ ) p θ ( x, y, a, t ) dθda = (cid:90) R (cid:90) π a i ν sin ( θ ) p θ ( x, y, a, t ) dθda = ν n i ( x, y, a, t ) , where n = n, j (1)0 = j (1) and j (2)0 = j (2) .In what follows, we apply the parabolic scaling of space and time to the moment equation (26), toderive a set of non-dimensional equations. Let L, T, ν and N be scale factors for the length, time,velocity and the particle density, respectively. The parabolic scales of space and time are given byˆ x = (cid:16) εLν T (cid:17) xL , ˆ y = (cid:16) εLν T (cid:17) yL , ˆ t = ε tT , (29)17or arbitrary small ε >
0. Then, the dimensionless parameters are as follows.ˆ ν = νν , ˆ n = nN , ˆ j (1) = j (1) N ν , ˆ j (2) = j (2) N ν , ˆ n i = n i N , ˆ j i (1) = j (1) i N ν , ˆ j i (2) = j (2) i N ν , i = 1 , , . . . , ˆ α i = T α i , ˆ A i = T A i , ˆ B i = LB i , ˆ B i = LB i , i = 1 , , . . . . (30)Denoting ˆ w = (ˆ n, ˆ j (1) , ˆ j (2) , ˆ n , ˆ j (1)1 , ˆ j (2)1 ) (cid:62) , and by Assumption 2, we derive the following system of dimensionless moments from the dimen-sional equations (28): ε ∂ ˆ w∂ ˆ t + ε ∂∂ ˆ x P ˆ w + ε ∂∂ ˆ y P ˆ w = ε Q ˆ w + R ˆ w. (31)Here, the matrices P , P , Q and R are defined by partitioning into four 3 × P = (cid:32) P P (cid:33) , P = (cid:32) P P (cid:33) , Q = (cid:32) Q Q (cid:33) , R = (cid:32) R R S S + R (cid:33) , where is a zero matrix of dimension 3 and P = ˆ ν , P = ˆ ν , Q i = B i ˆ B i ˆ ν ˆ B i ˆ ν ˆ B i ,R i = − ˆ α i
00 0 − ˆ α i , S i = ˆ A i A i
00 0 ˆ A i , i = 0 , . To apply the regular perturbation method for w, we setˆ w = ˆ w + ε ˆ w + ε ˆ w + · · · , where ˆ w i = (cid:16) ˆ n i , ˆ j (1) i , ˆ j (2) i , ˆ n i , ˆ j (1) i , ˆ j (2) i (cid:17) (cid:62) , i = 0 , , , . . . . Substituting ˆ w into the dimensionless moment system (31) and collecting (cid:15) i terms, for i = 0 , , ε : R ˆ w = 0 (32) ε : R ˆ w = −Q ˆ w + ∂∂ ˆ x P ˆ w + ∂∂ ˆ y P ˆ w (33) ε : R ˆ w = −Q ˆ w + ∂∂ ˆ x P ˆ w + ∂∂ ˆ y P ˆ w + ∂∂ ˆ t ˆ w . (34)18y Assumption 3, (32) has a unique solution ˆ w of the formˆ w = (ˆ n , , , , , (cid:62) , where ˆ n is nonzero. The second equation (33) yields − ˆ α ˆ j (1)1 − ˆ α ˆ j (1)11 − ˆ α ˆ j (2)1 − ˆ α ˆ j (2)11 ˆ A ˆ n ( ˆ A − ˆ α )ˆ j (1)11 + ˆ ν ˆ B ˆ n ( ˆ A − ˆ α )ˆ j (2)11 + ˆ ν ˆ B ˆ n = ˆ ν ∂ ˆ n ∂ ˆ x ˆ ν ∂ ˆ n ∂ ˆ y ; (35)and from the last two equalities of (35), it followsˆ j (1)11 = − ˆ ν ˆ B
2( ˆ A − ˆ α ) ˆ n and ˆ j (2)11 = − ˆ ν ˆ B
2( ˆ A − ˆ α ) ˆ n . Moreover, plugging ˆ j (1)11 and ˆ j (2)12 into the second and third equalities in (35), we obtainˆ j (1)11 = − ˆ ν α ∂ ˆ n ∂ ˆ x + ˆ ν ˆ α ˆ B α ( ˆ A − ˆ α ) and ˆ j (2)11 = − ˆ ν α ∂ ˆ n ∂ ˆ y + ˆ ν ˆ α ˆ B α ( ˆ A − ˆ α ) . (36)Noticing that the right hand side of (34) is in the image of R and (1 , , , , , (cid:62) is in the kernel of R , the right hand side of (34) must be orthogonal to (1 , , , , , (cid:62) by the Fredholm AlternativeTheorem, which yields ∂∂ ˆ t ˆ n + ∂∂ ˆ x ˆ j (1)1 + ∂∂ ˆ y ˆ j (2)1 = 0 . (37)By substituting the results in (36) into (37), we obtain the following equation for ˆ n : ∂ ˆ n ∂ ˆ t = 12 ∂∂ ˆ x (cid:32) ˆ ν ˆ α ∂ ˆ n ∂ ˆ x − ˆ ν ˆ α ˆ B ˆ α ( ˆ A − ˆ α ) ˆ n (cid:33) + 12 ∂∂ ˆ y (cid:32) ˆ ν ˆ α ∂ ˆ n ∂ ˆ y − ˆ ν ˆ α ˆ B ˆ α ( ˆ A − ˆ α ) ˆ n (cid:33) . (38)Similarly, we can derive the evolution equation for ˆ n which solves (38).For n ( x, y, t ) = n ( x, y, t ) + εn ( x, y, t ) + O ( ε ) , if the terms in O ( ε ) are ignored, (38) for theoriginal (dimensional) variable n is transformed into ∂n∂t = 12 ∂∂x (cid:18) ν α ∂n∂x − ν α B α ( A − α ) n (cid:19) + 12 ∂∂y (cid:18) ν α ∂n∂y − ν α B α ( A − α ) n (cid:19) . (39)For the spatial domain [0 , L ] × [0 , L ], assuming that the population of bacteria is conserved intime and there is no flux along the boundary, we can impose the following boundary conditions: D ∂n∂x (0 , y, t ) = χ (0 , y ) n (0 , y, t ) and D ∂n∂x ( L , y, t ) = χ ( L , y ) n ( L , y, t ) D ∂n∂y ( x, , t ) = χ ( x, n ( x, , t ) and D ∂n∂y ( x, L , t ) = χ ( x, L ) n ( x, L , t ) , (40)where D = ν α and for i = 1 , , χ i ( x, y ) = ν α B i α ( A − α ) . .1 Application to a population of E. coli bacteria
We now compute the coefficients of the macroscopic equation (39) for
E. coli bacteria. As wediscussed in Section 2.1, in a two-dimensional space, the internal state of
E. coli evolves accordingto the following ODE system: dadt = f ( a, S , S ) + ν cos( θ ) f ( a, S , S ) + ν sin( θ ) f ( a, S , S ) , where f ( a, S , S ) = ατ a N a ( a − a )( a − ,f ( a, S , S ) = N a ( a − (cid:16) γ γ ∂ x S S + 11 + γ ∂ x S S (cid:17) ,f ( a, S , S ) = N a ( a − (cid:16) γ γ ∂ y S S + 11 + γ ∂ y S S (cid:17) . The constant terms of the Taylor expansions of f and f are zero. However, in Assumption 3, wesaw that these constant terms must be non-zero. To fix this issue, we make a change of coordinate,ˆ a = a − a , and obtain the following new internal dynamics of ˆ a : d ˆ adt = ˆ f (ˆ a, S , S ) + ν cos( θ ) ˆ f (ˆ a, S , S ) + ν sin( θ ) ˆ f (ˆ a, S , S ) , (41)where ˆ f (ˆ a, S , S ) = pN ˆ a (ˆ a + q )(ˆ a + q − , ˆ f (ˆ a, S , S ) = N (ˆ a + q )(ˆ a + q − (cid:16) γ γ ∂ x S S + 11 + γ ∂ x S S (cid:17) , ˆ f (ˆ a, S , S ) = N (ˆ a + q )(ˆ a + q − (cid:16) γ γ ∂ y S S + 11 + γ ∂ y S S (cid:17) . We let q = a and p = ατ a (42)represent the adapted value and the the speed of adaptation, respectively.Next, we transform the tumbling rate, discussed in (5), into the new coordinate ˆ a as follows: λ (ˆ a ) = λ + r (ˆ a + q ) H , where r = 1 τ a H . (43)All the model parameters N, p, q, r, and H are assumed to be positive constants and are as givenin Table 1.Let Assumption 1 hold. In the following two lemmas, we provide sufficient conditions that lead toAssumption 2. Lemma 3 (Shallow condition) . Let c = min { q, − q } . If for any ( x, y ) ∈ [0 , L ] × [0 , L ] , and any θ ∈ [0 , π ) (cid:12)(cid:12)(cid:12) cos( θ ) (cid:16) γ γ ∂ x S S + 11 + γ ∂ x S S (cid:17) + sin( θ ) (cid:16)
11 + γ ∂ y S S + 11 + γ ∂ y S S (cid:17)(cid:12)(cid:12)(cid:12) ≤ cpν , (44) and | ˆ a (0) | ≤ c , then | ˆ a ( t ) | ≤ c for all t ≥ . roof. To show | ˆ a ( t ) | ≤ c , it suffices to show that d ˆ a/dt < >
0) at ˆ a = c (respec-tively, ˆ a = − c ). Using (44) and c < − q into (41), we obtain the desired result.Note that the inequality (44) holds if either the adaptation rate p is sufficiently large or γ , S , and S are chosen so that the LHS of (44) is sufficiently small. See the examples given in Section 6 formore details. Lemma 4.
Let
L, T, ν , and N be scale factors for the length, time, velocity, and cell density,respectively, as introduced in (29) and (30) . We define dimensionless quantities as follows: (cid:92) (cid:16) ∇ x S i S i (cid:17) = ν ε ∇ x S i S i , ( i = 1 , , ˆ N = T N, ˆ γ = γ, ˆ p = p, ˆ q = q and ˆ r = T r.
Then, under the shallow condition (44) , for any i ≥ , ˆ j (1) i ˆ n ≤ C (1) i ε i , ˆ j (2) i ˆ n ≤ C (2) i ε i , and ˆ n i ˆ n ≤ D i ε i , for some constants C (1) i = O (1) , C (2) i = O (1) , and D i = O (1) . The proof of Lemma 4 can be completed as proved in [15]; thus we omit the proof.The shallow condition (44) guarantees that the higher moments n i , j (1) i , and j (2) i , i ≥
2, are of order ε , O ( ε ). Indeed, we can close the moment equations (28) by considering the higher moments asthe error terms of O ( ε ) . Simple calculations show that the Taylor coefficients of ˆ f , ˆ f , and ˆ f are A = 0 , A = N pq ( q − , B = N q ( q − (cid:16) γ γ ∂ x S S + 11 + γ ∂ x S S (cid:17) ,B = N q ( q − (cid:16) γ γ ∂ y S S + 11 + γ ∂ y S S (cid:17) , α = λ + rq H , α = rHq H q , that satisfy Assumption 3. Then, with Assumption 1, the shallow condition for the stimuli, andthe internal dynamics (41), a population of E. coli , n ( x, y, t ), solves the following equation ∂n∂t = ∇ · (cid:16) D ∇ n − χ (cid:16) γ γ ∇ x S S + 11 + γ ∇ x S S (cid:17) n (cid:17) , (45)where the diffusion coefficient D and the advection constant χ are D = ν λ + rq H ) > , χ = rN Hq H ( q − ν λ + rq H )( N pq ( q − − λ − rq H ) > . For the spatial domain [0 , L ] × [0 , L ], the boundary conditions (40) becomes D ∂n∂x (0 , y, t ) = χV (0 , y ) n (0 , y, t ) and D ∂n∂x ( L , y, t ) = χV ( L , y ) n ( L , y, t ) ,D ∂n∂y ( x, , t ) = χV ( x, n ( x, , t ) and D ∂n∂y ( x, L , t ) = χV ( x, L ) n ( x, L , t ) , (46)21here V ( x, y ) = γ γ ∂ x S S + 11 + γ ∂ x S S and V ( x, y ) = γ γ ∂ y S S + 11 + γ ∂ y S S . In the following lemma, we present sufficient conditions that guarantee the existence and uniquenessof solutions of (45) with the boundary condition (46). The proof is followed by Lemma 1 due tothe assumptions on V and V . Lemma 5.
Let V ( x, y ) and V ( x, y ) be continuous on Ω := [0 , L ] × [0 , L ] , and n ( x, y ) bea smooth non-negative function. If V ( x, y ) = V ( x ) and V ( x, y ) = V ( y ) , then (45) with theboundary condition (46) and the initial condition n ( x, y,
0) = n ( x, y ) admits a unique solution n ( x, y, t ) in the form of (cid:80) ∞ n =1 X n ( x ) Y n ( y ) T n ( t ) . Moreover, n ( x, y, t ) is uniformly bounded in x, y and t. In a similar way to explaining the direction of bacterial migration in Section 3.2, we exploreproperties of the steady state solution of the advection-diffusion equation (45) and predict thedirection of bacteria. To do this, we choose S , S and γ so that V ( x, y ) and V ( x, y ) satisfy theconditions given in Lemma 5.To compute the steady state solution of the advection-diffusion equation (45) with zero flux bound-ary conditions (46), we let the flux at x direction and the flux at y direction be zero, i.e., J x ( x, y ) := D ∂n∂x − χV ( x, y ) n = 0 ,J y ( x, y ) := D ∂n∂y − χV ( x, y ) n = 0 , which yield (cid:32) ∂∂x log n ∂∂y log n (cid:33) = χD (cid:32) V ( x, y ) V ( x, y ) (cid:33) . (47)Note that this equation cannot be satisfied for any arbitrary V and V . Since the LHS is a gradient,a necessary and sufficient condition for the equation to hold is ∂V ∂y = ∂V ∂x . (48)Note that if V and V satisfy the conditions of Lemma 5, they automatically satisfy (48). Underthis condition, the steady state solution can be obtained by simple integration of (47):Φ( x, y ) = Φ( c , c ) exp (cid:26) χD (cid:18)(cid:90) xc V ( z, y ) dz + (cid:90) yc V ( c , z ) dz (cid:19)(cid:27) , (49)where ( c , c ) ∈ [0 , L ] × [0 , L ] are chosen such that Φ( c , c ) is a positive constant. Similar towhat we discussed in Section 3.2, if ∂V /∂x ≤ ∂V /∂y ≤
0, then the signs of V and V
22t the initial point ( x , y ) can determine the direction of the motion of bacteria. We let γ ∗ bethe bifurcation value that determines the right/left direction (i.e., V ( x , y , γ ∗ ) = 0) and γ ∗ bethe bifurcation value that determines the up/down direction (i.e., V ( x , y , γ ∗ ) = 0). Then, threescenarios are possible: (i) for max { γ ∗ , γ ∗ } < γ , the bacteria move to the northeast and accumulatein A := { x < x < L , y < y < L } ; (ii) for min { γ ∗ , γ ∗ } < γ < max { γ ∗ , γ ∗ } the bacteriaeither move to the southeast and accumulate in A := { x < x < L , < y < y } or move tothe northwest and accumulate in A := { < x < x , y < y < L } ; (iii) for γ < min { γ ∗ , γ ∗ } , thebacteria move to the southwest and accumulate in A := { < x < x , < y < y } .In the following section, we consider three sets of stimuli, which their corresponding V and V satisfy the condition of Lemma 5 (and hence (48)). For each set we find the bifurcation valueswhich determine the direction of bacteria. To validate our two-dimensional macroscopic approximation (45), we run a Monte-Carlo simulationfor microscopic equation (24). Our numerical experimental set up is very similar to that of Section4, which we generalize to a two-dimensional space as follows. Note that since this work is motivatedby [4], we choose a computational setting to be qualitatively similar to the experimental set up of[4] as well.
Spatial Domain.
A channel of area of 400 µm by 1600 µm ( x ∈ [0 , , y ∈ [0 , Stimuli.
Along the two sides of the channel x = 0 and x = 400, two opposing chemical signals S ( x, y ) and S ( x, y ), which respectively represent the concentrations of MeAsp and serin at ( x, y ),flow and diffuse across the channel. Three sets of stimuli will be considered in Sections 6.1– 6.3,below. Initial Condition. At t = 0 (sec), an ensemble of 100,000 agents is located in the center of thechannel ( x = 200 and y = 800). Boundary Condition.
We assume zero flux condition on the boundaries, i.e., when a cell reachesa boundary, we relocate the cell to stay inside the domain.
Simulation Duration.
We simulate the bacterial behavior for t ∈ [0 , t = 200.The distributions of the solutions are displayed by using histograms with 2500 equal-sized bins. Tosolve the advection-diffusion equation (45) with boundary conditions (46), we use an explicit finitedifference method. The summary of input data is given in Table 2 (see Appendix A.3), and moredetails can be also found in Section 4.In what follows, we show some numerical results for three different choices of the stimuli com-binations: Linear–Linear in Section 6.1, Exponential–Exponential in Section 6.2, and Linear × Exponential–Linear × Exponential in Section 6.3. We will show that (i) for some γ ∗ , when γ > γ ∗ ,the bacteria move to the the gradient of increasing MeAsp and when γ < γ ∗ , the bacteria moveto the gradient of increasing serine; and (ii) under the condition of Lemma 3, the Monte-Carloagent-based simulations and the numerical solutions of (45) agree well.23 .1 Chemotaxis in response to two linear gradients Let S ( x, y ) = 0 . x + 130 and S ( x, y ) = − . x + 20 be two opposing linear gradients for MeAspand serine, respectively. Note that the stimuli are constant with respect to y . In this case, V ( x, y ) = V ( x ), as defined in Section 4.1, and V ( x, y ) = 0. Therefore, the condition of Lemma5 and hence (48) hold and the bacteria only move to the right or left (no up or down movement).Furthermore, the bifurcation value is equal to γ ∗ ≈ . γ = 1 . γ = 0 . t = 0 (left), t = 60 (middle), and t = 200 (right). Figure 5(c)(respectively,Figure 6(c)) displays the corresponding CMCs in x − direction and y − direction. (a) Monte-Carlo simulation for γ = 1 . p = 1(b) Numerical solutions of (45) for γ = 1 . p = 1(c) CMC x (left) and CMC y (right) Figure 5: (a) and (b): Comparisons of the Monte-Carlo simulations and numerical solutions of(45) in response to two linear gradients, when γ = 1 .
5. In this case the bacteria move to the right,the gradient of increasing MeAsp. (c): Comparisons of the corresponding CMCs.24 a) Monte-Carlo simulation for γ = 0 . p = 1(b) Numerical solutions of (45) for γ = 0 . p = 1(c) CMC x (left) and CMC y (right) Figure 6: (a) and (b): Comparisons of the Monte-Carlo simulations and numerical solutions of (45)for two linear gradients in (21), γ = 0 .
5, and p = 1. In this case the bacteria move to the left, thegradient of increasing serine. Plots in (a) and (b) are displayed only for ( x, y ) ∈ [0 , × [600 , γ ∗ ≈ .
985 in Section4.1, these figures confirm that the chemotactic preference of bacteria depends on the relativeabundances of receptors, i.e., when γ = 1 . > γ ∗ , the bacteria move to the gradient of increasingMeAsp (CMC x > γ = 0 . < γ ∗ , the bacteria move to the gradient ofincreasing serine (CMC x < S and S are independent of y , the bacteria move in the y -direction very slightly, as evidencedby CMC y ≈
0. Thus, although we run all the simulations on the domain [0 , × [0 , , × [600 , .2 Chemotaxis in response to two exponential gradients We assume that bacteria are exposed to two opposing exponential gradients S ( x, y ) = 130 e . x and S ( x, y ) = 8 e − . x − . In this case, V ( x, y ) = V ( x ), as defined in Section 4.2, and V ( x, y ) = 0. Therefore, condition (48)holds and the bacteria only move to the right or left (no up or down movement). Furthermore, thebifurcation value is equal to γ ∗ = 1, as computed in Section 4.2.In Figures 7 and 8, we compare the results of the Monte-Carlo simulation with numerical solutionof (45) and their corresponding CMCs. From the plots, we can see that (45) captures the behaviorof individuals well. Recalling the bifurcation value γ ∗ = 1 of the ratio of Tar to Tsr in Section 4.2,as expected, the individuals travel to the right when γ = 1 . > γ ∗ as in Figure 7 and move to theleft when γ = 0 . < γ ∗ as in Figure 8. (a) Monte-Carlo simulation for γ = 1 . p = 0 . γ = 1 . p = 0 . x (left) and CMC y (right) Figure 7: (a) and (b): Comparisons of Monte-Carlo simulation and numerical solutions of (45) inresponse to two exponential gradients when γ = 1 .
1. Plots in (a) and (b) are displayed only for( x, y ) ∈ [0 , × [600 , a) Monte-Carlo simulation for γ = 0 . p = 0 . γ = 0 . p = 0 . x (left) and CMC y (right) Figure 8: (a) and (b): Comparisons of Monte-Carlo simulation and numerical solutions of (45) inresponse to two exponential gradients when γ = 0 .
9. Plots in (a) and (b) are displayed only for( x, y ) ∈ [0 , × [600 , In Sections 6.1 and 6.2, we used two opposing gradients, independent of y , to reproduce chemotaxisexperiments in the literature. In what follows, we assume that two opposing gradients MeAsp ( S )and serine ( S ) satisfy S ( x, y ) = (0 . x + 130) e . y − and S ( x, y ) = ( − . x + 20) e − . y − . (50)Note that each gradient increases toward the corners (0 ,
0) and (400 , V ( x, y ) = V ( x ), as defined in Section 4.1,and V ( x, y ) = 0 . γ − γ +1 . Therefore, condition (48) holds. Furthermore, the bifurcation valuesare γ ∗ ≈ . γ ∗ = 1. Therefore, three scenarios occur: (i) for27 > . < γ < γ < .
985 the bacteria move to southwest. As expected, the plots in Figure 9 showthat bacteria accumulate toward the corner (400 , γ = 1 . >
1. Also, the solution of(45) agrees well with the result of the Monte-Carlo simulation. (a) Monte-Carlo simulation for γ = 1 . p = 1 (b) Numerical solutions of (14) for γ = 1 . p = 1(c) CMC x (left) and CMC y (right) Figure 9: (a) and (b): Comparisons of Monte-Carlo simulation and numerical solutions of (45) forgradients (50) for γ = 1 .
5. (c): Comparisons of the corresponding CMCs.
In this work, we studied the movement of a population of
E. coli bacteria in response to two stimuliin a one- and a two-dimensional environment. Experimental results [4] show that the bacterialchemotactic preference to serine and MeAsp depends on the ratio of their chemoreceptors, namely γ = Tar / Tsr. In a shallow-gradient regime, we analytically found a threshold γ ∗ that determinesthe bacterial preference, i.e., if γ > γ ∗ , the bacteria move toward the gradient of MeAsp, and if γ < γ ∗ , the bacteria move toward the gradient of serine. We examined our results in an environmentwhere one stimulus is dominant everywhere and observed that in such a situation, a bigger force γ ∗ might be needed to change the preference of the bacteria.We started with a microscopic model for a population of bacteria carrying a one-dimensionalinternal dynamics. Indeed the microscopic equation is the forward Fokker-Planck equation of astochastic model which describes bacterial chemotaxis [47]. Then, we approximated the microscopicFokker-Planck equation by a macroscopic advection-diffusion equation which is more tractablemathematically. We compared the numerical solution of the advection-diffusion equation with a28onte-Carlo simulation of the bacterial chemotaxis to validate the approximation in a shallow-gradient regime.In [48], the authors found that E. coli cells respond to the gradient of chemoattractant not only bybiasing their own random-walk swimming pattern through the intracellular pathway, but also byactively secreting a chemical signal into the extracellular medium, possibly through a communica-tion signal transduction pathway. The extracellular signaling molecule is a strong chemoattractantthat attracts distant cells to the food source. They showed that cell-cell communication enhancesbacterial chemotaxis toward external attractants. Incorporating such chemoattractant into micro-scopic model is one of the main areas of our future investigation. This cell-cell communication canbe modeled as an external force to each cell and described by an extra term into the LHS of (7).In this work, we only considered a one-dimensional internal dynamics. To obtain the internaldynamics of
E. coli in response to multiple stimuli, we applied the heterogeneous MWC model (1)[17, 37, 49], which can capture the total activity level of bacterium affected by the stimuli andmathematically is tractable. In this model, all receptors within the cluster are assumed to turn onand off simultaneously, and therefore, only the total kinase activity and total methylation level areconsidered. However, in a mixed-receptor cluster, it was found that receptor methylation dynamicsis ligand specific. Hence, a local adaptation model, such as the Ising-type model, can better explainthe adaptation dynamics of the mixed-receptor cluster, see e.g., [49] and [50]. Such models requirehigher dimensional equations to describe the internal dynamics. In our future works, we generalizeour result to two-dimensional internal dynamics and for each receptor Tar and Tsr, we will considerseparate activity levels a and a instead of a in (1) and separate methylation dynamics dm /dt and dm /dt instead of (2).In this work, we showed that the advection-diffusion approximation is valid under the shallow-gradient condition. However, we numerically observed that even if the shallow-gradient conditiondoes not hold, some of our results remain valid. For example, Figure 10 shows that under thecondition of Section 4.1, the behavior of the bacteria does not change even when the adaptationrate p does not satisfy the sallow-gradient condition (gray region). We also observed that p does notaffect the preference of bacteria. Further future directions include consideration of a more generalclass of stimuli by relaxing the shallow-gradient condition as in [51, 52] and allowing time-varyingstimuli as in [7, 22]. Figure 10: CMC x of the steady state(19) for S and S in (21) for x = 200.For ( γ, p ) in the dark red (respectively,blue) region, CMC x becomes positive(respectively, negative) as shown in thecolor bar. For ( γ, p ) in the dark grey re-gion, the shallow condition (13) is notsatisfied. The dotted line represents γ ≈ . Acknowledgement
The authors would like to thank Professor Eduardo Sontag for sharing the Matlab codes for one-dimensional space (used in [15]) and Professor Hans Othmer for helpful discussions. This work ispartially supported by the University of Iowa Old Gold Fellowship.
References [1] H. C. Berg and D. A. Brown. Chemotaxis in
Escherichia coli analysed by three-dimensionaltracking.
Nature , 239(5374):500-504, 1972.[2] R. M. Macnab and D. E. Koshland. The gradient-sensing mechanism in bacterial chemotaxis.
Proc. Natl. Acad. Sci. , 69(9):2509-2512, 1972.[3] N. Vladimirov and V. Sourjik. Chemotaxis: how bacteria use memory.
Biological chemistry ,390(11):1097–1104, 2009.[4] Y. Kalinin, S. Neumann, V. Sourjik, and M. Wu. Responses of
Escherichia coli bacteria totwo opposing chemoattractant gradients depend on the chemoreceptor ratio.
J. Bacteriol. ,192(7):1796-1800, 2010.[5] M. J. Tindall, P. K. Maini, S. L. Porter, and J. P. Armitage. Overview of mathematicalapproaches used to model bacterial chemotaxis II: bacterial populations.
Bull Math Biol ,70(6):1570, 2008.[6] W. Alt. Biased random walk models for chemotaxis and related diffusion approximations.
J.Math. Biol. , 9(2):147-177, 1980.[7] Y. Tu, T. S. Shimizu, and H. C. Berg. Modeling the chemotactic response of
Escherichia coli to time-varying stimuli.
Proc. Natl. Acad. Sci. , 105(39):14855-14860,2008.[8] M. P. Edgington and M. J. Tindall. Mathematical Analysis of the
Escherichia coli
ChemotaxisSignalling Pathway.
Bull Math Biol , 80(4):758-787, 2018.[9] R. Erban and H. G. Othmer. From individual to collective behavior in bacterial chemotaxis.
SIAM J. Appl. Math. , 65(2):361-391, 2004.[10] E. F. Keller and L. A. Segel. Model for chemotaxis.
J. Theor. Biol. , 30(2):225-234,1971.[11] R. Erban and H. G. Othmer. From signal transduction to spatial pattern formation in
E. coli :a paradigm for multiscale modeling in biology.
Multiscale Model. Simul. , 3(2):362-394, 2005.[12] C. Xue and H. G. Othmer. Multiscale models of taxis-driven patterning in bacterial popula-tions.
SIAM J. Appl. Math. , 70(1):133-167, 2009.[13] H. G. Othmer and A. Stevens. Aggregation, blowup and collapse: the ABCs of generalizedtaxis,
SIAM J. Appl. Math , 57(4):1044-1081, 1997.[14] K. J. Painter, P. K. Maini, H. G. Othmer. Development and applications of a model for cellularresponse to multiple chemotactic cues.
J. Math. Biol. , 41(4):285–314, 2000.3015] Z. Aminzare and E. D. Sontag. Remarks on a population-level model of chemotaxis: advection-diffusion approximation and simulations. arXiv preprint arXiv:1302.2605 , 2013.[16] F. Menolascina, R. Rusconi, V. I. Fernandez, S. Smriga, Z. Aminzare, E. D. Sontag, and R.Stocker. Logarithmic sensing in
Bacillus subtilis aerotaxis.
NPJ Syst Biol Appl , 3:16036, 2017.[17] B. Hu and Y. Tu. Behaviors and strategies of bacterial navigation in chemical and nonchemicalgradients.
PLoS Comput. Biol. , 10(6):e1003672, 2014.[18] H. Salman and A. Libchaber. A concentration-dependent switch in the bacterial response totemperature.
Nature cell biology , 9(9):1098, 2007.[19] M. Demir, C. Douarche, A. Yoney, A. Libchaber, and H. Salman. Effects of population densityand chemical environment on the behavior of
Escherichia coli in shallow temperature gradients.
Physical Biology , 8(6):063001, 2011.[20] Y. Yang and V. Sourjik. Opposite responses by different chemoreceptors set a tunable prefer-ence point in
Escherichia coli pH taxis.
Molecular Microbiology , 86(6):1482-1489, 2012.[21] Y. V. Kalinin, L. Jiang, Y. Tu, and M. Wu. Logarithmic sensing in
Escherichia coli bacterialchemotaxis.
Biophysical Journal , 96(6):2439-2448, 2009.[22] L. Jiang, Q. Ouyang, and Y. Tu. Quantitative modeling of
Escherichia coli chemotactic motionin environments varying in space and time.
PLoS Comput. Biol. , 6(4):e1000735, 2010.[23] H. G. Othmer, S. R. Dunbar, and W. Alt. Models of dispersal in biological systems.
J. ofMath. Biol. , 26:263-298, 1988.[24] H. C. Berg and L. Turner. Chemotaxis of bacteria in glass capillary arrays.
Escherichia coli ,motility, microchannel plate, and light scattering.
Biophysical Journal , 58(4):919-930, 1990.[25] G. H. Wadhams and J. P. Armitage. Making sense of it all: bacterial chemotaxis.
Nat. Rev.Mol. Cell Biol. , 5(12):1024-1037, 2004.[26] M. Welch, K. Oosawa, S.-L. Aizawa, and M. Eisenbach. Phosphorylation-dependent bindingof a signal molecule to the flagellar switch of bacteria.
Proc. Natl. Acad. Sci. , 90(19):8787-8791,1993.[27] A. Bren, M. Welch, Y. Blat, and M. Eisenbach. Signal termination in bacterial chemotaxis:CheZ mediates dephosphorylation of free rather than switch-bound CheY.
Proc. Natl. Acad.Sci. , 93(19):10090-10093, 1996.[28] K. Lipkow, S. S. Andrews, and D. Bray. Simulated diffusion of phosphorylated CheY throughthe cytoplasm of
Escherichia coli . J. Bacteriol. , 187(1):45-53, 2005.[29] K. Lipkow. Changing cellular location of CheZ predicted by molecular simulations.
PLoSComput. Biol. , 2(4), 2006.[30] W. R. Springer and D. E. Koshland. Identification of a protein methyltransferase as the cheRgene product in the bacterial sensing system.
Proc. Natl. Acad. Sci. , 74(2):533-537, 1977.3131] J. B. Stock and D. E. Koshland. A protein methylesterase involved in bacterial sensing.
Proc.Natl. Acad. Sci. , 75(8):3659-3663, 1978.[32] D. Bray and R. B. Bourret. Computer analysis of the binding reactions leading to a trans-membrane receptor-linked multiprotein complex involved in bacterial chemotaxis.
Mol. Biol.Cell , 6(10):1367-1380, 1995.[33] T. C. Terwilliger, J. Y. Wang, and D. E. Koshland. Kinetics of receptor modification. Themultiplymethylated aspartate receptors involved in bacterial chemotaxis.
Journal of BiologicalChemistry , 261(23):10814-10820, 1986.[34] S. A. Simms, A. M. Stock, and J. B. Stock. Purification and characterization of the s-adenosylmethionine: glutamyl methyltransferase that modifies membrane chemoreceptor pro-teins in bacteria.
Journal of Biological Chemistry , 262(18):8537-8543, 1987.[35] N. Vladimirov, L. Løvdok, D. Lebiedz, and V. Sourjik. Dependence of bacterial chemotaxison gradient shape and adaptation rate.
PLoS Comput. Biol. , 4(12):e1000242, 2008.[36] D. Clausznitzer, O. Oleksiuk, L. Løvdok, V. Sourjik, and R. G. Endres. Chemotactic responseand adaptation dynamics in
Escherichia coli . PLoS Comput. Biol. , 6(5), 2010.[37] B. A. Mello and Y. Tu. An allosteric model for heterogeneous receptor complexes:under-standing bacterial chemotaxis responses to multiple stimuli.
Proc. Natl. Acad. Sci. ,102(48):17354-17359, 2005.[38] S. Neumann, C. H. Hansen, N. S. Wingreen, and V. Sourjik. Differences in signalling by directlyand indirectly binding ligands in bacterial chemotaxis.
The EMBO Journal , 29(20):3484-3495,2010.[39] J. Monod, J. Wyman, and J.-P. Changeux. On the nature of allosteric transitions: a plausiblemodel.
J. Mol. Biol. , 12(1):88-118, 1965.[40] V. Sourjik and H. C. Berg. Receptor sensitivity in bacterial chemotaxis.
Proc. Natl. Acad.Sci. , 99(1):123-127, 2002.[41] T. S. Shimizu, N. Delalez, K. Pichler, and H. C. Berg. Monitoring bacterial chemo-taxis byusing bioluminescence resonance energy transfer: absence of feedback from the flagellar motors.
Proc. Natl. Acad. Sci. , 103(7):2093-2097, 2006.[42] B. A. Mello and Y. Tu. Effects of adaptation in maintaining high sensitivity over a wide rangeof backgrounds for
Escherichia coli chemotaxis.
Biophysical Journal , 92(7):2329-2337, 2007.[43] O. Shoval, L. Goentoro, Y. Hart, A. Mayo, E. Sontag, and U. Alon. Fold-change detection andscalar symmetry of sensory input fields.
Proc. Natl. Acad. Sci. , 107(36):15995-16000, 2010.[44] M. D. Lazova, T. Ahmed, D. Bellomo, R. Stocker, and T. S. Shimizu. Response rescaling inbacterial chemotaxis.
Proc. Natl. Acad. Sci. , 108(33):13870-13875, 2011.[45] O. Shoval, U. Alon, and E. Sontag. Symmetry invariance for adapting biological systems.
SIAM J. Appl. Math. , 10(3):857-886, 2011. 3246] T. Kapitula and K. Promislow.
Spectral and dynamical stability of nonlinear waves . Springer,2013.[47] D. W. Stroock. Some stochastic processes which arise from a model of the motion of a bac-terium.
Z. Wahrscheinlichkeitstheor. verw. Geb. , 28(4):305-315, 1974.[48] Z. Long, B. Quaife, H. Salman, and Z. N. Oltvai. Cell-cell communication enhances bacterialchemotaxis toward external attractants.
Scientific reports , 7(1):12855, 2017.[49] J. E. Keymer, R. G. Endres, M. Skoge, Y. Meir, and N. S. Wingreen. Chemosensing in
Escherichia coli : two regimes of two-state receptors.
Proc. Natl. Acad. Sci. , 103(6):1786-1791,2006.[50] B. Hu and Y. Tu. Precision sensing by two opposing gradient sensors: how does
Escherichiacoli find its preferred pH level?
Biophysical Journal , 105(1):276-285, 2013.[51] M. Rousset and G. Samaey. Individual-based models for bacterial chemotaxis in the diffusionasymptotics.
Math. Models Methods Appl. Sci. , 23(11):2005-2037, 2013.[52] A. Gosztolai and M. Barahona. Cellular memory enhances bacterial chemotactic navigationin rugged environments.
Communications Physics , 3(1):1-10, 2020.[53] R. G. Endres and N. S. Wingreen. Precise adaptation in bacterial chemotaxis through assis-tance neighborhoods.
Proc. Natl. Acad. Sci. , 103(35):13040-13044, 2006.[54] G. Lan, S. Schulmeister, V. Sourjik, and Y. Tu. Adapt locally and act globally: strategy tomaintain high chemoreceptor sensitivity in complex environments.
Molecular systems biology ,7(1):475, 2011.[55] H. C. Berg and P. M. Tedesco. Transient response to chemotactic stimuli in
Escherichia coli . Proc. Natl. Acad. Sci. , 72(8):3235-3239, 1975.
A Appendix
A.1 Proof of Lemma 1
The equation of our interest is n t = (cid:0) Dn x − χV ( x ) n (cid:1) x = Dn xx − χV ( x ) n x − χV (cid:48) ( x ) n with boundary conditions Dn x (0 , t ) = χV (0) n (0 , t ) , Dn x ( L, t ) = χV ( L ) n ( L, t ) . Assume n ( x, t ) = ϕ ( x ) ψ ( t ). Then, it is satisfied ψ (cid:48) ( t ) ψ ( t ) = Dϕ (cid:48)(cid:48) ( x ) − χV ( x ) ϕ (cid:48) ( x ) − χV (cid:48) ( x ) ϕ ( x ) ϕ ( x ) =: − λ. To show that the solution n ( x, t ) is bounded, we prove that if λ exists, it is non-negative.33onsider the following eigenvalue problem: L ϕ ( x ) := Dϕ (cid:48)(cid:48) ( x ) − χV ( x ) ϕ (cid:48) ( x ) − χV (cid:48) ( x ) ϕ ( x ) = − λϕ ( x ) (EP)satisfying Dϕ (cid:48) (0) = χV (0) ϕ (0) , Dϕ (cid:48) ( L ) = χV ( L ) ϕ ( L ) . (BC)Putting (EP) into the Sturm-Liouville operator, we have L p ( x ) := ddx (cid:16) p ( x ) ddx (cid:17) + q ( x ) = − λσ ( x ) , (SL)where p ( x ) = e − (cid:82) χV ( x ) D dx > , q ( x ) = − χV (cid:48) ( x ) D p ( x ) ≥ , σ ( x ) = 1 D p ( x ) > . By Sturm-Liouville’s Theory, the problem (SL)-(BC) is naturally posed on H , where H ([0 , L ]) = (cid:8) u ∈ H ([0 , L ]) : Du x (0 , t ) = χV (0) u (0 , t ) , Du x ( L, t ) = χV ( L ) u ( L, t ) (cid:9) , and L p ( x ) is self-adjoint in the inner product < u, v > := (cid:90) L u ( x ) v ( x ) dx. Moreover, the eigenvalues and the corresponding normalized eigenfunctions of (SL)-(BC) satisfythe following properties:(a) All the eigenvalues are real, simple, and satisfy λ < λ < λ < · · · and lim n →∞ λ n = ∞ . (b) Each eigenfunction ϕ n ( x ) has n simple zeros in the open interval (0 , L ) . (c) < ϕ n , ϕ m > = δ nm .(d) { ϕ n ( x ) } ∞ n =0 forms a complete orthonormal basis of L ([0 , L ]).(e) The smallest eigenvalue λ is non-negative and satisfies L p ( x ) ϕ ( x ) = − λ σ ( x ) ϕ ( x ) , ⇒ < L p ( x ) ϕ ( x ) , ϕ ( x ) > = (cid:90) L ddx (cid:16) p ( x ) ddx ϕ ( x ) (cid:17) ϕ ( x ) (cid:124) (cid:123)(cid:122) (cid:125) =: I + q ( x ) ϕ ( x ) dx (cid:124) (cid:123)(cid:122) (cid:125) =: I = − (cid:90) L λ σ ( x ) ϕ ( x ) dx. ( ∗ )For simplicity we replace ϕ ( x ) and ddx by u ( x ) and (cid:48) , respectively. Then, by integration byparts, we have I = p ( x ) u ( x ) u (cid:48) ( x ) (cid:12)(cid:12)(cid:12) L − (cid:90) L p ( x )( u (cid:48) ( x )) dx, I = − χD (cid:90) L V (cid:48) ( x ) e − (cid:82) χV ( x ) D dx u ( x )= − χD V ( x ) e − (cid:82) χV ( x ) D dx u ( x ) dx (cid:12)(cid:12)(cid:12) L + χD (cid:90) L V ( x ) (cid:16) − χV ( x ) D (cid:17) e − (cid:82) χV ( x ) D dx u ( x ) dx + χD (cid:90) L V ( x ) e − (cid:82) χV ( x ) D dx u ( x ) u (cid:48) ( x ) dx = − χD V ( x ) p ( x ) u ( x ) (cid:12)(cid:12)(cid:12) L − (cid:90) L χ V ( x ) D p ( x ) u ( x ) dx + (cid:90) L χV ( x ) D p ( x )2 u ( x ) u (cid:48) ( x ) dx. Note that ( ∗ ) = I + I . Hence,( ∗ ) = p ( L ) u ( L ) (cid:16) u (cid:48) ( L ) − χD V ( L ) u ( L ) (cid:17) − p (0) u (0) (cid:16) u (cid:48) (0) − χD V (0) u (0) (cid:17) + (cid:90) L p ( x ) (cid:16) − u (cid:48) ( x ) − χ V ( x ) D u ( x ) + 2 χV ( x ) D u ( x ) u (cid:48) ( x ) (cid:17) dx, where the first two terms on the right hand side disappear due to (BC), and the integrandof the integral is non-positive since p ( x ) > u (cid:48) − χVD uu (cid:48) + χ V D u = (cid:16) u (cid:48) − χVD u (cid:17) ≥ . Therefore, from ( ∗ ), we arrive at λ = − < L p ( x ) ϕ , ϕ > (cid:82) L σ ( x ) ϕ ( x ) dx ≥ . . P a r a m e t e r s i n E . c o l i i n t e r n a l d y n a m i c s E q u a t i o n P a r a m e t e r D e s c r i p t i o n V a l u e R e f e r e n c e s M W C m o d e l ( ) N N u m b e r o f r ece p t o r s i n a c l u s t e r , c o m p o s e d o f T a r a nd T s r [ , ] r F r a c t i o n o f r ece p t o r T a r t o M e A s p a r F r a c t i o n o f r ece p t o r T s r t o s e r i n e a α F r eee n e r g y p e r a dd e d m e t h y l a t i o n g r o up . [ , , , ] m R e f e r e n ce m e t h y l a t i o n l e v e li n t h e f r eee n e r g y [ , , , ] K A D i ss o c i a t i o n c o n s t a n t o f M e A s p t o t h e a c t i v e r ece p t o r T a r . µ M [ , , , ] K A D i ss o c i a t i o n c o n s t a n t o f s e r i n e t o t h e a c t i v e r ece p t o r T s r m M [ , , ] K I D i ss o c i a t i o n c o n s t a n t o f M e A s p t o t h e i n a c t i v e r ece p t o r T a r µ M [ , , ] K I D i ss o c i a t i o n c o n s t a n t o f s e r i n e t o t h e a c t i v e r ece p t o r T s r µ M [ , ] A d a p t a t i o n m o d e l ( ) a A d a p t a t i o n l e v e l . [ , ] τ a A d a p t a t i o n t i m e v a r i e s b [ ] R un a nd T u m b l e m o t i o n ( ) λ R o t a t i o n a l d i ff u s i o n . r a d s − [ , ] H H ill c o e ffi c i e n t o f m o t o r ’ s r e s p o n s ec u r v e [ , ] τ R un a v e r ag e t i m e . s [ , ] ν R un v e l o c i t y . µ m s − [ , ] T r a n s f o r m e d i n t e r n a l d y n a m i c s ( ) q a . [ ] p α τ − a v a r i e s b γ R a t i o b e t w ee n T a r a nd T s rr ece p t o r s , r r − v a r i e s a T r a n s f o r m e d t u m b li n g r a t e ( ) r ( τ a H ) − [ ] a I n t h i s w o r k , w e a r e i n t e r e s t e d i n t h e r a t i oo f r a nd r s a t i s f y i n g r + r = . I n s t e a d o f t h e r a n g e o f r a nd r , w e p r e s e n tt h e r a n g e o f γ = r / r . b B a c t e r i a l a d a p t a t i o n t i m e v a r i e s , a nd i t d e p e nd s o n t h e s t r e n g t h o f s i g n a l s [ ]. I n t h i s w o r k , w e v a r y p b y c h oo s i n g d i ff e r e n t τ a b e t w ee n . nd s i n [ , , ]. T a b l e : P a r a m e t e r s u s e d i n i n t r a ce ll u l a r s i g n a li n g p a t h w a y o f E . c o l i . A n o v e r v i e w o f nu m e r i c a l s i m u l a t i o n s A b r i e f d e s c r i p t i o n o f M o n t e - C a r l o s i m u l a t i o n : I n ao n e - d i m e n s i o n a l ( r e s p ec t i v e l y , t w o - d i m e n s i o n a l ) c h a nn e l, w e l o c a t e a n e n s e m b l e o f , e n t s i n t h ece n t e r o f t h ec h a nn e l x = ( r e s p ec t i v e l y , ( x , y ) = ( , )) a tt i m e t = . A t e a c h t i m e s t e p , t h e i nd i v i du a l s c h oo s e a d i r ec t i o n + r - ( r e s p ec t i v e l y , ( c o s ( θ ) , s i n ( θ )) , θ ∈ [ , π )) a t r a nd o m , a nd m o v e i n t h a t d i r ec t i o n w i t h a c o n s t a n t s p ee d ν > . A t e a c h t i m e s t e p , t h e i n t e r n a l d y n a m i c s o f e a c h i nd i v i du a l a r ec o m pu t e db y E u l e r m e t h o d . A tt h ee nd o f e a c h t i m e s t e p , w ec h oo s e a nu m b e r b e t w ee n nd r a nd o m l y a nd c o m p a r e t h e nu m b e r w i t h t h e p r o b a b ili t y o f c h a n g e f r o m r un t o t u m b l e i n i n t e r v a l o f l e n g t h d t , n a m e l y λ ( a ) d t . I f t h e t u r n o cc u r s , t h ece ll m o v e s i n t h e o pp o s i t e d i r ec t i o n w i t h a p r o b a b ili t y o f . ( r e s p ec t i v e l y , r o t a t e s b y θ ∈ [ , π ) , w h e r e θ i s c h o s e n a t r a nd o m ) . I f a ce lli s l o c a t e d o u t s i d e t h e s p a t i a l d o m a i n , w e r e l o c a t e t h ece ll b y a pp l y i n g t h e b o und a r y c o nd i t i o n . E x p r e ss i o n V a l u e M o n t e - C a r l o s i m u l a t i o n N u m e r i c a l p a r t i a l d i ff e r e n t i a l e q u a t i o n s a Sp a t i a l d o m a i n ( x ) D : ≤ x ≤ ( µ m ) D : ≤ x ≤ ( µ m ) a nd ≤ y ≤ ( µ m ) T i m e d o m a i n ( t ) ≤ t ≤ ( s ec ) b I n i t i a l d a t a A n e n s e m b l e o f , e n t s p o s e s i n t h ece n t e r o f e a c hd o m a i n . D : e x p { − ( x − ) / ( ε ) } / √ π ε , ε = − D : e x p { − (( x − ) + ( y − ) ) / ( ε ) } / π ε , ε = − . B o und a r y c o nd i t i o n s N o flu x b o und a r y c o nd i t i o n s N o flu x b o und a r y c o nd i t i o n s : ( ) f o r D a nd ( ) f o r D Sp a t i a l s t e p s i ze ( ∆ x , ∆ y ) D : ∆ x = . ( µ m ) D : ∆ x = ∆ y = . ( µ m ) D : ∆ x = . ( µ m ) c D : ∆ x = ∆ y = . ( µ m ) c T i m e s t e p s i ze ( ∆ t ) ∆ t = . ( s ec ) d ∆ t = . ( s ec ) f o r D a nd ∆ t = . ( s ec ) f o r D c L i ga nd f un c t i o n ( S , S ) · D u a lli n e a r g r a d i e n t s i nS ec t i o n s . nd . : S ( x ) = . x + nd S ( x ) = − . x + · D u a l e x p o n e n t i a l g r a d i e n t s i nS ec t i o n s . nd . : S ( x ) = e . x a nd S ( x ) = e − . ( x − ) · M i x e d o pp o s i n gg r a d i e n t s i nS ec t i o n . : S ( x , y ) = ( . x + ) e . ( y − ) a nd S ( x , y ) = ( − . x + ) e − . ( y − ) T a r / T s rr a t i o ( γ ) A d a p t a t i o n e ( p ) · S ec t i o n s . nd . : γ = . nd γ = . w i t h p = . ( D ) , ( D ) . · S ec t i o n s . nd . : γ = . nd γ = . w i t h p = . ( D ) , . ( D ) · S ec t i o n . : γ = . nd p = a F i n i t e d i ff e r e n ce m e t h o d i s u s e d . b T h e s o l u t i o n s o f t h e M o n t e - C a r l o s i m u l a t i o n a nd a d v ec t i o n - d i ff u s i o n e q u a t i o n ss h o w n i nS ec t i o n s . t o6 . b ec o m e s t a t i o n a r y a t t = . c T h e t i m e s t e p s i ze a nd t h e s p a ce s t e p s i ze i n t h e fin i t e d i ff e r e n ce f o r m u l a s a t i s f y C o u r a n t - F r i e d r i c h s - L e w y ( C F L ) c o nd i t i o n a nd v o n N e u m a nn s t a b ili t y a n a l y s i s , r e s p ec t i v e l y . I t i s c o nfi r m e d t h a t u s i n g s m a ll e r s t e p s i ze s d o e s n o t a ff ec t o u rr e s u l t s a s l o n ga s t h e s t a b ili t y c o nd i t i o n s a r e s a t i s fi e d . d W ec h oo s e s m a ll v a l u e f o r ∆ t t o s o l v e t h e d y n a m i c s o f m e t h y l a t i o nb y E u l e r m e t h o d . e T h e v a l u e s o f p a nd γ s a t i s f y t h e s h a ll o w c o nd i t i o n . T a b l e : I npu t d a t a u s e d i nnu m e r i c a l s i m u l a t i o nn