A Positivity Preserving Centrul-Upwind Scheme for Chemotaxis and - PowerPoint PPT Presentation

A Positivity Preserving Centrul-Upwind Scheme for Chemotaxis and Haptotaxis Models Alina Chertock North Carolina State University chertock@math.ncsu.edu joint work with Alexander Kurganov

Outline • Chemotaxis model: The Keller-Segel Model • Numerical methods • Derivation of a new positivity-preserving numerical scheme • Numerical examples • Related models: – Chemotactic Bacteria Patterns in Semi-Solid Medium – Chemotactic Bacteria Patterns in Liquid Medium – Haptotaxis 1

Chemotaxis active orientation of cells and organisms along chemical gradients • Numerous examples in animal and insect ecology, biological and biomedical sciences: – animals and insects rely on an acute sense of smell for conveying information between members of the species; – bacterial infection: invades the body and may be attacked by movement of cells towards the source as a result of chemotaxis; – development of cancer: very much related to the ability of cancerous cells to move, and thus spread faster that healthy cells. • There are typically two kinds of patterns: – traveling waves (e.g., periodic swarm rings or band dynamics); – aggregate formation. 2

3

PDE Based Chemotaxis Models • typically two- or three-dimensional; • highly nonlinear; • described by time-dependent systems of PDEs, consisting of three distinct sets of terms: – reaction terms – model the interaction of different components (e.g., growth of cells, release of chemoattractant, etc.); – diffusion terms – model the random motion of each component; – chemotaxis terms – model the directed motion of one or more components in response to concentration gradient of another component (the chemoattractant). 4

Chemotaxis: the Keller-Segel Model [Patlak (1953) and Keller & Segel (1970,71)] � ρ t + ∇· ( χρ ∇ c ) = ∆ ρ x = ( x, y ) T ∈ Ω , t > 0 αc t = ∆ c − c + ρ α = 1 : parabolic case, α = 0 : parabolic-elliptic case Initial conditions: ρ ( x, y, t = 0) = ρ 0 ( x, y ) , c ( x, y, t = 0) = c 0 ( x, y ) Flux boundary conditions: ∇ ρ · n = ∇ c · n = 0 , x ∈ ∂ Ω , t > 0 • ρ ( x, y, t ) — cell density, • c ( x, y, t ) — chemoattractant concentration, • χ — chemotactic sensitivity constant. 5

Analytical Results [Herrero, Medina and Vel´ azquez (1997), Horstmann (2003, 2004), Perthame (2007), ...]: Recent review: [D. Horstmann; 2003, 2004] Recent book: [B. Perthame; 2007] • 1-D case – there are global smooth and unique solutions; • 2-D case – global existence depends on a threshold: – initial mass lies below the threshold → solutions exist globally; – initial mass lies above the threshold → solutions blow up in finite time; • Various regularizations. 6

Numerical Methods finite-difference [Tyson, Stern & LeVeque (2000)], finite element [Marocco (2003), Saito (2007)], finite volume [Filbet (2006)] methods, ... Difficulties: • highly nonlinear chemotaxis and reaction terms; • diffusion terms, which have infinite speed of propagation in the context of the solution; • patterns are sought on large domains across which a disturbance must propagate for some time without hitting the boundary; • local linear instability. Approaches: • to treat all terms simultaneously – typically need to be an implicit method due to diffusion and reaction terms; • to apply fractional step methods to treat each term separately; 7

Numerical Methods finite-difference [Tyson, Stern & LeVeque (2000)], finite element [Marocco (2003), Saito (2007)], finite volume [Filbet (2006)] methods, ... Difficulties: • highly nonlinear chemotaxis and reaction terms; • diffusion terms, which have infinite speed of propagation in the context of the solution; • patterns are sought on large domains across which a disturbance must propagate for some time without hitting the boundary; • local linear instability. Approaches: • to treat all terms simultaneously – typically need to be an implicit method due to diffusion and reaction terms; • to apply fractional step methods to treat each term separately; Goal: to develop a computationally efficient and stable method which can capture sharp gradients – difficult! 7

Keller-Segel Model – Numerical Example � ρ t + ∇· ( χρ ∇ c ) = ∆ ρ c t = ∆ c − c + ρ • Square domain Ω = [ − 1 2 , 1 2 ] × [ − 1 2 , 1 2 ] . • Initial conditions: ρ ( x, y, 0) = 1000 e − 100( x 2 + y 2 ) , c ( x, y, 0) = 500 e − 50( x 2 + y 2 ) . • Nuemann boundary cinditions. 8

Keller-Segel Model – Numerical Example � ρ t + ∇· ( χρ ∇ c ) = ∆ ρ c t = ∆ c − c + ρ • Square domain Ω = [ − 1 2 , 1 2 ] × [ − 1 2 , 1 2 ] . • Initial conditions: ρ ( x, y, 0) = 1000 e − 100( x 2 + y 2 ) , c ( x, y, 0) = 500 e − 50( x 2 + y 2 ) . • Nuemann boundary cinditions. According to theoretical results [Harrero, Vel´ azquez (1997)], both ρ - and c -components of the solution are expected to blow up at the origin in finite time. 8

Na¨ ıve Finite-Difference Scheme � ρ t + ( χρc x ) x + ( χρc y ) y = ρ xx + ρ yy c t = c xx + c yy − c + ρ  H y 2 − H y H x 2 ,k − H x   dρ j,k j + 1 j − 1 j,k + 1 j,k − 1  2 ,k  + D 2 2 = − − 0 ρ j,k ∆ x ∆ y dt   dc j,k  = D 2  0 c j,k − c j,k + ρ j,k dt where 2 ,k = χρ j +1 ,k + ρ j,k · c j +1 ,k − c j,k H x j + 1 2 ∆ x 2 = χρ j,k +1 + ρ j,k · c j,k +1 − c j,k H y j,k + 1 2 ∆ y 0 ρ j,k = ρ j +1 ,k − 2 ρ j,k + ρ j − 1 ,k + ρ j,k +1 − 2 ρ j,k + ρ j,k − 1 D 2 (∆ x ) 2 (∆ y ) 2 9

Small Times 10

Later Times t=4.4 ⋅ 10 −5 4 5x 10 4 3 2 1 0 −0.5 0 0.5 11

Later Times t=10 −4 6 x 10 8 6 4 2 0 −2 −0.5 0 0.5 12

Keller-Segel Model – Properties � ρ t + ( χρ c x ) x + ( χρ c y ) y = ρ xx + ρ yy c t = c xx + c yy − c + ρ Denote u := c x and v := c y and rewrite the Keller-Segel system   ρ t + ( χρ u ) x + ( χρ v ) y = ρ xx + ρ yy   u t − ρ x = u xx + u yy − u    v t − ρ y = v xx + v yy − v This is a system of convection-diffusion-reaction equations: U t + f ( U ) x + g ( U ) y = ∆ U + R ( U ) U := ( ρ, u, v ) T , f ( U ) := ( χρu, − ρ, 0) T , g ( U ) := ( χρv, 0 , − ρ ) T , R ( U ) := (0 , − u, − v ) T . 13

Keller-Segel Model – Properties U t + f ( U ) x + g ( U ) y = ∆ U + R ( U )           ∆ ρ 0 ρ χρu χρv         −   + − ρ + 0 = ∆ u u u v 0 − ρ ∆ v v t x y The Jacobians of f and g are:     χu χρ 0 χv 0 χρ ∂ f ∂ g   ,   − 1 0 0 0 0 0 ∂ U = ∂ U = 0 0 0 − 1 0 0 Their eigenvalues are: � � � u 2 − 4 ρ 1 , 2 = χ λ f λ f u ± 3 = 0 , 2 χ � � � 1 , 2 = χ v 2 − 4 ρ λ g λ g v ± , 3 = 0 2 χ 14

Keller-Segel Model – Properties � � � u 2 − 4 ρ 1 , 2 = χ λ f λ f u ± 3 = 0 , 2 χ � � � v 2 − 4 ρ 1 , 2 = χ λ g λ g v ± 3 = 0 , 2 χ The key (new) observation: the “purely” convective system U t + f ( U ) x + g ( U ) y = 0 is • hyperbolic (real e-values) if both χu 2 ≥ 4 ρ and χv 2 ≥ 4 ρ χ min( u 2 , v 2 ) < 4 ρ • elliptic (complex e-values) if Notice that the ellipticity condition is satisfied in generic cases, for example, when u = c x = 0 and ρ > 0 . The operator splitting approach may not be applicable! 15

Semi-Discrete Central-Upwind Scheme Central-upwind schemes were developed for multidimensional hyperbolic systems of conservation laws in 2000–2007 by Kurganov, Lin, Noelle, Petrova, Tadmor, ... Central-upwind schemes are Godunov-type finite-volume projection- evolution methods: • at each time level a solution is globally approximated by a piecewise polynomial function, • which is then evolved to the new time level using the integral form of the conservation law system. 16

Semi-Discrete Central-Upwind Scheme [Kurganov, Tadmor, Levy, Petrova, Noelle, Lin, Balbas, ...] U t + f ( U ) x + g ( U ) y = ∆ U + R ( U ) Divide the domain into cells: C j,k := [ x j − 1 2 ] × [ y k − 1 2 ] 2 , x j + 1 2 , y k + 1 � � 1 Denote: U j,k ( t ) := U ( x, y, t ) dxdy ∆ x ∆ y C j,k Evolve in time: H y 2 − H y H x 2 ,k − H x d j + 1 j − 1 j,k + 1 j,k − 1 2 ,k + Λ h 2 dt U j,k = − − j,k + R j,k ∆ x ∆ y j,k = U j +1 ,k − 2 U j,k + U j − 1 ,k + U j,k +1 − 2 U j,k + U j,k − 1 Λ h (∆ x ) 2 (∆ y ) 2 R j,k = (0 , − u j,k , − v j,k ) T . 17

Godunov-Type Schemes for Conservation Laws: projection-evolution methods       U E , W , N , S   ( t )       j,k     H x 2 ,k ( t ) j + 1 { U j ( t ) } → � a ± U ( · , t ) → 2 ,k ( t ) → → { U j,k ( t +∆ t ) } j + 1     H y     2 ( t )     j,k + 1  b ±  2 ( t ) j,k + 1 via either a fully-discrete scheme � � � � j,k − ∆ t − ∆ t n +1 n H x 2 ,k − H x H x 2 − H x +Λ h = U j,k + R j,k U j + 1 j − 1 j,k + 1 j,k − 1 j,k ∆ x 2 ,k ∆ y 2 or a semi-discrete scheme H y 2 − H y H x 2 ,k − H x d j + 1 j − 1 j,k + 1 j,k − 1 2 ,k + Λ h 2 dt U j,k = − − j,k + R j,k ∆ x ∆ y 18