Monte Carlo method for kinetic chemotaxis model and its applications - - PowerPoint PPT Presentation

Monte Carlo method for kinetic chemotaxis model and its applications on traveling pulse and pattern formation

Graduate School of Simulation Studies, University of Hyogo Shugo YASUDA

Stochastic Dynamics Out of Equilibrium @ IHP 11/04/2017

In collaboration with Benoît Perthame Laboratoire Jacques-Louis Lions

Plan of Talk

• 1. General Introduction
• 2. Kinetic Chemotaxis Model
• 3. Monte Carlo Method
• 4. Application 1: Traveling pulse
• 5. Application 2: Pattern formation
• 6. Concluding remarks

2

Collective dynamics of bacteria

Pattern formation by Budrene and Berg, Nature 349, 630 (1991).

2~4 µm

Homepage of H. C. Berg

http://www.rowland.harvard.edu/labs/bacteria/

Run-and-Tumble Bacteria

Introduction

• E. Coli

“Run”: Flagella rotate counterclockwise “Tumble”: flagella rotate clockwise Bacteria communicate via chemical cues

Motivation

• Multiscale mechanism and mathematical hierarchy

in the collective dynamics of bacteria.

– Relation between macroscopic phenomena, individual motions, and internal states

• Simulation method

– Extensible (modeling) and Scalable (computation)

• Applications

– Traveling pulse, Pattern formations, ....

4

Objective of study

• Development of a Monte Carlo method for

chemotactic bacteria based on a kinetic chemotaxis model.

• Applications on traveling pulse and pattern

formation.

– Validity of the MC method via comparisons to the theoretical and experimental results. – A new theoretical result on the instability analysis of a kinetic chemotaxis equation.

5

Why kinetic model?

• Mesoscopic modeling involving the individual

dynamics (multiscale nature)

– H. G. Othmer, S. R. Dunbar, and W. Alt (1988); R. Erban and H. G. Othmer (2004); Y. Dolak and C. Schmeiser (2005); N. Bellomo, A. Bellouquid, J. Nieto, and J. Soler (2007); etc..

• Mathematical hierarchy

– T. Hillen and H. G. Othmer (2000), (2002); F. A.C.C. Chalub, P. Markowich, B. Perthame, and C. Schmeiser (2004); F. James and N. Vauchelet (2013); G. Si, M. Tang, and X. Yang (2014); B. Perthame, M. Tang, and N. Vauchelet (2016).

• Development of experimental technologies

– J. Saragosti, V. Calvez, N. Bournaveas, B. Perthame, A. Buguin, and P. Silberzan (2011); C. Emako, C. Gayrard, A. Buguin, L. Almeida, and N. Vauchelet (2016).

6

Plan of Talk

• 1. Introduction
• 2. Kinetic Chemotaxis Model
• 3. Monte Carlo Simulation
• 4. Application 1: Traveling pulse
• 5. Application 2: Pattern formation
• 6. Concluding remarks

7

Schematic of kinetic modeling

8

2~4 µm bacterium Chemical cues << 1 µm Foods Secretions

Biased random motions searching for the chemical attractants

Kinetic description for bacterial density with the velocity distribution function 𝑔(𝑢, 𝒚, 𝒘) Continuum description for chemical cues, 𝑇(𝑢, 𝒚)

Schematic of kinetic modeling

9

Stiff response

From Block SM, Segall JE, Berg HC,

• J. Bacteriol 154, 312 (1983)

ramp up ramp down

Individual Motions of Bacteria

10

Run-and-Tumble motion

e.g., E-coli

Stochastic process

• 1. Tumbling at some rate 𝜇 .
• 2. Reorientation followed by

some PDF 𝐿(𝑤, 𝑤-).

• 3. Cell division/extinction with

some rate 𝑠. Bacterial density 𝑔(𝑢, 𝑦, 𝑤) changes during the stochastic process.

Homepage of H. C. Berg

http://www.rowland.harvard.edu/labs/bacteria/

𝜖1𝑔 + 𝑤 3 𝜖4𝑔 = 6 𝑈 𝑤, 𝑤- 𝑔 𝑢, 𝑦, 𝑤-

• 𝑒𝑤- − 6 𝑈 𝑤-, 𝑤 𝑔 𝑢, 𝑦, 𝑤
• 𝑒𝑤- + 𝑠𝑔(𝑢, 𝑦, 𝑤)

Kinetic Chemotaxis model with growth term

11

Gain Term Lost Term

Cell division

𝑈 𝑤, 𝑤- = 𝜇 𝑤- 𝐿(𝑤, 𝑤-) 6 𝐿 𝑤, 𝑤- 𝑒𝑤 = 1

• 𝑤′

Searching for foods and chemical cues along their trajectory

Transient kernel

Scattering Kernel

• Tumbling rate

– Stiff response function

12

𝜇 𝑤- = 1 2 𝜔? 𝐸 log 𝑂 𝐸𝑢 E

F-

+ 𝜔G 𝐸 log 𝑇 𝐸𝑢 E

F-

𝜔 𝑌 = 𝜔I − 𝜓 tanh 𝑌 𝜀

• Mean tumbling rate 𝜔I
• Modulation parameter 𝜓G,?
• Stiffness parameter 𝜀PQ

Temporal variation along the trajectory

Nutrient-Poor Nutrient-Rich

Scattering Kernel

• Reorientation (e.g., von Mises distribution)

– Reorientation angle q – Constant Speed 𝒘 = 𝑊

I.

– Standard deviation s

• Uniform scattering 𝐿 =

Q STU

V W as 𝜏 → ∞.

𝐿 𝒘, 𝒘- = exp − 1 − cos 𝜄 𝜏a 2𝜌𝑊

I a𝜏a 1 − 𝑓P a dW 𝜄 𝒘′ 𝒘

Basic equations

• Kinetic chemotaxis

– Modulation of tumbling frequency, for example, – PDF of reorientation angle, for example,

14

ˆ K(ˆ e, ˆ e0) = exp ⇣

1ˆ e·ˆ e0 σ2

⌘ 2πσ2 ⇣ 1 − e 2

σ2

⌘

|e0| = 1

(von Mises distribution)

Basic equations

• Reaction-Diffusion equations of chemical cues

15

ˆ ρ(ˆ t, ˆ x) = 1 4π Z

|ˆ e0|=1

ˆ f(ˆ t, ˆ x, ˆ e0)dΩ(ˆ e0)

Parameters

• Mean run length (or Knudsen number)
• Stiffness and modulation in response function
• Other Parameters

16

Plan of Talk

• 1. Introduction
• 2. Kinetic Chemotaxis Model
• 3. Monte Carlo method
• 4. Application 1: Traveling wave
• 5. Application 2: Pattern formation
• 6. Concluding remarks

17

Simulation method

• Monte Carlo method for chemotactic bacteria coupled

with a finite volume scheme for chemical cues.

• Motions of bacteria calculated by MC particles.
• Macroscopic quantities are calculated based on a

lattice-mesh system.

• Similar to the DSMC method for the Boltzmann

equation of gases.

18

Lattice System and MC Particles

• Motions of bacteria by Monte Carlo particles
• Macroscopic quantities on a lattice-mesh system

19

ˆ x0 ˆ x1 ˆ xn ˆ xn+1 ˆ xIx−1 ˆ xIx

Specular Specular Periodic non-flux ∂xS = ∂xN = 0 non-flux ∂xS = ∂xN = 0

x y z

superscript n: time step, subscript i: lattice site, and subscript (l): index of particle

i i + 1

ˆ rn

(l)

Calculation of Chemical Cues

• Finite volume scheme on the lattice mesh

20

Population densities are calculated from the numbers of MC particles in each lattice site.

Monte Carlo Method

• 0. Initialization: Distribute particles according to 𝑔

e

f0(𝒇

i).

• 1. Move particles in a time-step size ∆𝑢.
• 2. Calculation of local concentration of chemical cues.
• 3. Tumbling of each particle by a probability

𝜇 e(𝒇 i′ k )Δ𝑢̂ .

• 4. Reorientation angle by 𝐿 𝒇

i, 𝒇 i- .

• 5. Division by a probalibity 𝑠̂∆𝑢̂ .
• 6. Return to 1.

21

Monte Carlo Method

• 0. Initialization:

MC particles are distributed according to 𝑔 e

f I(𝒇

i).

– Calculate particle number in the i th lattice site 𝜈f

I via

– where 𝑥I is the uniform weight of a single MC particle – In each lattice site, particles are randomly distributed. – Velocity of particle is determined by the PDF 𝑔 e

f I(𝑓̂)/𝜍

i𝑗

22

Monte Carlo Method

• 1. Movement: Particles move with their velocities in a

time step size ∆𝑢.

– Particles beyond the boundaries are removed and new

• nes are inserted following the boundary conditions.

– Count the numbers of simulation particles in each lattice site 𝜈ftuQ (𝑗 = 0, … , 𝐽4 − 1).

23

𝒔 i(k), 𝒇 i(k): position and velocity

• f l th particle

Monte Carlo Method

• 2. Calculation of chemical cues at each lattice

site.

– with 𝜍 if

tuQ = 𝑥I𝜈ftuQ/[ 𝑀0∆𝑦 3𝜍I].

24

Monte Carlo Method

• 3. Tumbling of the l th particle by a probability

𝜔 }0∆𝑢̂ Ψ

• (𝒇

i k

t ),

25

Temporal variation along the pathway 𝒔 i(k)

t → 𝒔

i(k)

tuQ.

𝐸 log 𝑂 𝐸𝑢 E

𝒇

⟹ log𝑂(k)

tuQ − log𝑂(k) t

Δ𝑢

Monte Carlo Method

• 3. Tumbling of the lth particle by a probability

𝜔 }0∆𝑢̂ Ψ

• (𝒇

i k

t ),

– 𝑇 e(k)

t , 𝑂

• (k)

t : sensed by the lth MC particle at 𝒔

i(k)

t

– Calculated by the interpolation, – The particles that stay at the same lattice site after ∆𝑢 passes can sense the gradients.

26

Monte Carlo Method

• 4. Reorientation angle by a probability

𝐿

• (𝒇

i k

tuQ, 𝒇

i(k)

t ).

– Reorientation angle 𝜄 (for von Mises distribution)

• 5. Cell divisions (or deaths) with a probability 𝑠̂∆𝑢̂.
• 6. Return to step 1 (Movement Step).

27

𝜄 𝒇 i(k)

t

𝒇 i(k)

tuQ

Monte Carlo Method

28

• First order accuracy in time and space under the

assumption of the law of large numbers.

–

• B. Perthame and S. Yasuda, “Self-organized pattern formation of run-and-tumble chemotactic

bacteria: Instability analysis of a kinetic chemotaxis model”, hal-01494963 (2017).

Plan of Talk

• 1. Introduction
• 2. Kinetic Chemotaxis Model
• 3. Monte Carlo Simulation
• 4. Application 1: Traveling pulse
• 5. Application 2: Pattern formation
• 6. Concluding remarks

29

• Literature on traveling pulse

Vwave=4.1 µm/s by J. Saragosti, V. Calvez, N. Bournaveas, B. Perthame, A. Buguinn, and P. Silberzan, PNAS 108, 16235(2011)

Problem and parameter setting

31

• Initial condition and geometry
• Parameter setting
• Mean tumbling frequency 𝜔I = 3.0 [1/s] (𝜔

}I

PQ = 0.00833)

• Modulation of the response 𝜓̂G = 0.2 and 𝜓̂? = 0.6.
• Stiffness of the response δ = 0.125 [1/s].
• Division rate 𝑠̂ = 0.0067.
• Chanel length 𝑀

} = 18.

• Total particle number 56640.
• ∆𝑦

i = 0.025, ∆𝑢̂ = 0.005. (∆𝑢̂ < 𝜔 }I

PQ)

𝑦 i=0 𝑦 i=𝑀 } Initially accumulated at 𝑦 i=0 Specular reflection for bacteria at 𝑦 i = 0, 𝑀 }. Non-flux of chemical cues at 𝑦 i = 0, 𝑀 }.

𝜖4𝑇 e = 𝜖4𝑂

• = 0

𝜖4𝑇 e = 𝜖4𝑂

• = 0

𝑇 e = 0 and 𝑂

• = 1.

Other sides are periodic.

• J. Saragosti, V. Calvez, N. Bournaveas, B. Perthame,
• A. Buguin, and P. Silberzan, PNAS 108, 16235(2011)

Movie on the bacterial motions

32

x y Boundary condition Specular for x and periodic for y and z.

𝑓̂4

Time progress of population density

33

• Fig1. Time progress of population density of bacteria along the channel. (a) the snap

shots and (b) super position of the density profiles in the moving frame 𝑦 i∗ with a constant wave velocity Vwave=4.0 µm/s. (In experiment Vwave=4.1 µm/s.)

5 10 15 10 20 30 40

ˆ ρ ˆ x ˆ t = 25 ˆ t = 50 ˆ t = 75 ˆ t = 100

• 2
• 1

1 10 20 30 40

ˆ ρ ˆ x∗ ˆ Vwave = 0.16

(a) (b)

Application to the traveling wave

Comparison to the asymptotic analysis

34

• Diffusion scaling, a new reference time 𝑢′I
• Small modulation and small division rate

35

I[|r log F|] = Z 1 ζ tanh(δ−1|r log F|ζ)dζ

Chemotaxis Random walk Proliferation

Diffusion limit

Keller-Segel type

Comparison of Kinetic and Continuum

• Non-proliferation

– 𝑠 = 0

• Tumbling frequency

– 𝜁 = 0.02, 0.013, 0.01, 0.005, and 0.001

• Other Parameter

– 𝜚? = 72, 𝜚G = 24, 𝜀PQ = 0.2, – 𝑏 = 24, 𝑑 = 120, 𝐸Ž = 𝐸? = 3.84

36

MC vs. Continuum

• Snapshot of Population density

37

7 8 9 10 5 10 15 20 25

ˆ ρ ˆ x = 0.02 0.013 0.01 0.005 0.001 → 0

• Fig. 1 Comparison of the snapshots of population density of bacteria at t=0.5 between

various Knudsen numbers. No proliferation r=0

MC vs. Continuum

• Traveling speed
• Fig. 2 The convergence of the traveling speed in the continuum limit. The right-arrow

shows the result of the analytic formula obtained for the sign response function in the continuum limit (PLoS Comput. Biol. 6, e1000890 (2010) ).

10

2

10

3

14 15 16 17 18 19

ˆ Vwave/ ˆ δ−1=0.2 ˆ δ−1=1.0 ˆ δ−1 → ∞ → 0 1/

38

Effect of the stiffness and modulation

• Population density profile
• 1

1 25 50 75 100 125

ˆ ρ ˆ x∗

ˆ δ−1=0.15 0.2 0.25 0.5 1.0

(a)

• 1

1 10 20 30 40

ˆ ρ ˆ x∗

ˆ χN=0.5 0.6 0.7

(b)

Variation in stiffness 𝜀 ePQ Variation in modulation 𝜓̂?

• Fig. 4 The effect of the variations in stiffness and modulation parameters on the

population density profile in the moving frame 𝑦 i∗ .

Nutrient-Poor Nutrient-Rich

response function 39

Effect of the stiffness and modulation

• Traveling speed

0.5 1 0.14 0.15 0.16

ˆ Vwave ˆ δ−1 (a)

0.2 0.4 0.6 0.8 1 0.05 0.1 0.15 0.2 0.25

ˆ χN (b)

Present Analytic

Variation in stiffness 𝜀 ePQ Variation in modulation 𝜓̂?

• Fig. 5 The effect of the variations in stiffness and modulation parameters on the

traveling speed. The analytic formula is obtained for the sign response function in the continuum limit (PLoS Comput. Biol. 6, e1000890 (2010) ).

(𝜀 ePQ → ∞, 𝜔 }I

PQ → 0 )

40

Effect of the stiffness and modulation

• Velocity distribution at the peak of the wave
• 1
• 0.5

0.5 1 0.3 0.4 0.5 0.6 0.7 0.8

p(ˆ ex) ˆ ex

ˆ δ−1=0.125 0.15 0.2 0.25 0.5 1.0

(a)

ˆ δ−1increases

• 1
• 0.5

0.5 1 0.3 0.4 0.5 0.6 0.7 0.8

p(ˆ ex) ˆ ex

ˆ χN=0.5 0.55 0.6 0.65 0.7 0.8

(b)

ˆ χN increases

Variation in stiffness 𝜀 ePQ Variation in modulation 𝜓̂?

• Fig. 7 The effect of the variations in stiffness and modulation parameters on the

PDF of the velocity at the peak of the wave 𝑦 i∗ = 0.

41

Remarks on application 1

• Reproduce the experimental result.
• Recover the Keller-Segel equation in the continuum limit.
• Importance of the kinetic model for a small (but finite)

value of e.

• An orthogonal effect of the stiffness d and modulation c
• n the profile of population density and traveling speed.

42

Plan of Talk

• 1. Introduction
• 2. Kinetic Chemotaxis Model
• 3. Monte Carlo Simulation
• 4. Application 1: Traveling pulse
• 5. Application 2: Pattern formation
• 6. Concluding remarks

43

Basic equation

• Kinetic chemotaxis model with a population growth term
• Only one chemical attractant
• Biased Tumbling, Uniform scattering

– 𝐿 𝑌 = 1 − 𝐺[𝑌], – 𝐺 0 = 0, 𝐺- 0 > 0.

44

∂tf(t, x, v) + v · rf =1 k ⇢ 1 4π Z

V

K[Dt log S|v0]f(v0)dΩ(v0) K[Dt log S|v]f(v)

• + P[ρ]f(v)

−d∆S(t, x) + S(t, x) = ρ(t, x)

t ≥ 0, x ∈ R, v ∈ V ⊂ R : |v| = 1

Basic equation

• Growth term 𝑄[𝜍] : Saturated at 𝜍 = 1
• Stationary uniform solution

45

( Division at the rate 𝑄[𝜍]) (Extinction at the rate |𝑄 𝜍 |)

f(t, x, v) = S(t, x) = ρ(t, x) = 1

P[0] = 1, P[ρ] > 0, for 0 < ρ < 1, P[ρ] < 0, for ρ > 1, P[ρ] ' 1 ρ, for ρ ' 1.

Linear instability Condition

• The uniform solution is linearly unstable if the stiffness
• f the response function 𝐺-[0] is sufficiently large as
• In addition, the unstable eigenmodes are bounded and no

high frequency oscillations exist.

46

F 0[0] k > inf

λ

" 1 + k

kλ arctan(kλ) − 1

# (1 + dλ2)

• B. Perthame & S. Yasuda, “self-organized pattern formation of run-and-tumble chemotactic

bacteria: Instability analysis of a kinetic chemotaxis model”, hal-01494963 (2017).

Linear instability analysis

• Perturbation around the uniform state,
• Fourier transform on 𝒚 and Moment on 𝒘,

47

f(t, x, v) = 1 + g(x, v)eµt, S(t, x) = 1 + Sg(x)eµt, ρ(t, x) = 1 + ρg(x)eµt, λ: wave vector

ˆ g(λ, v) = 1 − k + i F 0[0]

1+dλ2 λ · v

1 + kµ + iλ · v ˆ ρg(λ)

ˆ ρ(λ) = 1 2 Z 1

−1

⇣ 1 − k + i F 0[0]λv

1+dλ2

⌘ (1 + kµ1 − ikλ(µ2 + v))) (1 + kµ1)2 + k2λ2(µ2 + v)2 dvˆ ρg(λ) µ1 = Re(µ), µ2 = Im(µ)/λ

• For non-trivial solution 𝜍

i“;

• No solutions at 𝜇 → ∞ for the first equation.
• 𝜈a = 0 always satisfies the second equation.

Linear instability analysis

48

α = 1 − k kλ , β = F 0[0] k(1 + dλ2), ξ = kλ 1 + kµ1

✓ α − β ξ ◆ [arctan(ξ(µ2 + 1)) − arctan(ξ(µ2 − 1))] − µ2β log ✓ 1 + 4µ2 ξ−2 + (µ2 − 1)2 ◆ = 2 − 2β µ2β [arctan(ξ(µ2 + 1)) − arctan(ξ(µ2 − 1))] + 1 2 ✓ α − β ξ ◆ log ✓ 1 + 4µ2 ξ−2 + (µ2 − 1)2 ◆ = 0

Linear instability analysis

49

• No solutions at 𝜇 → ∞ for the first equation.

α = 1 − k kλ , β = F 0[0] k(1 + dλ2), ξ = kλ 1 + kµ1

✓ α − β ξ ◆ [arctan(ξ(µ2 + 1)) − arctan(ξ(µ2 − 1))] + µ2β log ✓ξ−2 + (µ2 − 1)2 ξ−2 + (µ2 + 1)2 ◆ = 2 − 2β

l The RHS converges to 2 and the second term of LHS is always non-positive. l The first term of LHS converges to zero.

• When 𝜊 converges to a finite value or diverges, this is obvious because 𝛽, 𝛾 → 0.
• When 𝜊 → 0,

l No eigenmodes exist in the large-oscillation limit. | arctan(ξ(µ2 + 1)) − arctan(ξ(µ2 − 1))| =

• arctan

✓ 2ξ 1 + ξ2(µ2

2 − 1)

◆

• <
• arctan

✓ 2ξ 1 − ξ2 ◆

• = |2ξ + O(ξ2)|

Linear instability analysis

• Under an assumption 𝜈a = 0.

– 𝜈Q =

— ˜ − Q ™ > 0 ↔ 0 < 𝜊 < 𝑙𝑚.

• Instability condition

50

F 0[0] k > " 1 + k

kλ arctan(kλ) − 1

# (1 + dλ2)

(αξ − β)arctan(ξ) ξ = 1 − β

Kinetic Instability Diagram

51

0.5 1 1.5 2 5 10 15 20

A B C D

Kinetic Instability Diagram. Chemotaxis-induced instability takes place when the parameter value of (

• ž I

™

,

Ÿ ™) exceeds the critical line for each value of 𝑙.

Instable Stable ?

Monte Carlo results

52

A B C D

Slightly above the critical line (yellow) 𝑙 = 1

Monte Carlo results

53

A B C D

Slightly below the critical line (green) 𝑙 = 2

Monte Carlo results

54 0.5 1 1.5 2 5 10 15 20

A B C D

Instable Stable

Time progress x √ k ρ

(a)

Time progress x √ k ρ

(b)

Time progress x √ k ρ

(c)

Time progress x √ k ρ

(d)

Time progress x √ k ρ

(e)

Time progress x √ k ρ

(f)

Stationary periodic No patterns

Monte Carlo results

55

10 10

1

10

2

10

• 5

10

• 4

10

• 3

10

• 2

10

• 1

10 10

1

10 10

1

10

2

10

• 5

10

• 4

10

• 3

10

• 2

10

• 1

10 10

1

|ˆ ρ(l)|2/k l √ k (a)

A for k=1.0 B for k=0.1 C for k=1.0 Continuum

|ˆ ρ(l)|2/k l √ k (b)

A for k=2.0 B for k=1.0 C for k=2.0 D for k=1.0

Power spectra of population density

Periodic patterns No patterns

• The unstable frequencies remain bounded as in the Turing instability.
• Neither growth nor damping at high oscillations in the kinetic results.

Concluding remarks

• A Monte Carlo method for run-and-tumble

chemotactic bacteria.

• Chemotaxis-induced instability condition in a

kinetic chemotaxis equation with growth term.

• The validity of the MC method is strengthened

via the comparison with the experimental (from a literature) and theoretical results.

• Future works

– Applications; Traveling waves, 2D pattern – Development; Internal states (or Memories)

56

Thank you very much

Acknowledgements Financial supports: IHP RIP program, CNRS, and JSPS

References [1] S. Yasuda, “Monte Carlo simulation for kinetic chemotaxis model: An application to the traveling population wave”, J. Comput. Phys. 330, 1022–1042 (2017). [2] B. Perthame and S. Yasuda, “Self-organized pattern formation of run-and-tumble chemotactic bacteria: Instability analysis of a kinetic chemotaxis model”, hal-01494963 (2017).

57