==== Front J Math Biol J Math Biol Journal of Mathematical Biology 0303-6812 1432-1416 Springer Berlin Heidelberg Berlin/Heidelberg 33216181 1547 10.1007/s00285-020-01547-1 Article Travelling wave solutions in a negative nonlinear diffusion–reaction model Li Yifei 1 Heijster Peter van peter.vanheijster@wur.nl 12 Marangell Robert 3 Simpson Matthew J. 1 1 grid.1024.70000000089150953School of Mathematical Sciences, Queensland University of Technology, Brisbane, QLD Australia 2 grid.4818.50000 0001 0791 5666Biometris, Wageningen University and Research, Wageningen, The Netherlands 3 grid.1013.30000 0004 1936 834XSchool of Mathematics and Statistics, University of Sydney, Sydney, NSW Australia 20 11 2020 20 11 2020 2020 81 6 1495 1522 21 3 2019 4 2 2020 22 8 2020 © The Author(s) 2020Open AccessThis article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.We use a geometric approach to prove the existence of smooth travelling wave solutions of a nonlinear diffusion–reaction equation with logistic kinetics and a convex nonlinear diffusivity function which changes sign twice in our domain of interest. We determine the minimum wave speed, c∗, and investigate its relation to the spectral stability of a desingularised linear operator associated with the travelling wave solutions. Keywords Nonlinear diffusionTravelling wave solutionsGeometric methodsPhase plane analysisSpectral stabilityMathematics Subject Classification 92C1792D2535K5735B35Wageningen Universityissue-copyright-statement© Springer-Verlag GmbH Germany, part of Springer Nature 2020 ==== Body Introduction Invasion processes have been studied with mathematical models, especially partial differential equations (PDEs), for many years; see, for example, Murray (2002) and references therein. These models describe, for instance, how cells are transported to new areas in which they persist, proliferate, and spread (Mack et al. 2000). To incorporate information about individual-level behaviours in invasion processes, lattice-based discrete models are widely used (Deroulers et al. 2009; Johnston et al. 2017, 2012; Simpson et al. 2010c). In these discrete models, individual agents are permitted to move, proliferate and die on a lattice, and the average density of agents is related to PDE descriptions obtained using truncated Taylor series in the continuum limit (Anguige and Schmeiser 2009; Codling et al. 2008). The macroscopic behaviour described by the PDEs in terms of expected agent density reflects the individual microscopic behaviour. Travelling wave solutions are of particular interest among the macroscopic behaviours arising from these continuum models, as they reflect various modes of microscopic invasive behaviours. One famous model exhibiting travelling wave solutions is the Fisher–KPP equation (KPP refers to Kolmogorov, Petrovsky, Piskunov) proposed in 1937 to study population dynamics with linear diffusion and logistic growth (Fisher 1937; Kolmogorov et al. 1937). The existence and stability of travelling wave solutions of the Fisher–KPP equation has been widely studied, see, for instance, Aronson and Weinberger (1978), Fisher (1937), Harley et al. (2015), Kolmogorov et al. (1937), Larson (1978), Murray (2002) and Sherratt (1998). The Fisher–KPP equation can be derived as a continuum limit of a discrete model under the assumption that the population of cells can be treated as a uniform population without any differences in subpopulations (Bramson et al. 1986). However, differences between individual and collective behaviour have been observed in cell biology and ecology in practice. For instance, in cell biology, isolated cells called leader cells are more motile than the grouped cells, called follower cells (Poujade et al. 2007). Also, contact interactions lead to different motility rates between isolated cells and grouped cells in the migration of breast cancer cells (Simpson et al. 2010c, 2014), glioma cells (Khain et al. 2011), would healing processes (Khain et al. 2007) and the development of the enteric nervous system (Druckenbrod and Epstein 2007). In ecology, the population growth rate of some species decreases as their populations reach small sizes or low densities (Courchamp et al. 1999). This phenomenon is usually referred to as the Allee effect (Allee and Bowen 1932). To describe the invasion process and reflect the difference between collective and individual behaviour, Johnston and coworkers introduced a discrete model considering birth, death and movement events of agents that are isolated or grouped on a simple one-dimensional lattice (Johnston et al. 2017). A discrete conservation statement describing δUj, which is the change of the occupancy of a lattice site j during a time step τ, gives 1 δUj=Pmi2[Uj-1(1-Uj)(1-Uj-2)+Uj+1(1-Uj)(1-Uj+2)-2Uj(1-Uj-1)(1-Uj+1)]+Pmg2[Uj-1(1-Uj)+Uj+1(1-Uj)-Uj(1-Uj-1)-Uj(1-Uj+1)]-Pmg2[Uj-1(1-Uj)(1-Uj-2)+Uj+1(1-Uj)(1-Uj+2)-2Uj(1-Uj-1)(1-Uj+1)]+Ppi2[Uj-1(1-Uj)(1-Uj-2)+Uj+1(1-Uj)(1-Uj+2)]+Ppg2[Uj-1(1-Uj)+Uj+1(1-Uj)]-Ppg2[Uj-1(1-Uj)(1-Uj-2)+Uj+1(1-Uj)(1-Uj+2)]-Pdi[Uj(1-Uj-1)(1-Uj+1)]-PdgUj+Pdg[Uj(1-Uj-1)(1-Uj+1)]. Here, Uj represents the probability that an agent occupies lattice site j, thus, 1-Uj represents the probability that lattice site j is vacant (Simpson et al. 2010a). Pmi and Pmg represents the probability per time step that isolated or grouped agents, respectively, attempt to step to a nearest neighbour lattice site; Ppi and Ppg represents the probability per time step that isolated or grouped agents, respectively, attempt to undergo a proliferation event and deposit a daughter agent at a nearest neighbour lattice site; Pdi and Pdg represents the probability per time step that isolated or grouped agents, respectively, die, and are removed from the lattice. See Fig. 1a for a schematic of the lattice-based discrete model. To obtain a continuous description, Johnston and coworkers treat Uj as a continuous function, U(x, t), and divide (1) by the time step τ. Next, they expand all terms in (1) in a Taylor series around x=jΔ, where Δ is the lattice spacing, and neglect terms of O(Δ3) (Simpson et al. 2010a). As Δ→0 and τ→0 with the ratio Δ2/τ held constant (Codling et al. 2008; Simpson et al. 2010a), they obtain a nonlinear diffusion–reaction equation 2 ∂U∂t=∂∂xD(U)∂U∂x+RU, where 3 DU=Di1-4U+3U2+Dg4U-3U2, is the nonlinear diffusivity function, and 4 RU=λgU1-U+λi-λg-Ki+KgU1-U2-KgU, is the kinetic term. Furthermore, the parameters are given by Dg=limΔ,τ→0PmgΔ22τ,Di=limΔ,τ→0PmiΔ22τ,λg=limτ→0Ppgτ,λi=limτ→0Ppiτ,Kg=limτ→0Pdgτ,Ki=limτ→0Pdiτ, where we require that Ppi,Ppg,Pdi,Pdg are O(τ) (Simpson et al. 2010a). Here, U(x, t) denotes the total density of the agents at position x∈R and time t∈R+; Di≥0 and Dg≥0 are diffusivities of the isolated and grouped agents, respectively; λi≥0 and λg≥0 are the proliferation rates of isolated and grouped agents, respectively; Ki≥0 and Kg≥0 are the death rates of isolated and grouped agents, respectively (Johnston et al. 2017). Note that this particular form (2) was proposed by Johnston et al. (2017). This was one of the first studies that proposed a nonlinear diffusion–reaction model to a mean-field description of a lattice-based stochastic model incorporating agent movement, proliferation and death. Previous work leading to nonlinear diffusion equations only considered the movement of agents and thus did not involve kinetic terms (Johnston et al. 2012; Anguige and Schmeiser 2009).Fig. 1 a One possible time step of the lattice-based discrete model of Johnston et al. (2017): a new grouped agent (agent E) is born and the grouped agent B moves from lattice site 5 to lattice site 4 to become an isolated agent. Pink circles represent isolated agents with birth rate Ppi, death rate Pdi and motility rate related to Pmi; cyan circles represent grouped agents with birth rate Ppg, death rate Pdg and motility rate Pmg. b presents a diffusivity function D(U), given by (3) (cyan curve) satisfying Di>4Dg which makes D(U) change sign twice on (0, 1), and the kinetic term R(U), given by (5) (orange curve) which is positive on (0, 1) and zero at end points U=0 and U=1 (colour figure online) Fig. 2 a The evolution of a Heaviside initial condition to a smooth travelling wave solution obtained by simulating (2) with (3) and (5) with parameters Di=0.25, Dg=0.05 and λ=0.75. We use a finite difference method with space step δx=0.1, time step δt=0.01 and no-flux boundary conditions. Notice that D(U)=0 at α=0.5 and β≈0.83. b The position of the wave L(t), measured by the left-most leading edge point where U is smaller than 10–5, indicating that the solution is travelling at a constant speed c = 0.864. c The wave speed as a function of the initial condition U(x,0)=1/2+tanh-η(x-40)/2. Notice that as η grows to infinity this initial condition limits to the Heaviside initial condition used for the simulation in (a), and the wave speed converges to c≈0.864. The minimum wave speed c∗=2λDi≈0.866 (11) (colour figure online) In this manuscript, we study the effect that aggregation, which is modelled with a nonlinear diffusivity function that goes negative (Simpson et al. 2010b), has on the dynamics of the continuous PDE model. Therefore, we assume that Di>4Dg such that D(U) given by (3) is convex and changes sign twice in our domain of interest (additionally, see Sect. 4.2 for a short discussion related to the other case). For simplicity, we furthermore assume equal proliferation rates, λ=λi=λg, and no agent death, Ki=Kg=0. This way, the kinetic term simplifies to a logistic term 5 RU=λU1-U, and DU has a sign condition: 6 DU>0forU∈0,α∪β,1,DU<0forU∈α,β, where the interval where D(U)<0 is centred at U=2/3, and α,β are given by 7 α=23-Di2+4Dg2-5DiDg3Di-Dg,β=23+Di2+4Dg2-5DiDg3Di-Dg, with 1/3<α<2/3 and 2/3<β<1, see Fig. 1b. That is, we have negative diffusion for U∈(α,β). The relation that Di is larger than Dg indicates that isolated agents are more active than grouped agents, which agrees with the experimental observation that leader cells are more motile than follower cells (Poujade et al. 2007; Simpson et al. 2014). Ferracuti et al. (2009) showed the existence of travelling wave solutions for a range of positive wave speeds for (2) with general convex D(U) that changes sign twice on (0, 1) and R(U) given by (5) based on the comparison method introduced by Aronson and Weinberger (1978). Related studies proved the existence of travelling wave solutions for a similar range of speeds for nonlinear diffusion–reaction equations with different D(U) and different R(U): Malaguti and Marcelli (2003) studied (2) with a logistic kinetic term and a nonlinear diffusivity function satisfying D(0)=0andD(0)>0for allU∈(0,1]. Maini et al. (2006) studied (2) with a logistic kinetic term and a nonlinear diffusivity function satisfying 8 D(U)>0in(0,θ)andD(U)<0inU∈(θ,1), for some given θ∈(0,1) and with D(0)=D(θ)=D(1)=0. In addition, Maini et al. (2007) studied (2) with (8) and a bistable kinetic term satisfying R(0)=R(ϕ)=R(1)=0,R(U)<0inU∈(0,ϕ)andR(U)>0inU∈(ϕ,1). A travelling wave solution of (2) is a solution that travels with constant speed c>0 and constant wave shape, and that asymptotes to 1 as x→-∞ and to 0 as x→∞ (i.e. the roots of R(U)). We only consider positive wave speeds since (2) with (3) and (5) is monostable with a Fisher–KPP imprint, that is, U≡1 is a PDE stable solution of (2), while U≡0 is a PDE unstable solution (in an appropriate function space which will be introduced in Sect. 3). Hence, to study travelling wave solutions we introduce the travelling wave coordinate z=x-ct, where z∈R and c>0, and write (2) in its travelling wave coordinate 9 ∂U∂t=∂∂zD(U)∂U∂z+c∂U∂z+R(U). A travelling wave solution is now a stationary solution to (9), that is, ∂U/∂t=0 (Sandstede 2002). In other words, a travelling wave solution is a solution to the second-order ordinary differential equation (ODE) 10 ddzD(u)dudz+cdudz+R(u)=0, with asymptotic boundary conditions limz→-∞u=1 and limz→∞u=0. In this manuscript, we show the following result: Theorem 1 Model (2) with (3) and (5) and Di>4Dg supports smooth monotone nonnegative travelling wave solutions for 11 c≥2λDi=:c∗. This theorem agrees with the result of Ferracuti et al. (2009), and because of the specific nonlinear diffusivity function, we can further extend their results. Moreover, instead of the comparison method used by Ferracuti et al. (2009), we use a geometric approach to prove the existence of travelling wave solutions. This geometric approach has the advantage that it can also be used to study shock-fronted, discontinuous travelling wave solutions (Wechselberger and Pettet 2010; Harley et al. 2014a, b). While shock-fronted travelling wave solutions are not the focus in this manuscript, we show in the final section that they do exist for (5) with different D(U), see Fig. 10a in Sect. 4.3. The lower bound c∗ in Theorem 1 is often called the minimum wave speed as it represents the monotone nonnegative travelling wave solutions with the lowest wave speed (Murray 2002). Numerical simulations show that (2) with (3) and (5) indeed support smooth travelling wave solutions even though the nonlinear diffusivity function goes negative. Moreover, the speed relates to the initial condition, and the wave speed converges to the minimum wave speed c∗ as the initial condition limits to the Heaviside initial condition, see Fig. 2. We will also show the connection between the existence of smooth monotone nonnegative travelling wave solutions, the spectrum of a desingularised linearised operator associated with the travelling wave solutions, and the minimum wave speed c∗. This manuscript is organised as follows. We prove Theorem 1 in Sect. 2 by using desingularisation techniques (Aronson 1980) and detailed phase plane analysis which have not been applied to (2) before. In Sect. 3, we determine the spectral properties of a desingularised linearised operator associated with the travelling wave solutions and show how the minimum wave speed c∗ is related to absolute instabilities (Sandstede 2002; Kapitula and Promislow 2013; Sherratt et al. 2014). Some interesting results for different nonlinear diffusivity functions with the same kinetic term (5) are discussed in Sect. 4. Here, we also discuss the implications of the analytical results for the discrete model. Note that throughout the manuscript all theoretical results are supported by high-quality numerical simulations of the continuum PDE model. Remark 1 Many essential mathematical questions related to, for instance, well-posedness, remain open for PDEs with forward–backward diffusion, i.e. models like (2) with nonlinear diffusivity functions that change sign. For instance, the well-studied Perona–Malik model (Perona and Malik 1990) from image analysis with forward–backward diffusion, but without a kinetic term, is ill-posed (Weickert 1988). See also Höllig (1983). The ill-posedness of these PDEs with forward–backward diffusion can often be addressed by adding a small regularisation term, like a viscous regularisation term (Novick-Cohen and Pego 1991) or a nonlocal Cahn–Hilliard-type regularisation term (Pego and Penrose 1989). For the Perona–Malik model this was done, with another type of regularisation term, by Barenblatt et al. (1993). Interestingly, different regularisations can have different singular limits, in particular, when shock solutions are formed (see also Sect. 4.3). This is particularly interesting when you realise that most numerical schemes introduce some artificial regularisation. In other words, different numerical schemes can correctly yield different solutions (Witelski 1995). Also, recall that in the derivation of the continuum limit higher order terms were ignored. These higher order terms potentially have a regularising effect and can shed light on the “right” type of regularisation. Since we are constructing smooth solutions in this manuscript, we do not address the question of well-posedness of (2). Existence of travelling wave solutions Transformation and desingularisation We use a dynamical systems approach to analyse the second-order ODE (10) whose solutions that asymptote to limz→-∞u=1 and limz→∞u=0 correspond to travelling wave solutions of (2). Upon introducing p:=D(u)du/dz, (10) can be written as a singular system of first-order ODEs 12 D(u)dudz=p,D(u)dpdz=-cp-D(u)R(u). Travelling wave solutions of (2) now correspond to heteroclinic orbits of (12) connecting (1, 0) to (0, 0). Note that p>0 if du/dz<0 and D(u)<0. Thus, while we expect that the derivative of a travelling wave solution is always negative, p is not necessarily always negative.The nullclines of system (12) are given by p=0 and -cp-D(u)R(u)=0 with the constraint that D(u)≠0. However, D(u) vanishes when u=α and u=β (7), and system (12) is thus undefined, or singular, along the lines u=α and u=β (Simpson and Landman 2007). These lines are sometimes called walls of singularities (Pettet et al. 2000; Wechselberger and Pettet 2010; Harley et al. 2014a). Trajectories can potentially still cross through these walls at special points, sometimes referred to as holes in the wall (Pettet et al. 2000; Wechselberger and Pettet 2010; Harley et al. 2014a), when, in addition to D(u)=0, the right hand sides of the singular system also vanish (and if the holes in the wall are of the correct type (Wechselberger 2005; Wechselberger and Pettet 2010; Harley et al. 2014a)). These holes in the wall, and the trajectories crossing them, can often be linked to folded singularities and canard solutions upon embedding the singular system into higher-dimensional singularly perturbed systems with folded critical manifolds, we refer to Szmolyan and Wechselberger (2001), Wechselberger (2005), Wechselberger and Pettet (2010) and Harley et al. (2014a), and references therein, for more details on this now well-established theory. For system (12) the holes in the wall are (α,0) and (β,0). To remove the singularities, we desingularise system (12) by introducing a stretched variable ξ satisfying D(u)dξ=dz (Aronson 1980; Murray 2002; Sánchez-Garduño and Maini 1994; Harley et al. 2014a). Subsequently, system (12) becomes 13 dudξ=p,dpdξ=-cp-D(u)R(u). Fig. 3 a The phase plane of system (12) with parameters Di=0.25, Dg=0.05, λ=0.75 and c=0.866. The vertical dashed lines are the walls of singularities u=α and u=β and the solid blue lines are nullclines. Red arrows show the orientation of the trajectories. b The phase plane of system (13) for the same parameter values and red lines are nullclines. For u in between α and β, the orientation of the trajectories is opposite compared to (a), while the orientation is the same for u<α and u>β (colour figure online) Here we see that the desingularisation changes the independent variable z in a nonlinear fashion, but it does not change the dependent variables (u, p). Consequently, the (u, p) phase planes of (12) and (13) will have the same trajectories but the “time” it takes to evolve along such a trajectory is different. In particular, when D(u)>0, dξ/dz>0 and therefore trajectories on the phase planes of (12) and (13) have the same orientation. In contrast, when D(u)<0, dξ/dz<0 and trajectories on the two phase planes are in the opposite direction, see Fig. 3. Therefore, heteroclinic orbits of (12) connecting (1, 0) to (0, 0) crossing the holes in the walls (α,0) and (β,0), if they exist, are transformed and separated as heteroclinic orbits connecting (1, 0) to (β,0), (α,0) to (β,0) and (α,0) to (0, 0) of (13) and vice versa. Next, we will prove the existence of these heteroclinic orbits in system (13) for a range of wave speeds c, and then combine these heteroclinic orbits in system (13) as one global heteroclinic orbit in system (12). Phase plane analysis of the desingularised system We first study the desingularised system (13). It has nullclines p=0 and 14 p=-D(u)R(u)c. The intersections of the two nullclines give four equilibrium points: (0,0),(1,0),(α,0),(β,0). Lemma 1 The equilibrium points (1, 0) and (α,0) are saddles. The equilibrium point (0, 0) is a stable node if 15 c≥2D(0)R′(0)=2λDi=c∗, and a stable spiral otherwise. The equilibrium point (β,0) is a stable node if 16 c≥2D′(β)R(β), and a stable spiral otherwise. Proof The Jacobian of system (13) is 17 J(u,p)=01-F(u)-c,whereF(u):=dduD(u)R(u)=D′(u)R(u)+D(u)R′(u), with D(u)R(u) the pointwise product of D(u) and R(u) and where we, as usual, omit the dot. The Jacobian has eigenvalues and eigenvectors λ±=-c±c2-4F(u)2,E±=(1,λ±). For the equilibrium point (1, 0) this reduces to 18 λ1±=-c±c2-4D(1)R′(1)2,E1±=(1,λ1±). The eigenvalues λ1± are real and of opposite sign since D(1)=Dg>0 and R′(1)=-λ<0. Thus (1, 0) is a saddle. Similarly, the Jacobian of the equilibrium point (α,0) has eigenvalues and eigenvectors 19 λα±=-c±c2-4D′(α)R(α)2,Eα±=(1,λα±). Knowing that D′(α)<0 and R(α)>0, λα+ is real and positive and λα- is real and negative. Thus (α,0) is a saddle. The Jacobian of the equilibrium point (0, 0) has eigenvalues and eigenvectors 20 λ0±=-c±c2-4D(0)R′(0)2,E0±=(1,λ0±). The eigenvalues λ0± are real and negative if (15) holds since D(0)=Di>0 and R′(0)=λ>0. Thus the equilibrium point (0, 0) is a stable node if (15) holds. Otherwise, λ0± are complex-valued with negative real parts and (1, 0) is a stable spiral. Similarly, the Jacobian of equilibrium point (β,0) has eigenvalues and eigenvectors 21 λβ±=-c±c2-4D′(β)R(β)2,Eβ±=(1,λβ±). The eigenvalues λβ± are real and negative if (16) holds since D′(β)>0 and R(β)>0. Thus the equilibrium point (β,0) is a stable node if (16) holds. Otherwise, λβ± are complex-valued with negative real parts and (β,0) is a stable spiral. □ Lemma 2 For Di>4Dg, the thresholds of conditions (15) and (16) are ordered as 22 c∗>2D′(β)R(β). Proof The right hand side of (22) is given by 2D′(β)R(β)=23λ(Di-Dg)β(1-β)(β-α). Since c∗=2λDi, proving relation (22) is equivalent to proving Di>3(Di-Dg)β(1-β)(β-α), which is equivalent to proving 23 DiDi-Dg>3β(1-β)(β-α). Knowing that 2/3<β<1 and 0<β-α<2/3 gives 3β(1-β)(β-α)<2/3. Since Di>4Dg, we have that Di/(Di-Dg)>1 since Di>Di-Dg. Hence, (23) holds and thus (22) holds. □ For c0. Similarly, under condition (16), the least negative slope of the stable eigenvectors of (β,0) is λβ+, see (21). This gives λβ+-χ(β)<0. Thus, both eigenvectors have slopes that are more negative than nullcline (14) at (β,0). Therefore, the trajectory moving in (β,0) with decreasing u initially lies under the nullcline (14) for u>β, while they lie above the nullcline (14) for u<β, see also Fig. 4.Fig. 4 A qualitative phase plane of system (13). The three dashed lines are u=α, u=β and u=1. The blue lines are the nullclines p=0 and p=-D(u)R(u)/c. Region R1 is bounded by p=0, u=α and a straight line l1 with negative slope passing through (0, 0). Region R2 is bounded by p=0, u=α and a straight line l2 with negative slope passing through (β,0). Region R3 is bounded by p=0, u=1 and l2 (colour figure online) Next, we consider the region R1 bounded by p=0, u=α and a straight line l1 through (0, 0) with a negative slope μ1. We aim to prove that for c≥c∗, there always exists a slope μ1 so that no trajectories in region R1 can cross through its boundaries. Trajectories starting on p=0 have negative vertical directions since du/dξ=p=0 and dp/dξ=-D(u)R(u)<0 for u∈(0,α). Thus, trajectories in R1 cannot cross through p=0. Trajectories starting on u=α with negative p values point into region R1 since du/dξ=p<0 and dp/dξ=-cp>0. Trajectories starting on l1 satisfy p=μ1u, and they point into R1 only if dpdu|p=μ1u=-c-D(u)R(u)μ1u≤μ1,foru∈(0,α). After rearranging and recalling that μ1<0, we obtain 25 μ1(μ1+c)≤-D(u)R(u)u=-λD(u)(1-u),foru∈(0,α). Lemma 3 For c≥c∗, there exists a μ1 such that inequality (25) is valid for any u∈(0,α). Proof Proving inequality (25) is equivalent to proving 26 μ1(μ1+c)≤-λsupu∈(0,α)D(u)(1-u). The left hand side of inequality (26) is minimal when μ1=-c/2. Setting μ1=-c/2 and substituting into inequality (26) gives a lower bound 27 c1=2λsupu∈(0,α]D(u)(1-u), such that (26) holds for c≥c1. The right hand side of (27) gives 2λsupu∈(0,α)D(u)(1-u)=2λD(0)=2λDi, since D(u) and (1-u) are both decreasing functions on u∈(0,α). Thus, c1=c∗. Hence, for c≥c∗, inequality (26) is valid for μ1=-c/2. □ Knowing that for c≥c∗ inequality (25) is valid, trajectories on l1 with μ1=-c/2 point into region R1. Thus, based on the Poincaré-Bendixson theorem (Jordan and Smith 1999), the observation that the derivative of u is negative in the region R1 (preventing the existence of a homoclinic orbit) and the absence of fixed points in the interior of R1 (preventing the existence of a limit cycle), the trajectory leaving from the equilibrium point (α,0) with decreasing u and decreasing p must connect with the equilibrium point (0, 0) without going negative in u. Similarly, we consider the region R2 bounded by p=0, u=α and a straight line l2 through (β,0) with a negative slope μ2, and the region R3 bounded by p=0, u=1 and l2. Trajectories starting on p=0 have positive vertical directions for u∈(α,β) since du/dξ=p=0 and dp/dξ=-D(u)R(u)>0 and they have negative vertical directions since for u∈(β,1), du/dξ=0 and dp/dξ=-D(u)R(u)<0. Trajectories starting on u=α with positive p point into region R2 since du/dξ=p>0 and dp/dξ=-cp<0. Similarly, trajectories starting on u=1 with negative p point into region R3. In addition, requiring the existence of a slope μ2 such that trajectories starting on l2 point into regions R2 and R3 leads to the condition 28 μ2(μ2+c)≤-D(u)R(u)u-β=-3(Di-Dg)(u-α)R(u),foru∈(α,1). Lemma 4 For c≥c∗, there exists a μ2 such that inequality (28) is valid for any u∈(α,1). Proof The proof of Lemma 4 is analogous to the proof of Lemma 3 and we will omit some of the details. Again, there exists a lower bound c2=23(Di-Dg)supu∈(α,1)(u-α)R(u), such that (28) holds for c≥c2. Next, we show that c223(Di-Dg)supu∈(α,1)(u-α)R(u). This is equivalent to proving Di/(Di-Dg)>3u(1-u)(u-α) for u∈(α,1). Noticing that u-α<2/3, and u(1-u)≤1/4, we obtain 3u(1-u)(u-α)<1/2. Subsequently, we have DiDi-Dg>1>12>3u(1-u)(u-α), since Di>4Dg by assumption. Thus, c22D′(β)R(β), but that only the travelling wave solutions with c≥c∗ (11) have nonnegative densities. The minimal wave speed for the Fisher–KPP equation is closely related to the onset of absolute instabilities.1 Roughly speaking, absolute instabilities imply that perturbations to a travelling wave solution (in an appropriate Sobolev space that will be discussed further on) will grow for all time and at every point in space (Sherratt et al. 2014). These instabilities are related to the absolute spectrum of the linear operator associated with the travelling wave solution and is fully determined by the asymptotic behaviour (z→±∞) of the travelling wave solution (Kapitula and Promislow 2013; Sandstede 2002). Note that the absolute spectrum is, strictly speaking, not part of the spectrum of the linear operator. However, it gives an indication on how far the essential spectrum can be shifted to the left upon using a weighted Sobolev space (Kapitula and Promislow 2013; Sandstede 2002). Consequently, if parts of the absolute spectrum lie in the right half plane, then the essential spectrum cannot be fully weighted into the open left half plane, and the associate solution is hence absolutely unstable.2 The travelling wave solutions of (2) with (3) and (5) as constructed in Sect. 2 asymptote to 0 and 1 and the nonlinear diffusivity function D(U) is positive near U=0 and U=1, see (6). That is, near these points (2) with (3) and (5) has a Fisher–KPP imprint and we therefore expect that the minimal wave speed c∗ of (2) is also closely related to the onset of absolute instabilities. In other words, we expect that the travelling wave solutions of (2) with (3) and (5) are absolutely unstable for 2D′(β)R(β)0 and R′(1)=-λ<0, see Fig. 6. That is, all travelling wave solutions of (2) with (3) and (5) have unweighted essential spectrum in the right half plane.Fig. 6 a The unweighted essential spectrum and the absolute spectrum of the linear operator T~ for c>c∗. The boundary of the unweighted essential spectrum is determined by the dispersion relations of A+ (dashed blue curve) and A- (solid blue curve) and the green region is the interior of the unweighted essential spectrum. The solid red line is the absolute spectrum σabs+ (35), while the dashed red line is the absolute spectrum σabs+ (35). b The unweighted essential spectrum is, for a weight ν=c/(2D(0)) with c≥c∗, shifted to the rightmost boundary of the absolute spectrum σabs+ (colour figure online) From (33) we get that the absolute spectrum at +∞ is given by 35 σabs+=Λ∈R|Λ<-c24D(0)+R′(0)=-c24Di+λ=:K+. Similarly, from (34) we get that the absolute spectrum at -∞ is given by 36 σabs-=Λ∈R|Λ<-c24D(1)+R′(1)=-c24Dg-λ=:K-. That is, σabs- is always fully contained in the open left half plane including the origin, while σabs+ is only fully contained in the open left half plane including the origin for c≥c∗=2λDi, see Fig. 6. The essential spectrum in the weighted space Hν1(R) is determined by the operator Tν(Λ)qs:=D(u^)ddz-A(z;Λ)+D(u^)νIqs=0, or T~ν(Λ)qs:=ddξ-A(ξ;Λ)+D(u^)νIqs=0, see Kapitula and Promislow (2013), and the weighted asymptotic matrices are A+ν(Λ)=A+(Λ)+D(0)νI=D(0)ν1D(0)(Λ-R′(0))-c+D(0)ν, and A-ν(Λ)=A-(Λ)+D(1)νI=D(1)ν1D(1)(Λ-R′(1))-c+D(1)ν. Hence, the boundary of the essential spectrum in the weighted space is given by the dispersion relations Λ+ν=-D(0)k2+i(c-2D(0)ν)k+D(0)ν2-cν+R′(0),Λ-ν=-D(1)k2+i(c-2D(1)ν)k+D(1)ν2-cν+R′(1). These dispersion relations still form two parabolas opening leftward and the intersections with the real axis now depend on ν. We define the intersection of Λ+ν with the real axis as K+ν:=D(0)ν2-cν+R′(0), and the intersection of Λ- on the real axis as K-ν:=D(1)ν2-cν+R′(1). For 2D′(β)R(β)4Dg so that we can obtain a convex nonlinear diffusivity function D(U), given by (3), which changes sign twice in our domain of interest. Furthermore, the assumption of equal proliferation rates and zero death rates leads to a logistic kinetic term R(U), given by (5). The associated numerical simulations of (2) with (3) and (5), see Fig. 2, provided evidence of the existence of smooth monotone travelling wave solutions. To study these travelling wave solutions of (2), we used a travelling wave coordinate z=x-ct and looked for stationary solutions in the moving frame. Consequently, (2) was transformed into the singular second-order ODE (10) which we transformed into a singular system of first-order ODEs (12). To remove the singularities, we used the stretched variable D(u)dξ=dz and transformed (12) into system (13). Next, we analysed the phase plane of the desingularised system (13) and proved the existence of heteroclinic orbits connecting the equilibrium points (0,0),(α,0),(β,0) and (1, 0) for wave speeds c≥c∗, given by (11). Subsequently, based on the relation between the phase planes of (12) and (13), we proved the existence of a heteroclinic orbit in (12) connecting the equilibrium points (1, 0) and (0, 0) passing through (α,0) and (β,0), that are special points on the phase plane called a hole in the wall of singularities. That is, we proved the existence of smooth monotone travelling wave solutions of (2) for c≥c∗. In the end, we showed that the linear operator T~ (32), associated with the travelling wave solutions of (2), with wave speed c4Dg. Notice that c∗=2λDi, hence, the lowest speed for the travelling wave only relates to the diffusivity of individuals and is independent of the diffusivity of the grouped agents. That is, the diffusivity of grouped agents which is smaller than that of isolated agents (Di>4Dg) does not give restrictions for the lowest speed of the moving front. Consequently, we infer that the speed of invasion processes for organisms, for instance, cells, is mainly determined by the behaviour of individuals. Furthermore, the Fisher–KPP equation also has a minimum wave speed for the existence of smooth monotone travelling wave solutions (Kolmogorov et al. 1937; Fife 2013). Hence, a discrete mechanism of invasion processes considering the differences in individual and collective behaviours can lead to a macroscopic behaviour similar to that observed in the discrete mechanism with no differences in isolated and grouped agents.Fig. 7 a D(U) with Di=0.25 and two different Dg. b The corresponding phase planes of system (12) for λ=0.75, c=1, Di=0.25, Dg=0.2 and Dg=0.6, respectively. The two solid curves are the nullclines p=-D(u)R(u)/c with Dg=0.2 (blue curve) and Dg=0.6 (orange curve), respectively. The red dashed lines are the corresponding heteroclinic orbits representing travelling wave solutions in (2) (colour figure online) Smooth travelling wave solutions for positive D(U) If Di<4Dg, then the nonlinear diffusivity function D(U) is positive for U∈[0,1], see Fig. 7a. Thus the corresponding system of first-order ODEs (12) is not singular, and the nullcline p=-D(u)R(u)/c does not cross u-axis, see Fig. 7b. In other words, (0, 0) and (1, 0) are the only equilibrium points. Following the same method as applied in Sect. 2, we obtain the lower bound S1=supu∈(0,1)2D(u)R(u)u=supu∈(0,1)2λ(1-u)D(u), such that there exist smooth monotone travelling wave solutions of (2) for c≥S1. The origin is still a stable node for c≥2λDi:=S2 and S1≥S2. So, if S1≠S2, c≥S1 is only a sufficient condition because there may exist smooth monotone travelling wave solutions of (2) for wave speeds S2≤c