Open access peer-reviewed chapter - ONLINE FIRST

Exact Traveling Wave Solutions of One-Dimensional Parabolic-Parabolic Models of Chemotaxis

By Maria Vladimirovna Shubina

Submitted: December 13th 2019Reviewed: January 15th 2020Published: March 31st 2020

DOI: 10.5772/intechopen.91214

Downloaded: 26

Abstract

In this chapter we consider several different parabolic-parabolic systems of chemotaxis which depend on time and one space coordinate. For these systems we obtain the exact analytical solutions in terms of traveling wave variables. Not all of these solutions are acceptable for biological interpretation, but there are solutions that require detailed analysis. We find this interesting, since chemotaxis is present in the continuous mathematical models of cancer growth and invasion (Anderson, Chaplain, Lolas, et al.) which are described by the systems of reaction–diffusion-taxis partial differential equations, and the obtaining of exact solutions to these systems seems to be a very interesting task, and a more detailed analysis is possible in a future study.

Keywords

  • parabolic-parabolic system
  • exact solution
  • soliton solution
  • Patlak-Keller-Segel model
  • chemotaxis

1. Introduction

This chapter uses the publications of Shubina M.V.:

  1. Exact Traveling Wave Solutions of One-Dimensional Parabolic-Parabolic Models of Chemotaxis, Russian J Math Phys., Maik Nauka/Interperiodica Publishing (Russian Federation), 25(3), 383–395, 2018.

  2. The 1D parabolic-parabolic Patlak-Keller-Segel model of chemotaxis: The particular integrable case and soliton solution, J Math Phys., 57(9), 091501, 2016.

Chemotaxis, or the directed cell (bacteria or other organisms) movement up or down a chemical concentration gradient, plays an important role in many biological and medical fields such as embryogenesis, immunology, cancer growth, and invasion. The macroscopic classical model of chemotaxis was proposed by Patlak in 1953 [1] and by Keller and Segel in the 1970s [2, 3, 4]. Since then, the mathematical modeling of chemotaxis has been widely developed. This model is described by the system of coupled nonlinear partial differential equations. Proceeding from the study of the properties of these equations, it is concluded that the model demonstrates a deep mathematical structure. The survey of Horstmann [5] provides a detailed introduction into the mathematics of the Patlak-Keller-Segel model and summarizes different mathematical results; the detailed reviews also can be found in the textbooks of Suzuki [6] and Perthame [7]. In the review of Hillen and Painter [8], a number of variations of the original Patlak-Keller-Segel model are explored in detail. The authors study their formulation from a biological perspective, summarize key results on their analytical properties, and classify their solution forms [8]. It should be noted that interest in the Patlak-Keller-Segel model does not weaken and new works appear devoted to the study of various properties of equations and their solutions [9, 10, 11, 12] and the links below.

In this chapter we investigate a number of different models describing chemotaxis. The aim of this paper is to obtain exact analytical solutions of these models. For one-dimensional parabolic-parabolic systems under consideration, we present these solutions in explicit form in terms of traveling wave variables. Of course, not all of the solutions obtained can have appropriate biological interpretation since the biological functions must be nonnegative in all domains of definition. However some of these solutions are positive and bounded, and their analysis requires further investigation. Despite the large number of works devoted to the systems under consideration and their properties, as well as the properties of their solutions, it seems to us that the solutions obtained in this paper are new.

The Patlak-Keller-Segel model describes the space–time evolution of a cell density u(t,r)and a concentration of a chemical substance v (t,r). The general form of this model is:

utδ1uη1uϕv=0vtδ22vfuv=0,

where δ1>0and δ20are cell and chemical substance diffusion coefficients, respectively, and η1is a chemotaxis coefficient; when η1>0, this is an attractive chemotaxis (“positive taxis”), and when η1<0, this is a repulsive (“negative”) one [13, 14]. ϕvis the chemosensitivity function, and fuvcharacterizes the chemical growth and degradation. These functions are taken in different forms that correspond to some variations of the original Patlak-Keller-Segel model. We follow the reviews of Hillen and Painter [8] and of Wang [15] and consider the models presented therein.

This paper is concerned with one-dimensional simplified models when the coefficients δ1, δ2, and η1are positive constants, x,t0, u=uxt, and v=vxt.

2. Signal-dependent sensitivity model

Let us start with a model that allows nonnegative bounded solutions that may be of interest from a biological point of view. Now consider the “logistic” model, one of versions of signal-dependent sensitivity model [8] with the chemosensitivity functions ϕv=1+blnv+b, where b=const, and fuv=σ˜uβ˜v. In the review [5] one can see a mathematical analysis of this model. When b=0and β˜=0, the existence of traveling waves was established in [16, 17]. The replacements of tδ1tand uσσ˜δ1ugive δ1=1, α=δ2δ1, β=β˜δ1, and σ=±1. We also set η=η11+bδ1, 1+b>0, as well as ϕv=lnv+b. It should be noted that a sign of σmay effect on the mathematical properties of the system. So, σ=1corresponds to an increase of a chemical substance, proportional to cell density, whereas σ=1corresponds to its decrease. And as we shall see later, various solutions correspond to these two cases.

After the above replacements, the model reads:

utuxx+ηuvxv+bx=0vtαvxxσu+βv=0.E1

If we introduce the function υ=v+b, in terms of traveling wave variable y=xct, where c=const, this system has the form:

uy+cuηulnυy+λ=0αυyy+cυyβυ+βb+σu=0,E2

where u=uy, υ=υy, and λis an integration constant.

In this chapter we will consider the case of λ=0. Then Eq. (2) gives:

u=Cuecyυη,E3

Cuis a constant and we will examine the following equation for υ:

αυyy+cυyβυ+βb+σCuecyυη=0.E4

Since ηis a positive constant, we consider two cases: η=1[Eq. (4) is a linear nonhomogeneous equation] and η1.

  1. η=1

Let us begin with η=1. We introduce the new variable zand the new function w:

z=4σCuαc212ecy2w=4σCuαc2α24αυecy2αE5

and Eq. (4) becomes:

z2wzz+zwz+wz2ν2=Λz1α,E6

where ν2=1α21+4αβc2and Λ=4βbαc24σCuαc214. Eq. (6) is the Lommel differential equation [18, 19] with μ=11α, and we consider σCu>0. Since this is a linear inhomogeneous second-order differential equation, one can integrate it by the method of variation of parameters. We assume a solution in the form:

wz=CJzJνz+CYzYνz,

where Jνzand Yνzare Bessel functions and CJzand CYzare the functions of zthat satisfy the equations:

JνzCJzz+YνzCYzz=0JνzzCJzz+YνzzCYzz=Λz1α.

Considering that Wronskian WJνYνz=2πz, we obtain:

CJz=cJΛπ2z11αYνzdzCYz=cY+Λπ2z11αJνzdz,

where cJand cYare constants. If both of the numbers 1α±νare positive, the lower limits in the integrals may be taken to be zero. Then a particular integral of Lommel equation “proceeding in ascending powers of z” is sμ,νz[19]; if one considers a solution of Lommel equation “in the form of descending series,” one obtains the function Sμ,νz[19] [see Eq. (8)]. Thus, quoting Watson [19] “...and so, of Lommel’s two functions sμ,νzand Sμ,νz, it is frequently more convenient to use the latter.” Then the general solution of Eq. (6) has the form:

wz=CJJνz+CYYνz+ΛSμ,νz,E7

where CJand CYare constants,

Sμ,νz= sμ,νz+2μ1Γμν+12Γμ+ν+12sinπ2μνJνzcosπ2μνYνz,sμ,νz= zμ+1μ+12ν21F21μν+32μ+ν+32z24E8

are Lommel functions, and 1F2is the generalized hypergeometric function [18, 19]. Further, substituting the initial variable yand the function v[see Eq. (5)] into Eq. (7), we obtain a formal solution.

1. b=0

We first consider the case b=0. Then υ=v0and Cu>0. Eq. (6) becomes homogeneous, and for σ=1, its general solution is:

wz=CJJνz+CYYνz.E9

However one can check that the function u=uydiverges as cyfor all ν.

Consider now σ=1. For vyto be real, let α=2. Then Eq. (6) becomes the modified Bessel equation; the analysis of solution behavior at ±leads to suitable solutions for vyand uy:

vy=ecy4Kν2Cuc2ecy2uy=Cue5cy4Kν2Cuc2ecy2E10

with restrictions ν12and β0. So one can see that vy0as cyfor all ν12; vy0for ν<12and vyπ2c28Cu4for ν=12as cyand uy0as y±for all ν12. The curves of these functions are presented in Figures 1 and 2 , and the plots for c=5are thicker than for c=1. Thus, the solution obtained may be considered as a biologically appropriated one, and this requires further investigation.

2. b>0

Figure 1.

(a) v y ; c = 1 ; c = 5 ; C u = 18 . (b) v y ; c = 1 ; c = 5 ; C u = 2 .

Figure 2.

(a) u y ; c = 1 ; C u = 18 . (b) u y ; c = 5 ; C u = 18 .

Let us return to Eq. (6) with Λ0. The analysis of solution asymptotic forms at ±[18, 19] gives the following expressions for vyand uy:

vy+b=4βbαc24σCuαc212αecy2αSμ,ν4σCuαc2ecy2uy=Cu4βbαc24σCuαc212αecy1+12αSμ,ν4σCuαc2ecy2E11

with σCu>0and ν<1α. The latter condition leads to the requirement c24αβ<0. The vyb, uyβbσas cyand vy0and uy0as cy. Thus, one can see that for b>0, σ=1, and Cu>0, uy0is satisfied but vy<0. These functions are presented in Figures 3 and 4 . It should be noted that ν1αor β0because of the pole in Γfunction.

3. b<0

Figure 3.

v y ; c = 1 ; C u = 9 ; σ = 1 ; b = 0.1 .

Figure 4.

(a) u y ; c = 1 ; C u = 9 ; σ = 1 ; b = 0.1 . (b) u y ; c = 1 ; σ = 1 ; b = 0.1 ; α = 1 ; β = − 1 / 4 ; ν = 0 .

Using the analysis of Eq. (11), one can see that the condition b<0along with σ=1and Cu<0(σCu>0) leads to the fact that the function uyhas not changed, but vybecomes positive on all domains of definition. This function is presented in Figure 5 .

  1. η1

Figure 5.

v y ; c = 1 ; C u = − 9 ; σ = − 1 ; b = − 0.1 .

Let us return to Eq. (4) and rewrite it in terms of the variable ξ=ecyα:

ξ2υξξαβc2υ+σαCuc2ξαυη=αβbc2.E12

To integrate this equation, we use the Lie group method of infinitesimal transformations [20]. We find a group invariant of a second prolongation of one-parameter symmetry group vector of (12), and with its help, we transform Eq. (12) into an equation of the first order. It turns out that nontrivial symmetry group requires some conditions:

αβbc2=0,β=α2α+η+1c2αη+32E13

and we consider the case b=0. Thus, υ=v, and for:

z=v1ηαyw=vyvα+η1αE14

we obtain the Abel equation of the second kind:

wz1ηwαz+α+η1z1w2+αzαβc2+σαCuc2zα=0.E15

Then we find the solutions of Eq. (15) in parametric form [21] with the parameter t. Now we consider the case 2α+η1. A combination of substitutions leads to:

z=η+3η+1t2+2σαCuc222α+η1ϑttϑt2αw=z2α2t+22α+η+1η1η+3zα2+α1ηz,E16

where we take

ϑt>0and2α+η1ϑtt<0,E17

and Eq. (15) becomes an equation for the function ϑt. Solving it, for σCu>0, we obtain:

ϑt=C˜ϑ2σαCuc2η+32η+1t2F112η+32η+132η+1c22σαCut2+Cϑ,E18

where C˜ϑand Cϑare constants and 2F1is the hypergeometric Gauss function. Further we obtain the solutions of initial Eqs. (3)(4) in parametric form:

yt=αη+3c2α+η1lnϑtvt=C˜ϑη+322α+η121ηη+1t2+2σαCuc21η+1ϑt2α2α+η1ut=CuC˜ϑη+322α+η121ηη+1t2+2σαCuc21η+1ϑtαη+2α+22α+η1E19

where the constant C˜ϑis chosen so that 2α+η1C˜ϑ<0, which is consistent with Eq. (17). Using the asymptotic representation of hypergeometric Gauss function as t±[18], we can take:

Cϑ>C˜ϑπ2η+12σαCuc21η+1Γ1η+1Γη+32η+1E20

in order for y,v, and uto be real. Then one can see that all functions in Eq. (19) are continuous bounded ones for tand v,uare positive. Hence, one may try to biologically interpret the functions vyand uy, and this requires further investigation. In Figure 6 one may see the different curves vyfor η=0.1and different α. Figure 7 demonstrates vyand uyfor two values η:η=0.1and η=0.01, see Figure 7 . Further, for larger values of αand η, it seems more convenient to present the curves yt, vt, and utto analyze them (see Figures 8 10 ). One can see from Eq. (13) that β0when α2, and the case of β=0and α=2is presented in Figure 11 .

Figure 6.

v y ; η = 0.1 ; σα C u c 2 = 2 ; c = 1 ; C ϑ = 1.4 ; ∣ C ˜ ϑ ∣ = 1 .

Figure 7.

v y ; u y ; α = 0.4 ; σα C u c 2 = 2 ; c = 1 ; C ϑ = 1.35 ; ∣ C ˜ ϑ ∣ = 1 .

Figure 8.

y t ; σα C u c 2 = 2 ; c = 1 ; C ϑ = 1.4 ; ∣ C ˜ ϑ ∣ = 1 .

Figure 9.

v t ; σα C u c 2 = 2 ; c = 1 ; C ϑ = 1.4 ; ∣ C ˜ ϑ ∣ = 1 .

Figure 10.

u t ; σα C u c 2 = 2 ; c = 1 ; C ϑ = 1.4 ; ∣ C ˜ ϑ ∣ = 1 .

Figure 11.

v t ; u t ; α = 2 ; σα C u c 2 = 2 ; c = 1 ; C ϑ = 1.4 ; ∣ C ˜ ϑ ∣ = 1 .

3. Logarithmic sensitivity

The model with logarithmic chemosensitivity function ϕvlnvis also studied. For the case of fuv=vmu+β˜v, where β˜=const, an extensive analysis is performed in [15]. This survey is focused on different aspects of traveling wave solutions. When m=0, this model coincides with Eq. (1) for b=0. When β˜=0and m=1, the system was studied in [22, 23]. The complete analysis for β˜=0is performed in [15]. An existence of global solution is established in [24].

Now we consider the system with ϕv=lnvand fuv=σ˜vuβ˜v. Similarly, a replacement of tδ1tand uσσ˜δ1ugives δ1=1, η=η1δ1, α=δ2δ1, β=β˜δ1, and σ=±1. Then the model has the form:

utuxx+ηuvxvx=0vtαvxxσvu+βv=0.E21

Let us rewrite the system (21) in terms of the function υxt=lnvxt:

utuxx+ηuυxx=0υtαυxxαυx2+βσu=0,E22

Then in terms of the traveling wave variable y=xct, where c=const, Eq. (22) has the form:

uy+cuηuυy+λ=0αυyy+αυy2+cυyβ+σu=0,E23

where u=uy, υ=υy, and λis an integration constant. To integrate Eq. (23), we tested this system on the Painlevé ODE test. One can show that for η>0, it passes this test only if α=2with the additional condition λ=σcβ1+η2[25]. If we express uyas υyfrom Eq. (23), we obtain an equation only for υy; for α=2, it has the form:

2υyyy+3cυyy+c2+ηβυy+22ηυyυyy+22ηυy22ηυy3σλ=0.E24

For λ=σcβ1+η2, this equation can be linearized. It becomes equivalent to the following linear equation for F:

Fy+cF=0,whereFy=e2υ2υyy+cυyηυy2+ηβ2E25

that gives the equation for υy:

2υyy+cυyηυy2+ηβ2=CFe2υcy,E26

where CF=const. If we rewrite Eq. (26) in terms of the variable ξ=ecy2for the function Ψξ=eη2υ, we obtain an equation similar to Eq. (12) with zero right-hand side:

ξ2Ψξξη2β2c2Ψ+ηCFc2ξ2Ψ4η+1=0.E27

Using the result of the symmetry group analysis of Eq. (12), we can write the solution for β=0[see Eq. (19)]:

yt=2clnϑtvt=C˜ϑ22η+2ηt2+2ηCFc21η+2E28

where ϑtis given in Eq. (18) and uymay be expressed from Eq. (23). However one may see that vas t±, and this solution is unacceptable as a biological function.

Another possibility to solve this equation exactly is to put CFequal to zero. When CF=0, that means Fy=0, for β0; Eq. (27) can be linearized by ξ=eτ[21]. Its solution has three forms according to a sign of the expression D=2η2βc2+1. Since vshould be a nonnegative and bounded function as cy±, the only suitable solution is:

vy=ec2ηyCecD4y+C+ecD4y2ηE29

where C±are positive constants and β>0. Unfortunately, the corresponding solution for uyis alternating and has the form:

uy=σc2η+22η2(C21+DecD4y+C+21DecD4y4η2βc2CC+)CecD4y+C+ecD4y2η.E30

It is easy to see that σuyc2η+22η21±Das cy±. These functions are presented in Figures 12 and 13 .

Figure 12.

v y ; c = 1 .

Figure 13.

σu y ; c = 1 .

4. Linear sensitivity

Let us consider the system with linear function ϕvv. When fuv=uv, the system is called the minimal chemotaxis model following the nomenclature of [26]. This model is often considered with fuv=σ˜uβ˜v(σ˜and β˜are constants), and it is studied in many papers. As was proved in [27, 28], the solutions of this system are global and bounded in time for one space dimension. The case of positive σ˜and nonnegative β˜is studied in [29, 30, 31, 32, 33]. As we noted earlier, a sign of σ˜may effect on the mathematical properties of the system, which changes its solvability conditions [34].

Now we consider the linear chemosensitivity function ϕv=vand fuv=σ˜uβ˜v. The replacement of tδ1t, vη1δ1v, and uσσ˜η1δ12uleads to δ1=η1=1, α=δ2δ1, β=β˜δ1, and σ=±1. Then the system has the form:

utuxx+uvxx=0vtαvxx+βvσu=0.E31

This system reduces to the system of ODEs in terms of traveling wave variable y=xct, where c=const:

uy+cuuvy+λ=0αvyy+cvyβv+σu=0,E32

where u=uy, v=vy, and λis an integration constant. As shown in [35], this system passes the Painlevé ODE test only if α=2and β=0. Let us focus on this case.

It is convenient to solve Eq. (32) in terms of variable:

z=κcecy2,E33

where κ>0is an arbitrary constant. Then for vand u, we obtain the solutions in the form:

v=lncκzZν2zu=c2z2114vz2λc,where ν2=14λc3.E34

The function Zνzsatisfies the modified Bessel’s equation and can be present as a linear combination of Infeld’s and Macdonald’s functions.

Using the series expansion of the Infeld’s function, as well as theirs asymptotic behavior [36], one may obtain the following asymptotic forms for evνzand uνz:

z:evνz0;uνz0.E35
z0:evνz,0ν<12;κcC28ππ+22,ν=12;0,ν>12;E36
uνzc2ν12,E37

where the expression for ν=12agrees with Eq. (39).

So, the exact solution obtained has the form:

v=lnecy2A2Iνκcecy2+BKνκcecy22u=σvy2κ2ecy+λc,where ν2=14λc3,E38

where κ>0, A, and Bare arbitrary constants and the functions Iνand Kνare Infeld’s and Macdonald’s functions, respectively. This solution is not satisfactory from the biological point of view, since vyis an alternating function for any ν. However it seems interesting because of the following: in the case of ν=12and B=2+π2πin terms of ecy2, its form coincides with the well-known Korteweg-de Vries soliton.

Consider now the class of solutions with half-integer index ν=n+12, when Zνzcan be expressed in hyperbolic functions. The requirement of absence of divergence ufor finite zleads to the following form for Zn+12z:

Zn+12z=Czn+12dzdzncoshz+ζz,n=2k,Czn+12dzdznsinhz+ζz,n=2k+1;k=0,1;ζ=12ln2π,C=const.E39

At first let us consider the solutions obtained for evn+12and un+12as functions of z. We begin with n=0or ν=12. It is interesting to present the expressions for eu12zand u12z:

ev12z=κC2csech2z+ζE40
u12z=z2c2sech2z+ζ,E41

where Eq. (40) appears the one-soliton solution exactly the same as the well-known one of the Korteweg-de Vries equation. Returning to the variable y:

evecy2=κC2csech2κcecy2+12ln2πuy=σπB1κ2ecysinhκcecy2+π2Beκcecy22.E42

One can see that for σ=1(an increase of a chemical substance), the cell density uy0for B1πand that for B>0uyis the solitary continuous solution vanishing as y±, whereas for B<0uyhas a point of discontinuity. One can say that when B<0, we obtain “blow-up” solution in the sense that it goes to infinity for finite y, and this is true for different ν.

The expressions for n1become more complicated, and one can see the solitonic behavior of evn+12zand the curves for un+12zin Figures 14 and 15 .

Figure 14.

e v n + 1 2 z ; n = 0 , ...6 ; c = 1 .

Figure 15.

u n + 1 2 z ; n = 0 , ...5 ; c = 1 .

The explicit form of our solution in terms of the variable ycan be obtained by direct substitution of Eq. (33) into Eq. (39), where λc=c2nn+1. The resulting formulae are complicated and slightly difficult for analytic analysis; it seems to be more convenient to present the plots.

For n=0in the function ev12y, we have the “step” whose altitude depends on the values of velocity cand arbitrary constant κ. One may see that these curves become higher and shift to the right with different rates for the rising κ. The u12yis the positive function whose altitude and sharpness of peak depend on c (see Figures 16 and 17 ).

Figure 16.

e v n + 1 2 y ; n = 0 .

Figure 17.

u n + 1 2 y ; n = 0 .

For n1we can see that the solitonic behavior of evn+12yis retained for different values of cand κ; the curves become higher and more tight, and they shift to the right also with an increase of cand κ. For the cell density un+12y, the obtained solution has the negative section converging to zero for cy( Figures 18 21 ).

Figure 18.

e v n + 1 2 y ; n = 1 ; 3 ; 5 .

Figure 19.

e v n + 1 2 y ; n = 2 ; 4 ; 6 .

Figure 20.

u n + 1 2 y ; n = 1 ; 3 ; 5 .

Figure 21.

u n + 1 2 y ; n = 2 ; 4 ; 6 .

The curves for the concentration of the chemical substance vn+12yare presented in Figure 22 . Since vn+12yhas to be positive (nonnegative), we see that these functions do not satisfy this requirement in all domains of definition.

Figure 22.

v n + 1 2 y ; n = 0 , 1 , ...6 .

In conclusion it seems interesting to present the plots for evνyand uνyfor different values of ν( Figures 23 25 ). It is interesting to see that there are irregular solutions for evνy; however, the corresponding solutions for uνyare regular [see Eqs. (35)(37)].

Figure 23.

e v ν y ; ν = 1 / 5 ; 8 / 17 ; 7 ; 45 .

Figure 24.

u ν y ; ν = 7 ; 45 .

Figure 25.

u ν y ; ν = 1 / 5 ; 8 / 17 .

5. Conclusion

We investigate three different one-dimensional parabolic-parabolic Patlak-Keller-Segel models. For each of them, we obtain the exact solutions in terms of traveling wave variables. Not all of these solutions are acceptable for biological interpretation, but there are solutions that require detailed analysis. It seems interesting to consider the latter for the experimental values of the parameters and see their correspondence with experiment. This question requires further investigations.

How to cite and reference

Link to this chapter Copy to clipboard

Cite this chapter Copy to clipboard

Maria Vladimirovna Shubina (March 31st 2020). Exact Traveling Wave Solutions of One-Dimensional Parabolic-Parabolic Models of Chemotaxis [Online First], IntechOpen, DOI: 10.5772/intechopen.91214. Available from:

chapter statistics

26total chapter downloads

More statistics for editors and authors

Login to your personal dashboard for more detailed statistics on your publications.

Access personal reporting

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.

More About Us