검색
Article Search

JMB Journal of Microbiolog and Biotechnology

OPEN ACCESS eISSN 0454-8124
pISSN 1225-6951
QR Code

Article

Kyungpook Mathematical Journal 2023; 63(2): 141-154

Published online June 30, 2023 https://doi.org/10.5666/KMJ.2023.63.2.141

Copyright © Kyungpook Mathematical Journal.

Numerical Nonlinear Stability of TravelingWaves for a Chemotaxis Model

Min-Gi Lee

Department of Mathematics, Kyungpook National University, Daegu 41566, Republic of Korea
e-mail : leem@knu.ac.kr

Received: March 26, 2023; Revised: April 8, 2023; Accepted: April 20, 2023

We study the stability of traveling waves of a certain chemotaxis model. The traveling wave solution is a central object of study in a chemotaxis model. Kim et al. [8] introduced a model on a population and nutrient densities based on a nonlinear diffusion law. They proved the existence of traveling waves for the one dimensional Cauchy problem. Existence theory for traveling waves is typically followed by stability analysis because any traveling waves that are not robust against a small perturbation would have little physical significance. We conduct a numerical nonlinear stability for a few relevant instances of traveling waves shown to exist in [8]. Results against absolute additive noises and relative additive noises are presented.

Keywords: traveling waves, stability, chemotaxis

In this paper, we study the stability of traveling waves of a certain chemotaxis model. In a chemotaxis model, formation of traveling fronts is a central theme of the study: let a long tube be filled with a nutrient that is to be consumed by a certain micro-organism, such as bacteria. Let a certain amount of bacteria population be placed at the leftmost entrance of the tube. Conversion of nutrient to a bacteria population at some rate would take place and eventually results in the exhastion of nutrient at the place. If the bacteria is capable of moving from one place to another, then those that move right will find a new opportunity of further growth with a high nutrient density at the new place. As a proof that this does take place, traveling wave solutions where the nutrient front retreats, and that of the bacteria population propagates, comes with an important physical relevance.

Kim et al. [8] introduced a following chemotaxis model:

ut=Δ(γ(n)u)+r0nu,nt=ϵΔnnufor (t,x)(0,T)×Ω,(γ(n)u)ν=nν=0for (t,x)(0,T)×Ω,u(0,x)=u0(x),n(0,x)=n0(x)for xΩ.

Here Ωn is a simply connected bounded open set with smooth boundary, u is the population density, n is the nutrient density, r0>0 is the growth rate, and ϵ>0 is a small diffusion constant for the nutrient. Initial and boundary data are assigned as above, where zero flux boundary conditions are imposed for both u and n. ν is the outward normal vector. For the nonlinearty γ:[0,), we assume

γissmoothandboundedineverycompactsetof[0,)γ(n)γ0>0foreverynforsomeγ0>0.

This chemotaxis model has been suggested from the consideration of non-fickian flux law for diffusion, which is an independent and substantial subject of the random walk in a heterogeneous medium. See the series of studies [7, 5, 6]. The key difference from the classical Keller-Segel model [4] is to have a flux function as in the first equation of (1.1)

J=p(n,u),p(n,u):=γ(n)u.

As a consequence, migration of population is from the place where p is high to the place p is low. The bacteria are supposed to sense this pressure p by sensing the values of u and n at the place. The value of p is decided by the relative values of n and u. For instance, if p is decreasing in n and increasing in u then migration occurs from the place where the population is crowded, relative to the nutrient, to the place that is less so. This type of chemotaxis model has attracted considerable attention in recent years for a variety of applications (see [7, 1, 14]). Existence theorems on related models also have been actively studied by those authors in [12, 13, 3, 2, 11, 9, 10].

The measurements of this pressure p is done by measurements of n and u at the place by the bacteria and the value of p is decided by the relative values of n and u.

Consider one space dimensional Cauchy problem

ut=(γ(n)u)xx+r0nu,nt=ϵnxxnufor (t,x)(0,)×.

Kim et al. [8] proved the existence of traveling wave solutions of (1.3) for the case ϵ=0 that are of the form

u(t,x)=u¯(xct),n(t,x)=n¯(xct),ξ=xct, (u¯,n¯)(u*,0)asξ,(u¯,n¯)(0,n*)asξ.

Plugging in the ansatz (1.4) to (1.3), one obtains the system of ordinary differential equations

cu¯=(γ(n¯)u¯)+r0n¯u¯,c n¯=ϵnn¯u¯.

The existence of traveling waves are proved by phase space analysis of (1.6).

Because those traveling waves are constructed by solving the system (1.6), not the system (1.3), once existence is established, the stability against a small perturbation as a solution of (1.3) is subjected to a further study. Any of traveling waves that is not robust against a small perturbation would have little physical significance.

In general, nonlinear stability is a difficult problem, in particular around this non-steady state solution. The objective of this paper is thus to study numerical nonlinear stability of traveling waves: we numerically construct the traveling wave solutions of (1.6) and the computed data are prepared as an initial data of the system (1.3). In this preparation, a small perturbation is added. The system (1.3) with the perturbed initial data is then numerically solved until a finite time T>0. The precise test procedure will be specified in Section 4.

The remainder of this paper is organised as follows. In Section 2 we summarize the existence results of traveling waves presented in [8] for completeness of our discussion. In Section 3 we present numerical integration results done by python SciPy Library for a few relevant instances of traveling waves. In Section 4 we present the numerical stability results, and finally we give conclusions in Section 5.

In this section, the system of odes (2.1) is further discussed, which is a summary of what are presented in [8]. For notational simplicity overlines on variables are suppressed.

Since nu=cn for the case ϵ=0, the first equation of (2.1) is integrated once to have

cu=(γ(n)u)+r0cn+P0,cn=nu,

for some integration constant P0. By imposing the condition (u,n)(0,n*) as ξ and using new variables (p,n)=(γ(n)u,n), the system can be arranged in the form

γ(n)p=cpr0c(nn*)γ(n),γ(n)n=npc.

Here, we have used the assumption (1.2) that γ(n)γ0>0 for every n.

2.1. Equilibrium points and linearization

It is straightforward to verify that (2.1) has only two equilibrium points that are

(p0,n0)=(r0n*γ(0),0),(p1,n1)=(0,n*).

This implies that u* in (1.5) is not an independent parameter but is decided as

u*=r0n*.

Linearization around each equilibrium point is summarized below.

1. (p0,n0) is a saddle point:

Eigenvalues: λ1=cγ(0),λ2=r0n*c,Eigenvectors: X1= 1 0 X2= x 1 ,where x=r0c2γ(0)r0n*γ(0)+c2(n*γ(0)γ(0)-1).

2. If c>c*=2r0n*γ(n*) then (p1,n1) is a stable node with two real and negative eigenvalues:

Eigenvalues: Tworootsλ˜1andλ˜2ofλ2+cγ(n*)λ+r0n*γ(n*)=0,Eigenvectors: X˜1= 1 n* cγ(n* ) λ ˜1 X˜2= 1 n* cγ(n* ) λ ˜2 .

Phase space analysis [8] based on the discussion so far gives the following result:

Theorem 2.1. Fix n*>0 and c>c*. There exists a heteroclinic orbit that is an intersection of the unstable manifold of (p0,n0) and the stable manifold of (p1,n1). The heteroclinic orbit entirely lies in the first quadrant of pn-plane.

As stated in Theorem 2.1, the traveling wave exists as a heterocilinic orbit of (2.1) that is a saddle-to-stable node connection. If (p^,n^) is any point on the stable manifold of (p0,n0), integrating numerically (2.1) from (p^,n^) forward in time, capturing the saddle-to-stable node connection, is not a stiff problem.

Based on that, in this section we numerically compute the saddle-to-stable node connection. We specify the procedure below.

1. We first let δ>0 sufficiently small and set

(p^,n^)=(p0,n0)+δX2,

where X2 is the eigenvector in Section 2.1.

2. Let M>0 be sufficiently large. We use SciPy Library to integrate (2.1) on interval [0,M] with initial data

(p(0),n(0))=(p^,n^).

We present three instances of traveling waves. As discussed earlier, γ that is decreasing in n corresponds to a chemotaxis model where migration of population uccurs from a crowded place to uncrowded place. Note that the value x in the eigenvector X2 will be negative in this case because γ(0)<0. γ that is increasing in n also is relevant to a certain application. In this case, the sign of x in the eigenvector X2 is decided by n*γ(0)γ(0)1. In case x>0, one observes an overshoot of p profile before it converges to 0 as ξ. See Figure 1 (c).

Figure 1. Traveling waves of cases A, B, and C.

Considering the above, those three instances are specified below, including various wave speeds and parameters.

[A]γ:n11+n1.2,r0=0.7,n*=1.1,c=1.5665c*(1+0.3),[B]γ:n1+n1.5,r0=1.8,n*=1.0,c=3.5474c*(1+0.1),[C]γ:n1+n,r0=1.2,n*=2.1,c=5.5900c*.

Numerically captured heteroclinic orbits are presented in Figure 1. In every case, n has a retreating front and p has a propagating front. In wave C, we see the overshoot of p profiles. On the other hands, in Figure 1 (a) and (b) where we have x<0, we observe p decreasing all the time in ξ converging to 0 and n increasing all the time converging to n*.

Any traveling wave solution that is not robust under a perturbation has little physical significance. In this section, we conduct the numerical nonlinear stability tests on three captured data in the preceding section. The specification of stability test procedure is presented below.

4.1. Procedure of numerical stability test

1. Points on the heteroclinic orbit are collected, for j=0,,N1 with h=MN, that are

(ξj,p(ξj),n(ξj)),ξj=jh.

Note that the heteroclinic orbits have been computed on sufficiently large interval [0,M] so that profiles of orbits are sufficiently flat near end points. See Figure 1.

2. Preparation of (u,n) profiles. Since variables we integrate are p and n in heteroclinic orbit, we take for j=0,,N1

(ξj,φ(ξj)):=(ξj,u(ξj),n(ξj))=(ξj,p(ξj)γ(n(ξj)),n(ξj)).

3. Preparation of initial data U0. The system we will run is the cauchy problem (1.3) in the wholse space . In its approximation, computational domain has to be large enough so that within a certain time T>0 traveling wave phenomena can be observed. For this reason, we take an interval [0,M¯] with M¯ possibly larger than M as a computational domain for (1.3). The heteroclinic orbit data is translated and extended to U0 data as follows.

Fori=0,,N¯1xi=ih[0,M¯],N¯=M¯h,(xi,U0(xi))= (ξ0 ,φ(ξ0 )) ifi<i0 (ξii0 ,φ(ξii0 )) ifi0 i<i0 +N (ξN1 ,φ(ξN1 )) ifii0 +N

for i0 an integer shift factor.

4. Preparation of a perturbed initial data U˜0. For a given U0, we consider two types of absolute and relative perturbations: For i=0,,N¯1

U˜0(xi)=U0(xi)+η¯E(xi)  for absolute additive noise, or U˜0(xi)=U0(xi)(1+ηE(xi))  for relative additive noise

for some E with E=1. η¯,η>0 are small parameters. η¯ measures the level of absolute noises and η measures the level of relative noises.

5. Finally, we numerically solve (1.3) in [0,T]×[0,M¯]. The time interval is discretized as

Form=0,,K¯1tm=km,k=TK¯.

For the parabolic solver, we simply used the explicit scheme:

u˜(tm,xj)=u˜(tm1,xj)+kr0n˜(tm1,xj)u˜(tm1,xj)  +kh2(p(tm1,xj1)2p(tm1,xj)+p(tm1,xj+1),  wherep(tm,xj):=u˜γ(n˜)(tm,xj),andn˜(tm,xj)=n˜(tm1,xj)kn˜(tm1,xj)u˜(tm1,xj)  +ϵkh2(n˜(tm1,xj1)2n˜(tm1,xj)+n˜(tm1,xj+1)

for m=1K¯1 and j=1,,N¯2. Assuming M¯ is large enough, we assigned the Dirichilet boundary condition, that is for every m

(u˜,n˜)(tm,x0)=φ(ξ0),(u˜,n˜)(tm,xN¯1)=φ(ξN1).

4.2. Results

In Figure 2 and 3 are the tests of wave A against an oscillatory sinusoidal perturbation, i.e.,

Figure 2. Test of wave A against the sinusoidal absolute additive noise. For the computation, h=0.1, k=0.005, η¯=0.01 have been used.
Figure 3. Test of wave A against the sinusoidal absolute additive noise. For the computation, h=0.1, k=0.005, η¯=0.05 have been used.
E(x)=sin2πxM¯l

with l=10.

We first consider the severer absolute noises. Figure 2 presents snapshots of p=uγ(n) and n profiles taken at time t=0, 0.5, 1.0, and 5.0 with η¯=0.01. The value of η¯ is the 1% of the maximum (p,n). Interesting behaviors against the absolute additive noise are exhibited: Oscillations first seem to be stabilized as time proceeds so that by the time around t=1.0 we see almost plateu for both p and n, but this attenuated oscillation is then amplified in later time. This behavior is more drastically observed with absolute noise at η¯=0.05 in Figure 3.

This illustrates that the shape of traveling wave is maintained for a short time in a stable manner, but due to the reaction of the system, traveling wave can be oscillatory in later time.

Note that the two end states at ± have nearly 0 state for one of n and u, or the added error in Figure 2 and 3 in relative sense are huge. Furthermore, in the second equation (1.3) for n, currently we have no diffusion (ϵ=0), which is an extreme case. Considering those, behaviors seems not drastically bad in the sense that the traveling wave does not disappear.

In Figure 4, 5, and 6 are tests against relative additive noises of 10% levels. We conducted computations until t=5.0. What are in the last row of each of Figure 4, 5, and 6 are the snapshots at t=1.0, 2.0, 3.0, 4.0, and 5.0 from the left to right. We considered an ambitious white noise for these tests where E(xi) is a random variable valued in [-1,1] with a uniform probability. To add this white noise implies that the initial data are discontinuous and highly oscillatory.

Figure 4. Snapshots of u and n.
Figure 5. Snapshots of u and n.
Figure 6. Snapshots of u and n.

We see in Figure 4 (a) the farily clean traveling wave fronts of u against the relative additive white noises. White noise in n cannot be flatten out simply because the equation for n in (1.3) is an ode.

In Figure 6 (a) for wave C, we presented the profile of p in place of u, to illustrates the overshoot before the convergence to 0 as x. Figure 6 (a) clearly shows the recovery to the clean profile of p traveling with an overshoot.

In this paper, numerical stability of traveling waves of (1.3) against perturbations has been tested on three instances. The traveling waves themselves were computed by Python SciPy Library, and are assigned as initial data for the system (1.3) after adding a perturbation. We tested both the absolute and relative perturbations. Considering that the end states of traveling waves contain the zero state in one of u and n, the added absolute additive perturbations near the zero state are huge ones if measured in relative sense. Yet we observed fairly stable behaviors against the absolute additive noise in Figure 2 and 3. Against the relative additive noise, traveling waves seem to be highly robust as seen in Figure 4, 5, and 6.

  1. M. Alfaro, T. Giletti, Y. -J. Kim, G. Peltier and H. Seo, On the modelling of spatially heterogeneous nonlocal diffusion: deciding factors and preferential position of individuals, J. Math. Biol., 84(2022), 1-35.
    Pubmed CrossRef
  2. L. Desvillettes, P. Laurencot, A. Trescases and M. Winkler, Weak solutions to triangular cross diffusion systems modeling chemotaxis with local sensing, Nonlinear Analysis, 226(2023), 113153.
    CrossRef
  3. K. Fujie and J. Jiang, Boundedness of classical solutions to a degenerate Keller-Segel type model with signal-dependent motilities, Acta Appl. Math., 176(3)(2021), 1-36.
    CrossRef
  4. E. F. Keller and L. A. Segel, Traveling bands of chemotactic bacteria: a theoretical analysis, J. Theor. Biol., 30(2)(1971), 235-248.
    Pubmed CrossRef
  5. H. -Y. Kim, Y. -J. Kim and H. -J. Lim, Heterogeneous discrete kinetic model and its diffusion limit, Kinet. Relat. Models, 14(5)(2021), 749-765.
    CrossRef
  6. J. Chung, Y-J Kim and M-G Lee, Random walk with heterogeneous sojourn time, arXiv preprint, (2023).
  7. Y-J Kim and C. Yoon, Modeling bacterial traveling wave patterns with exact crossdiffusion and population growth, Discrete Contin. Dyn. Syst. - B, (2023).
    CrossRef
  8. Y-J Kim, M. Mimura and C. Yoon, Nonlinear diffusion for bacterial traveling wave phenomenon, Bull. Math. Biol., 85(5)(2023).
    Pubmed CrossRef
  9. W. Lyu and Z.-A. Wang, Global classical solutions for a class of reaction-diffusion system with density-suppressed motility, Electron. res. arch., 30(3)(2022), 995-1015.
    CrossRef
  10. Y. Tao and M. Winkler, Global solutions to a Keller-Segel-consumption system involving singularly signal-dependent motilities in domains of arbitrary dimension, J. Differential Equations, 343(2023), 390-418.
    CrossRef
  11. Z. -A. Wang and J. Zheng, Global boundedness of the fully parabolic Keller-Segel system with signal-dependent motilities, Acta Appl. Math., 171(2021), 1-19.
    CrossRef
  12. M. Winkler, Global generalized solvability in a strongly degenerate taxis-type parabolic system modeling migration-consumption interaction, Z. Angew. Math. Phys., 74(1)(2023).
  13. M. Winkler, Can simultaneous density-determined enhancement of diffusion and cross-diffusion foster boundedness in Keller-Segel type systems involving signaldependent motilities?, Nonlinearity, 33(12)(2020), 6590-6623.
  14. C. Yoon and Y-J Kim, Bacterial chemotaxis without gradient-sensing, J. Math. Biol., 70(2015), 1359-1380.