Co-operation, Competition and Crowding: A Discrete Framework Linking Allee Kinetics, Nonlinear Diffusion, Shocks and Sharp-Fronted Travelling Waves

Invasion processes are ubiquitous throughout cell biology and ecology. During invasion, individuals can become isolated from the bulk population and behave differently. We present a discrete, exclusion-based description of the birth, death and movement of individuals. The model distinguishes between individuals that are part of, or are isolated from, the bulk population by imposing different rates of birth, death and movement. This enables the simulation of various co-operative or competitive mechanisms, where there is either a positive or negative benefit associated with being part of the bulk population, respectively. The mean-field approximation of the discrete process gives rise to 22 different classes of partial differential equation, which can include Allee kinetics and nonlinear diffusion. Here we examine the ability of each class of partial differential equation to support travelling wave solutions and interpret the long time behaviour in terms of the individual-level parameters. For the first time we show that the strong Allee effect and nonlinear diffusion can result in shock-fronted travelling waves. We also demonstrate how differences in group and individual motility rates can influence the persistence of a population and provide conditions for the successful invasion of a population.

For P i m = P g m , P i p = P g p and P i d = P g d = 0, there is no co-operative or competitive mechanism. This gives F (C) = D = D i = D g , and R(C) = λC(1 − C), where λ = λ i = λ g . Therefore, Equation (2) reduces to the Fisher-Kolmogorov equation [24,27,29] (S1.1) As the source term is non-negative for all physical values of C ≥ 0, the agent population will always eventually reach the carrying capacity.
The Fisher-Kolmogorov equation has been studied extensively [18][19][20][21]24,25,27,29]. Here we present the key results in the context of examining the long time travelling wave solution. We seek right moving travelling waves in the co-ordinate z = x − vt, −∞ < z < ∞, where v is a constant wave speed [29]. Transforming Equation (S1.1) into the travelling wave co-ordinate gives With U = dC/dz, Equation (S1.2) can be expressed as a system of ordinary differential equations (ODEs), This system has two equilibrium points: (C, U ) = (0, 0), and (C, U ) = (1, 0). The linear stability of these equilibrium points can be analysed by examining the eigenvalues of the Jacobian at each equilibrium point. At (0, 0) the characteristic equation has solutions ξ = (−v ± √ v 2 − 4λD)/2D, implying that (0, 0) is a stable node provided that v > 2 √ λD, and a stable spiral (focus) otherwise, as λ and D are both positive. We therefore have a minimum wave speed condition, v * = 2 √ λD, that must be satisfied otherwise the solution trajectory will enter non-physical regions of the phase plane [29]. The Jacobian of the linearised system at (1, 0) has eigenvalues ξ = (−v ± √ v 2 + 4λD)/2D, implying that (1, 0) is a saddle point.
The phase plane and associated heteroclinic orbit for Equations (S1.3)-(S1.4) are shown in Figure S1(a). Details of the numerical techniques used to solve Equation (S1.2) and to generate the phase planes are given in the Methods (Main Document). Provided v ≥ 2 √ λD we observe a heteroclinic orbit between (1, 0) and (0, 0). The numerical solution of Equation (S1.2) and the numerical solution of Equation (S1.1), transformed into (C, U ) co-ordinates, are superimposed, showing a good match. This result is unsurprising, as Equation (S1.2) is solved using the minimum wave speed, v = v * = 2 √ λD, and the numerical solution of Equation (S1.1) evolves from a Heaviside initial condition, which is known to approach a travelling wave moving at the minimum wave speed [29]. The numerical solution of Equation (S1.1) at t = 25 and t = 50 is shown in Figure S1(b), confirming that the waveform does not change with time. To quantify the wave speed we calculate the time evolution of the leading edge, L(t) = x f such that C(x f , t) ≈ 1 × 10 −4 . If the solution of Equation (S1.1) forms a travelling wave, L(t) will tend to a straight line with slope v, as t → ∞. In Figure  S1(c), we observe that L(t) is approximately linear with slope v, and hence the solution of Equation (S1.1) moves with approximately constant speed at late time. Overall, these features suggest that the solution of Equation (S1.1) is a travelling wave. Case 2. Different motility rates, equal proliferation rates, no agent death.
If P i m ̸ = P g m the governing PDE contains a nonlinear diffusivity term. Since the agent birth rate is independent of agent type and agents do not die, we consider the same source term as for Case 1. Again, there are no competitive or co-operative mechanisms associated with birth/death but it could be either advantageous (P i m > P g m ) or disadvantageous (P i m < P g m ) for an individual to be isolated from the bulk population, if the goal for the population is to invade unoccupied space. In this parameter regime, Equation (2) simplifies to F (C) has different properties depending on the choice of P i m and P g m . To illustrate this, we present the (P i m , P g m ) parameter space in Figure S2(a), and highlight regions of different behaviour of F (C). If P i m > 4P g m , there will be an interval, 1/3 ≤ α < C < β ≤ 1, centred around C = 2/3, where F (C) < 0. Specifically, this interval is given by .
All parameter pairs that result in F (C) < 0, which we refer to as positive-negative-positive, correspond to the purple region in Figure S2(a), and an example F (C) curve is given in Figure S2(b). Parameter pairs that result in F (C) < 0 with F (1) = 0 correspond to the black line in Figure S2(a), and an example F (C) curve is given in Figure S2(b). We refer to this type of nonlinear diffusivity function as capacity-degenerate positive-negative. It is relevant for us to remark that nonlinear diffusivity functions with negative regions can lead to shocks in the solution of nonlinear diffusion equations without any source term [50,51]. Therefore, it is instructive to consider whether shock-fronted travelling waves exist with Fisher-Kolmogorov kinetics.
For specific parameter regimes, F (C) is degenerate at C = 0, that is, F (0) = R(0) = 0. This type of nonlinear diffusivity function, which we refer to as extinction-degenerate non-negative, leads to sharp-fronted travelling waves, provided that F (C) ≥ 0 for 0 ≤ C ≤ 1 [40,46,47]. For Equation (S2.1), this corresponds to P i m = 0. The parameter pairs that satisfy this condition correspond to the orange line in Figure S2(a), and a typical F (C) curve is given in Figure S2(b). The special case P i m = P g m leads to a constant diffusivity, and parameter pairs that satisfy lie along the cyan line in Figure S2(a). A typical F (C) curve for this case is presented in Figure S2(b). For all other parameter pairs F (C) > 0, which we refer to as strictly positive, and these parameter pairs correspond to the grey region in Figure S2 Positive F (C) (grey), corresponding to P i m = 1 and P g m = 0.5. Negative F (C) for an interval of C (purple), corresponding to P i m = 1 and P g m = 0.1. Negative F (C) for an interval of C with F (1) = 0 (black), corresponding to P i m = 1 and P g m = 0.1. F (0) = 0 (orange), corresponding to P i m = 0 and P g m = 1. F (2/3) = 0 (green), corresponding to P i m = 1 and P g m = 0.25. Constant F (C) (cyan), corresponding to P i m = 1 and P g m = 1. The white circles in (a) denote the parameter pairs used to generate the curves in (b). is shown in Figure S2(b).
We look for a right moving travelling wave solution of Equation (S2.1) in terms of the co-ordinate z = x − vt. Transforming Equation (S2.1) into travelling wave co-ordinates, we obtain Making the substitution U = dC/dz gives dC dz = U, We now consider the properties of the travelling wave solutions for several sub-cases within Case 2. Unlike the Fisher-Kolmogorov equation, the minimum wave speed is unknown and hence all phase planes presented in this section are generated with v obtained from the numerical solution of Equation (S2.1) at sufficiently late time.  Figure S3(a) and Figure S3(d), forms a heteroclinic orbit between (1, 0) and (0, 0). Interestingly, the waveform in Figure S3(e), with P g m > P i m , is relatively sharp near C = 0. If P g m > P i m , F (C) is concave up with a minimum value of P i m /2 at C = 0, for 0 ≤ C ≤ 1, whereas F (C) has a minimum value at C = 2/3 for P i m > P g m . This suggests that F (0) has considerable influence on the waveform.
The observed wave speed in Figure S3(a), v = 0.992, is close to the predicted minimum wave speed v * = 2 √ λD i = 1, whereas the observed wave speed in Figure S3(b), v = 0.584, is greater than the predicted minimum wave speed v * = 0.447. To determine whether v * provides an accurate prediction of the observed wave speed, we calculate the long time numerical solution of Equation (S2.1) and measure v for a suite of P i m and P g m values. Predicted minimum wave speeds and observed wave speeds are compared in Figure S4. In Figure S4(a), the predicted wave speed is accurate for all P g m values with P i m = 0.8. Interestingly, with P i m = 0.2, the predicted wave speed is only accurate for P g m ≤ 0.4. Setting P g m = 1 and varying P i m we observe, in Figure S4(b), that the prediction is accurate for P i m ≥ 0.5. Hence, it appears that for P i m ≥ 2P g Figure S4. Wave speed comparison for Case 2.1. Comparison of the minimum wave speed condition (solid) and the observed wave speed at sufficiently late time (crosses) for (a) constant P g m and a suite of P i m values; (b) constant P i m and a suite of P g m values. All results are obtained using P i p = P g p = 0.5, P i d = P g d = 0, δx = 0.1, δt = 0.01, ϵ = 10 −6 and the Heaviside initial condition.
the minimum wave speed condition is accurate. For P i m < 2P g m the grouped agents may have more successful movement events than the individual agents. Therefore, the dominant contribution to the invasion of the population may be attributed to the grouped agents, which could explain why the minimum wave speed, which depends on P i m , does not provide a good estimate of the observed wave speed in these parameter regimes.
Numerical solutions illustrating travelling wave behaviour for Equation (S2.1) with P i m = 0 are given in Figure S5. In the phase plane for both cases, Figure S5(a) and Figure S5(d), the solution trajectory tends to the origin with dU/dC large and negative. The corresponding numerical solutions of Equation (S2.1), presented in Figure S5(b) and Figure S5(e), approach a travelling wave solution with a sharp front near C = 0. This result is expected as the Heaviside initial condition results in the minimum wave speed that, for a degenerate diffusivity function, results in a sharp-fronted wave [46]. Sub-case 2.3. Positive-negative-positive nonlinear diffusivity function. In order for F (C) to change sign twice, that is, F (C) < 0 for 1/3 ≤ α < C < β ≤ 1 and F (C) ≥ 0 otherwise for 0 ≤ C ≤ 1, the parameters must lie within the purple region in Figure S2(a). In this situation, Equations (S2.4)-(S2.5) are undefined at C = α and C = β, and these singularities cannot be removed using a stretching transformation since R(α) ̸ = 0 and R(β) ̸ = 0. However, it is possible for dU/dz to be finite at C = α and C = β if U α and , (S2. 12) are both finite. This requires the numerator in the expressions (S2.11)-(S2.12) vanish at C = α and C = β, respectively. As such, U α and U β are obtained by solving the system . We note that, as R(C) ≥ 0 for 0 ≤ C ≤ 1, and F ′ (α) ≤ 0 for all possible α values, U α will be real-valued. Subsequently, we have a wave speed condition that v ≥ 2 √ F ′ (β)R(β), as F ′ (β) ≥ 0 for all possible β values. Ferracuti et al. [36] prove that the minimum wave speed, v * , is greater than a threshold value, which in turn is greater than max{R ′ (0)F (0), F ′ (β)R(β)}. Therefore, U β will also always be real-valued.
Applying L'Hopital's Rule to Equation (S2.5), we obtain , (S2. 16) which are finite provided that α ̸ = 2/3 and β ̸ = 2/3. For the system of Equations (S2.4)-(S2.5), we have two straight lines in the phase plane where dU/dz is infinite, at C = α and C = β. These kind of lines have previously been called walls of singularities for hyperbolic models related to chemotactic and haptotactic invasion [52]. For a smooth solution trajectory to exist between two equilibrium points on opposite sides of the wall of singularities, we require that the trajectory passes through the wall of singularities. This implies that the solution trajectory must pass through the wall of singularities at the special points, (α, U α ) and (β, U β ), known as holes in the wall [52,53]. Otherwise, a smooth heteroclinic orbit between (1, 0) and (0, 0) cannot exist, as lim C→α |U | → ∞ and lim C→β |U | → ∞. As U α and U β are real valued and the limits in Equations (20)- (21) are finite, the holes in the wall always exist for Fisher kinetics.
Ferracuti et al. [36] prove that travelling wave solutions exist for reaction-diffusion equations with positivenegative-positive F (C) and Fisher kinetics, however travelling wave profiles arising from the PDE are not presented. An upper bound on the minimum wave speed is stated as [36] (S2. 17) , (S2. 20) and the prime denotes ordinary differentiation with respect to C. Numerical solutions of Equation (S2.1) with P i m > 4P g m are presented in Figure S6. We superimpose the numerical solution of Equation (S2.1) in (C, U ) co-ordinates on the phase plane for the system (S2 in the wall of singularities, denoted using purple circles. Continuum models with negative diffusivity and no source terms have been relatively well studied, and exhibit shock behaviour across the region of negative diffusion [50,51]. Interestingly, our solution does not include a shock and is instead smooth through the region of negative diffusion. Numerical solutions of Equation (S2.1) are presented in Figures S6(b) and S6(e), which appear to take the form of travelling waves. The observed wave speeds, v = 0.864 and v = 0.456, in Figure S6(c) and Figure  S6(f), respectively, are well approximated by the upper bound on the minimum wave speed presented by Ferracuti et al. [36]. The bound provides values for the minimum wave speed of v * = 0.866 and v * = 0.447, respectively. We might expect that the observed wave speeds correspond to the minimum wave speeds since the initial conditions for the numerical solutions are given by a Heaviside initial condition. Sub-case 2.4. Capacity-degenerate positive-negative nonlinear diffusivity function. For the special case where P g m = 0, F (1) = 0. As F (C) is degenerate at C = 1, it is intuitive to expect there could be sharp-fronted travelling wave solutions, with the sharp front near C = 1, similar to the results in Figure S5. However, unlike for the parameter regimes in Figure S5, we have an interval 1/3 < C < 1 where F (C) < 0. To determine whether this negative diffusivity influences the presence of sharp fronts, we follow the approach of Maini et al. [38], who show that the existence of travelling waves for reaction-diffusion equations with capacity-degenerate positive-negative F (C) can be determined by considering the existence of travelling waves for [38]. Equations (S2.21)-(S2.22) have minimum travelling wave speeds v * 0 and v * 1 , respectively. Maini et al. [38] prove that sharp fronts in the travelling wave near The first condition is obviously satisfied, while the second can be determined through linear analysis of Equations (S2.21)-(S2.22) in travelling wave co-ordinates. Both equations have minimum wave speed conditions, v * 0 = v * 1 = 2 √ λD i , to obtain physically-relevant heteroclinic orbits, and hence travelling wave solutions with a sharp region near C = 1 do not exist.
Travelling wave behaviour for a parameter regime with F (1) = 0 is shown in Figure S8. The equilibrium point at (1, 0) is also a hole in the wall. The solution trajectory forms a heteroclinic orbit between (1, 0) and (0, 0), and moves through the region of C where F (C) < 0. Although F (1) = 0, we do not observe a solution trajectory corresponding to a sharp front, as we observed in Figure S5, where F (0) = 0. This result is consistent with the analysis of Maini et al. [38]. The numerical solution of Equation (S8), presented in , has a relatively steep front but is not sharp near C = 1. As L(t), presented in Figure S8(c), becomes linear as t increases and the waveform in Figure S8(b) are consistent, the numerical solution of Equation (S2.1) with F (1) = 0 appears to form a classic travelling wave.
Case 3. Equal motility rates, equal proliferation rates, equal death rates.
For P i d = P g d > 0, with P i m = P g m and P i p = P g p , there are no competitive or co-operative mechanisms. In this case, Equation (2) becomes and, with U = dC/dz, we obtain The source term in Equation (S3.1) is non-positive for all relevant C values if K > λ, and negative for C > (λ − K)/λ otherwise. Hence the population will never reach the original carrying capacity of unity. The new carrying capacity can be determined by considering the zeros of the source term, which occur at C = 0 and C = (λ − K)/λ. Introducing a new variable, C = λC/(λ − K), and rewriting Equation (S3.1) in terms of C we obtain Equation (S3.5) is the Fisher-Kolmogorov equation in terms of the new variable, C, with an intrinsic growth rate (λ − K). As such, the analysis performed for Case 1 is applicable here and we obtain information about the stability of the equilibrium points, as well as the minimum wave speed required for physically meaningful travelling wave solutions. If λ > K, the minimum wave speed is v * = 2 √ (λ − K)D. If K > λ, there is only one physically relevant equilibrium point, C = 0, and hence the population will tend to extinction and travelling wave solutions do not exist.
Travelling wave behaviour for two parameter regimes with λ > K are illustrated in Figure S9. The phase plane for K = 0.1, presented in Figure S9(a), displays qualitatively similar behaviour to the phase plane for K = 0.2, presented in Figure S9(d). Unsurprisingly, the unstable equilibrium point moves closer to zero as K approaches λ. The numerical solutions of Equation (S3.1), presented in Figures S9(b) and S9(e), have significantly different densities behind the wave fronts. However, both travelling wave fronts represent heteroclinic orbits between (C, U ) = ((λ − K)/λ, 0) and (C, U ) = (0, 0). Interestingly, the two travelling wave fronts have approximately the same support, even though the waveform is significantly shallower for the case with K = 0.2. Results in Figures S9(c) and (f) show that both solutions approach travelling waves as t increases, and that increasing the death rate reduces the wave speed. Case 4. Different motility rates, equal proliferation rates, equal death rates.
For P i p = P g p , P i d = P g d and P i m ̸ = P g m , the co-operative or competitive mechanism arises due to the difference in the rate of motility. In this case, the governing PDE is in travelling wave co-ordinates, and with the substitution U = dC/dz, we obtain The system of Equations (S4.3)-(S4.4) has equilibrium points (C, U ) = (0, 0) and (C, U ) = (S, 0), where S = (λ − K)/λ. Increasing K causes a decrease in the carrying capacity, S. If K > λ, the non-zero equilibrium point occurs at a negative C value and hence only one physically relevant equilibrium point exists, implying that the population will become extinct. Hence we only investigate the behaviour of parameter regimes where λ > K.
We introduce the variable C = C/S such that the agent density is scaled by the carrying capacity and, subsequently, the zeros of the source term occur at C = 0 and C = 1. This approach allows us to repeat the analysis for Case 2 with a different F (C). We transform Equation (S4.1) in terms of C to obtain If we define U = dC/dz, Equation (S4.5) corresponds to the system, The transformed nonlinear diffusivity function has different properties depending on S, D i and D g . To highlight this, Figure S10 shows the (P i m , P g m ) parameter space for three different S values and the qualitative behaviour of the corresponding F s (C) function. For S = 1, presented in Figure S10(a), we recover the nonlinear diffusivity function examined for Case 2, Decreasing S to 0.9, presented in Figure S10(b), we observe that the purple region again occurs for for ω < C ≤ 1, and hence F s (C) has only one zero in 0 ≤ C ≤ 1. This type of nonlinear diffusivity function is not observed for Case 2 and we refer to it as positive-negative. Specifically, this behaviour occurs when (16 Furthermore, this implies that for S < 2/3 there are no (P i m , P g m ) values that correspond to positive-negative-positive F s (C). An example of this (P i m , P g m ) parameter space is shown in Figure S10(c). For S < 1/3, F s (C) ≥ 0 for 0 ≤ C ≤ 1. Parameter pairs that correspond to extinction-degenerate nonnegative F s (C) (orange) and constant F s (C) (cyan) exist for all S values.
introducing agent death has the effect of reducing the density sufficiently far behind the wave front, and we observe that the non-zero equilibrium point now occurs at C = S ≤ 1.
Sharp-fronted travelling wave solutions of Equation (S4.1) with extinction-degenerate non-negative F s (C) and S = 0.9 are shown in Figures S11(d)-(f). Introducing agent death does not change the qualitative behaviour compared to the corresponding case with K = 0 (Case 2.2). Specifically, dC/dz approaches C = 0 with a non-zero value and hence the wave front is sharp near C = 0. Again, the density behind the wave front decreases such that C = S, corresponding to the non-zero equilibrium point.
Sub-case 4.3. Positive-negative-positive nonlinear diffusivity function. For positive-negative-positive F s (C), the analysis in Case 2.3 holds provided that λ > K. Specifically, the minimum wave speed condition proved by Ferracuti et al. [36] implies that real-valued holes in the wall will be present for the scaled Fisher-Kolmogorov equation with λ > K. In turn, this suggests that the smooth travelling wave solutions passing through the region of negative diffusivity observed for Case 2.3 will be present with non-zero K. Travelling wave behaviour for Equation (S4.1) with positive-negative-positive F s (C) is demonstrated in Figures S11(g)-(i). The travelling wave solution behaviour is similar to the behaviour in the corresponding case with K = 0 (Case 2.3), with the exception of the reduced carrying capacity.
Sub-case 4.4. Capacity-degenerate positive-negative nonlinear diffusivity function. The capacitydegenerate positive-negative diffusivity case, where F s (1) = R(S) = 0, F s (C) < 0 for ω < C < 1 and F s (C) ≥ 0 otherwise, might be thought to lead to travelling wave solutions with a sharp front near the carrying capacity density [38]. Similar to the approach for Case 2.4, we consider the conditions proposed by Maini et al. [38]. Again, we satisfy the condition that F s (1) = 0. However, the minimum wave speed for the transformation of Equation (S2.1) in 0 ≤ C < ω is the same as the minimum wave speed for the transformation of Equation (S2.1) in ω < C ≤ 1. As such, we do not expect that Equation (S4.1) will approach a travelling wave solution with a sharp front near C = 1. We present travelling wave behaviour for Equation (S4.1) with capacity-degenerate positive-negative F s (C) in Figures S11(j)-(l) and, as anticipated, observe that the travelling wave solution is a classic front.
While there are infinitely many solution trajectories out of the unstable node, we require that the solution trajectory passes through the hole in the wall,  Figure 1 (Main Document) indicate that restricting death events to isolated agents significantly change the behaviour of the agent population. This represents a co-operative mechanism, as there is a benefit to being in close proximity to another agent. In the case where P i d ̸ = 0 and P g d = 0, the source term can be expressed as an Allee effect [15] (S5.1)

Results in
is the intrinsic growth rate, and is the Allee parameter. It follows that Equation (2) becomes .
represents the strong Allee effect, A > 0 [15]. The strong Allee effect has bistable growth kinetics, namely, R(C) < 0 for 0 < C < A and R(C) > 0 for A < C < 1. For low densities there are significantly more isolated agents than grouped agents, which corresponds to negative growth if K i > λ i . If λ i > K i , and λ g > 0, R(C) > 0 for 0 < C < 1. There are two possibilities for this case: r > 0 and r < 0. If r > 0 then A < 0 and hence the growth rate is inhibited at low density, but remains positive, which corresponds to the weak Allee effect [15]. For the case where r < 0 and R(C) > 0, we obtain A > 1 for all parameter combinations. Interestingly, this implies that the growth rate is inhibited at high density, but remains positive. This situation does not correspond to either the weak or strong Allee effect, and we term this behaviour the reverse Allee effect. It is not possible to have a combination of parameters that results in r < 0 and 0 < A < 1 as all of our parameters are non-negative. Representative source terms showing the three types of Allee effect are compared with a logistic source term in Figure S13.
For P i m = P g m , we have linear diffusion in Equation (S5.4). Reaction-diffusion equations with linear diffusion and either weak or strong Allee kinetics have been well-studied [15,19,22,23,25,26,28,30]. We briefly present results here and interpret these in the context of examining the long time travelling wave solution. For additional details we refer the reader to [15]. We look for solutions in the travelling wave co-ordinate z = x − vt. The existence of such solutions has been examined previously and requirements for the initial conditions to converge to a travelling wave solution have been found for both the case where A < 0 and where 0 < A < 1 [23]. Transforming Equation (S5.4) into travelling wave co-ordinates we obtain √ v 2 + 4λ g D)/2D, which implies that the equilibrium point is a saddle point, as λ g is non-negative. A heteroclinic orbit between (1, 0) and (0, 0) exists for a unique wave speed [25]. The equilibrium point at (A, 0) has a characteristic equation with solutions ξ = (−v ± √ v 2 − 4λ g AD)/2D. As we are only concerned with physically realistic equilibrium points, that is, where 0 < A < 1, the equilibrium point (A, 0) will be a stable node provided that the minimum wave speed v > 2 √ λ g AD is satisfied, and a stable spiral otherwise. The spiral behaviour does not cause the solution trajectory to become non-physical and hence this wave speed condition is not required to obtain physically meaningful solutions.
Solutions that display travelling wave behaviour for the weak Allee effect are presented in Figures S14(a)-(c). There is a heteroclinic solution trajectory for Equations (S5.6)-(S5.7) between the two equilibrium points, and the numerical solution of Equation (S5.4) matches the solution trajectory in (C, U ) co-ordinates. The results in Figures S14(b)-(c) suggest that the numerical solution of Equation (S5.4) with the weak Allee effect approaches a travelling wave solution. Since the source term for the reverse Allee effect is qualitatively similar to the weak Allee effect, we do not present solutions here. Solution behaviour for the reverse Allee effect can be found in Figure S15.
A travelling wave solution of Equation (S5.4) in a parameter regime that results in a strong Allee effect is now considered. The phase plane for the system (S5.6)-(S5.7), presented in Figure S14(d), has three physically meaningful equilibrium points, and the equilibrium point at (0, 0) is unstable, unlike in Figure S14(d). However, there is still a heteroclinic orbit between (1, 0) and (0, 0). Unlike for the weak Allee effect, the wave speed that admits this solution trajectory is unique [25]. The numerical solution of Equation (S5.4) shows the solution approaches a travelling wave, although the approach is slower than for the weak Allee effect. This is intuitive, as the growth rate for the weak Allee effect is non-negative, while the strong Allee effect has regions of negative growth.
It is instructive to consider how v depends on A. We calculate the numerical solution of Equation (S5.4) for a range of A values with r = 1 for A ≤ 1, and r = −1 for A > 1, and use the numerical solution to calculate v at sufficiently late time. We consider r < 0 for A > 1 because, due to our parameters being non-negative, A < 1 for r > 0. The minimum wave speed for the travelling wave solution is known for [28]. Consequently, for A > 1/2, the population will tend to extinction because the travelling wave speed is negative [28]. For the case where A > 1 it is unclear whether there is a minimum wave speed condition. A comparison between the observed wave speed for each A value and the predicted minimum wave speed is given in Figure S16. The predicted wave speeds match the observed wave speeds well for A ≤ 1. For A > 1 we superimpose the wave speed prediction v = 2 √ −ArD, and observe that the predictions match the numerical wave speeds well. For the case A > 1 we require that λ i > K i , and hence the minimum wave speed condition is the same as for the weak Allee effect. For P i m ̸ = P g m and P g d = 0, there is a co-operative mechanism in terms of increased survival for agents in close proximity to other agents. In this parameter regime Equation (2) becomes where F (C) = D i (1 − 4C + 3C 2 ) + D g (4C − 3C 2 ). Note that F (C) is the same as in Case 2 and, as such, encodes the same four types of qualitative behaviour. To examine the long term travelling wave behaviour of Equation (S6.1), we transform Equation (S6.1) into travelling wave co-ordinates, z = x − vt, giving It is therefore of interest to determine whether travelling wave solutions can be found for each class of F (C).
Sub-case 6.1. Strictly positive nonlinear diffusivity function. Strictly positive F (C) corresponds to parameters in the grey region of Figure S2 . (S6.7) For v > 0 we require that the transformed source term is, on average, positive, which corresponds to [42][43][44][45] (S6.8) For the case with r > 0, Condition (S6.9) is equivalent to Interestingly, for A > 0, the threshold A value for the population to persist increases if P g m > P i m , and decreases otherwise. For example, if P i m = 0 then A < 6/11 leads to persistence, higher than the threshold A value in the case with constant F (C). Alternatively, as P i m → 4P g m , A → 2/7. This implies that populations where isolated agents are significantly more motile than grouped agents are more susceptible to extinction. This result is intuitive, as the parameter regime considered here describes a co-operative benefit, namely, a reduced death rate for agents in close proximity to other agents. Finally, for the reverse Allee case, where r < 0 and A > 1, Condition (S6.9) is always satisfied and the population persists.
Travelling wave solutions for the strong Allee effect with a strictly positive F (C) are shown in Figure S17. For the strong Allee effect, with parameters that correspond to A = 1/4, presented in Figures S17(a)-(c), we observe a heteroclinic orbit between (1, 0) and (0, 0). The numerical solution of Equation (S6.1) in this parameter regime approaches a travelling wave solution with v > 0. However, if we consider a parameter regime that corresponds to the strong Allee effect with A = 4/9, presented in Figures S17(d)-(f), we observe that, while a heteroclinic orbit between (1, 0) and (0, 0) exists, it corresponds to a negative wave speed. As a consequence, the population tends to extinction in a birth/death parameter regime that would otherwise result in the persistence of the population if the diffusivity is constant. As both the weak and reverse Allee effect are qualitatively similar to Fisher kinetics, numerical solutions are not presented here. Numerical solutions can be found in Figure S18.
Sub-case 6.2. Extinction-degenerate non-negative nonlinear diffusivity function. The case where F (0) = 0 corresponds to parameters along the orange line in Figure S2(a). Sánchez-Garduño and Maini [41] demonstrate that Condition (S6.8) must be satisfied for travelling wave solutions to have v > 0. Furthermore, there is a critical wave speed that results in a sharp-fronted travelling wave [41]. From the results obtained for Case 6.1, Condition (S6.8) is always satisfied for A < 0 or A > 1. For parameter regimes where 0 < A < 1 the choice of P i m and P g m influences whether Condition (S6.8) is satisfied. To obtain an extinction-degenerate diffusivity we require that P i m = 0. Hence (S6.9) implies that for A < 6/11 the wave speed will be positive. To obtain a positive wave speed with constant F (C), we require A < 1/2, which implies that the population is more likely to persist in an parameter regime that leads to extinction-degenerate non-negative F (C). Travelling wave behaviour for the strong Allee effect with extinction-degenerate non-negative F (C) is shown in Figure S19. The numerical solution of Equation (S6.1) with A = 1/4, in Figures S19(a)-(c), leads to a sharp-fronted travelling wave solution near C = 0 with v > 0. With A = 1/4, we expect to obtain v > 0. For a parameter regime that results in A = 4/7, we obtain a travelling wave solution of Equation (S6.1) with v < 0 (Figures S19(d)-(f)). Interestingly, the sharp front near C = 0 is not present for the strong Allee effect with v < 0, unlike with v > 0, where the wave front is smooth. We present travelling wave behaviour for both the weak Allee effect and the reverse Allee effect in Figure S20.
Sub-case 6.3. Positive-negative-positive nonlinear diffusivity function. A positive-negative-positive F (C), where there is an interval α < C < β where F (C) < 0, corresponds to parameter pairs highlighted in purple in Figure S2(a). Kuzmin and Ruggerini [37] examine reaction-diffusion equations with similar properties for the strong Allee effect, in the context of diffusion-aggregation models, and provide conditions for smooth travelling wave solutions to exist. For a solution with v > 0, we require A < α [37] and (S6.10) Furthermore, we require [37] (S6.11) 3 where Φ(y) = 8α 2 y + 4 √ 4α 2 y 2 − 2mα 3 y, , and m = min A suite of P g m values with P i m = 1, which correspond to 1/3 < α < 2/3, are considered for parameter regimes that result in A < α. Figures S21(a)-(c) show the parameter spaces, (A, α), that satisfy Condition (S6.10), Condition (S6.11) and Conditions (S6.10)-(S6.11) simultaneously, respectively. Orange regions represent parameter pairs where the condition is satisfied and grey regions represent parameter pairs where the condition is not satisfied. These results suggest that smooth travelling wave solutions should exist for certain choices of parameters. Interestingly, all parameter pairs that satisfy Condition (S6.10) also satisfy Condition (S6.11).
For Case 2.3 and Case 4.3, smooth travelling wave solutions that pass through holes in the wall of singularities for positive-negative-positive F (C) are obtained. The minimum wave speed bound presented by Ferracuti et al. [36] implies that the location of the holes in the wall occur are real-valued for the wave speed arising from the Heaviside initial condition. As such, to obtain smooth travelling wave solutions of Equation (S6.1) with positive-negative-positive F (C), we might expect that the wave speed satisfies v > 2 √ F ′ (β)R(β), such that the holes in the wall at C = β are real-valued.
Following the approach used for Case 2.3, it is simple to demonstrate that both the weak and reverse Allee effect have real-valued holes in the wall. As such, we observe heteroclinic orbits between (1, 0) and (0, 0) that pass through the holes in the wall, and present the corresponding travelling wave solutions in Figure S22. We now examine numerical solutions of Equation (S6.1) with the strong Allee effect. For parameter regimes that give rise to wave speeds that satisfy v > 2 √ F ′ (β)R(β), numerical travelling wave solutions could not be found. While the condition for real-valued holes in the wall is satisfied, the zeros of Equation (S6.4) are imaginary for a certain interval of C > β. This corresponds to a nullcline that is not real-valued for certain C values.
We now consider parameter regimes corresponding to the strong Allee effect with the additional restriction that v < 2 √ F ′ (C)R(C) for 2/3 < C ≤ 1. For all P i m and P g m that give rise to a positive-negative-positive F (C), holes in the wall at C = β do not exist and, as such, we do not expect to obtain smooth solutions. Interestingly, we observe travelling wave solutions with shocks such that the solution never enters the region α < C < β. An example of a shock-fronted travelling wave solution for the strong Allee effect with both shock-fronted travelling wave solutions arise in other kinds of models, including multispecies models of combustion [54] and haptotactic cell migration [53]. However, the models presented here are very different, and it is therefore of interest to determine the properties of the reaction-diffusion equation that lead to shock-fronted travelling wave solutions.
Sub-case 6.4. Capacity-degenerate positive-negative nonlinear diffusivity function. Capacity-degenerate positive-negative F (C), where F (1) = 0, arises if P g m = 0 and includes an interval 1/3 < C < 1 where F (C) < 0. For Case 2.4, despite the degenerate nature of the nonlinear diffusivity function at C = 1, we did not obtain solutions with a sharp front near C = 1. Instead, the solution passes through the region of negative diffusivity and a hole in the wall at C = 1/3 , leading to smooth travelling wave solutions. As such, we expect similar solutions for both the weak and reverse Allee effect due to the qualitatively similar behaviour of the R(C) function. It is of interest to examine whether smooth or shock-fronted travelling wave solutions arise from Equation (S6.1) for the strong Allee effect, as for the positive-negative-positive diffusivity examined for Case 6.3 no smooth travelling wave solutions could be found. As expected, smooth travelling wave solutions for both the weak and reverse Allee effects with capacitydegenerate positive-negative F (C) are obtained. The solution behaviour for both the weak and reverse Allee effects are presented in Figure S24. For the strong Allee effect, we examine a considerable number of parameter regimes and initial conditions and are unable to find travelling wave solutions. Case 7. Equal motility rates, different death rates.
Without the restriction that only isolated agents are able to undergo death events (P g d ̸ = 0), death events can be considered as either a co-operative mechanism (P i d > P g d ), such as group defence against predation, or a competitive mechanism (P i d < P g d ), where a population is more easily discovered and eradicated, compared to an isolated individual. In these parameter regimes, Equation (2) can be expressed as where If this is not satisfied, R(C) ≤ 0 for 0 ≤ C ≤ 1 and the population will tend to extinction. The corresponding ODE in travelling wave co-ordinates is and, making the substitution U = dC/dz, results in Introducing a new variable C = C/A 1 which, upon substitution into Equation (S7.1), results in . Equation (S7.6) is a reaction-diffusion equation with Allee kinetics in terms of the scaled variable C. Both the carrying capacity and Allee parameter are scaled by A 1 , which influences the maximum population density as well as the threshold density required for positive growth. Following the analysis for Case 5, the minimum wave speed for Equation (S7.6) with Interestingly, this implies that introducing grouped agent death at a rate that does not result in a population tending to extinction has no influence on the invasion speed of the population. Specifically, the condition for A < −1/2 in Case 5 corresponds to 3(λ i − K i ) > λ g . It can be shown that, with 3(λ i − K i ) > λ g , we require 3K g < λ g for A < −1/2. This implies that there is a range of K g values that result in a travelling wave with a minimum wave speed that is independent of both K g and λ g . Interestingly, this suggests that if a control is implemented that increases the death rate of grouped agents, there is a threshold value for the control to influence the invasion speed and the subsequent persistence of the population. Introducing a non-zero K g value for a parameter regime that results in the strong Allee effect with K g = 0 never changes the type of Allee effect. Hence it is possible to go from a weak Allee effect to a reverse Allee effect by introducing a non-zero K g value. Non-zero K g values correspond to a decreased benefit for grouped agents, which explains why the source term, previously a weak Allee effect, becomes the reverse Allee effect, corresponding to inhibited growth at high density.
For the strong Allee effect, corresponding to 0 [28]. This implies that for A 2 > A 1 /2, v < 0 and v > 0 otherwise. Furthermore, the same wave speed applies for −A 1 /2 < A 2 < 0 [28]. For both intervals, the minimum wave speed does depend on the K g value, and hence implementing any kind of partial eradication of the grouped agents will either reduce the speed of invasion or cause the extinction of the population.
Travelling wave behaviour for the weak and strong Allee effect and constant F (C) is shown in Figure S25. For both numerical solutions, calculated with K g = 0.1, the carrying capacity is reduced by approximately 27%. With the exception of K g , the parameters used to obtain the numerical solutions in Figures S25(a)-(c) are the same as in Figures S14(a)-(c) and we observe that, as expected, the wave speed is the same. This demonstrates that, while the carrying capacity is reduced, the population is able to invade vacant space at the same speed, even though a control measure for the grouped agents has been implemented. Results for the reverse Allee effect are presented in Figure S24. Case 8. Different motility rates, different death rates.
Setting P i m ̸ = P g m and P i d ̸ = P g d ̸ = 0 allows for significant flexibility in describing a combination of competitive and/or co-operative mechanisms, depending on the relevant motivation. In this case, Equation (2) can be expressed as . Note that, again, this simplification requires that λ g ≥ 2(K g + √ Kg(K i − λ i )) or λ i > K i , otherwise the population will tend to extinction. In travelling wave co-ordinates, Equation (S8.1) is and, making the substitution U = dC/dz, it corresponds to Introducing the variable C = C/A 1 , Equation (S8.1) can be written as . The transformed nonlinear diffusivity, F A (C), has the same characteristics as F s (C), presented in Figure S10, albeit in terms of the scaled Allee carrying capacity, A 1 . Here we examine the five types of F A (C) for A 1 ̸ = 1.
For the strong Allee effect, A 1 > A 2 = AA 1 , we can determine the threshold value for the persistence of the population, namely, .
Considering the two limiting cases, where D i = 0 and D i = 4D g , A takes on a value of (6A 2 1 − 12A 1 )/(9A 2 1 − 20A 1 ) and (18A 2 1 − 36A 1 + 20)/(27A 2 1 − 60A 1 + 30), respectively. These values reduce to 6/11 and 2/7 in the case that A 1 = 1, as in Case 6.1. To illustrate how the threshold value changes with A 1 , P i m and P g m , Figure S27 shows the maximum A 2 and A values for three different P i m and P g m combinations. The A 2 value corresponds to the persistence threshold for a given A 1 value. The A value can be interpreted as the highest proportion of a given A 1 value that will result in the persistence of the population. For example, in Figure  S27(a), we see that with P i m = 0 and A 1 = 0.5 we require A 2 < 0.194 for persistence. This corresponds to A < 0.388.
Travelling wave behaviour for Equation (S8.1) in a parameter regime corresponding to strictly positive F A (C) and the strong Allee effect is shown in Figures S28(a)-(c). This parameter regime leads to A 1 = 0.723 and A 2 = 0.2764, which is below the persistence threshold value of A 2 = 0.315 for this P i m and P g m combination, and hence the population persists.
Sub-case 8.2. Extinction-degenerate non-negative nonlinear diffusivity function. For extinctiondegenerate F A (C), P i m = 0. As such, the persistence threshold corresponds to (6A 2 1 − 12A 1 )/(9A 2 1 − 20A 1 ). For Case 6.2 we observe that sharp fronts for the strong Allee effect with a extinction-degenerate non-negative For a positive-negative-positive F A (C), there are exactly two zeros at C = α and C = β. In Case 6.3 the strong Allee effect does not give rise to smooth travelling wave solutions, even with real-valued holes in the wall at C = α and C = β. However, interestingly, shock-fronted travelling wave solutions arise from the Heaviside initial condition. Again, we are unable to find numerical travelling wave solutions of Equation (S8.1) in parameter regimes with real-valued holes in the wall. Shock-fronted travelling wave solutions of Equation (S8.1) are given in Figures S28(g)-(i) where the observed wave speed is v = 0.014 < 2 √ min{F ′ (C)R(C)} on the interval 2A 1 /3 < C < A 1 . Smooth travelling wave solutions obtained from the weak and reverse Allee effect are shown in Figure S29.
Sub-case 8.4. Capacity-degenerate positive-negative nonlinear diffusivity function. Capacity-degenerate positive-negative F A (C) requires P i m = 0 and, subsequently, F A (1) = 0. Furthermore F A (C) < 0 for ω < C < S. For Case 6.4 we found smooth travelling wave solutions for both the weak and reverse Allee effect with capacity-degenerate positive-negative F (C) but could not obtain stable solutions for the strong Allee effect. As F A (C) is qualitatively similar to the F (C) considered for Case 6.4 similar results are expected here.  Figure S30) Allee effects are obtained. As for Case 6.4, we consider a variety of parameter regimes corresponding to the strong Allee effect with capacity-degenerate positive-negative F A (C), as well as a number of initial conditions, but are unable to find long time travelling wave-type solutions.
Sub-case 8.5. Positive-negative nonlinear diffusivity function. For the case where F A (C) has exactly one zero on the interval 0 ≤ C ≤ 1 at C = ω, Maini et al. [39] examine the existence of travelling wave solutions, and provide the necessary conditions for existence, (S8.9) where F (ω) = 0 and 0 < ω < 1. For the strong Allee effect in this parameter regime, the third part of Condition (S8.9) corresponds to ,t ≥ 0, on the interval 0 ≤ C < ω, and equivalent to ,t ≥ 0, whereF (C) = −F (1 − C), on the interval ω < C ≤ A 1 . The final necessary and sufficient condition from Maini et al. [39] for the existence of travelling wave solutions is that the minimum wave speed for Equation (S8.11), v * 1 , is greater than, or equal to, the minimum wave speed for Equation (S8.12), v * 2 . On the interval 0 ≤ C < ω, Equation (S8.1) has a strictly positive F A (C), where F A (C) ≤ D i , and strong Allee kinetics. Hence, the minimum wave speed for Equation (S8.11) has an upper bound, v * 1 ≤ √ 2(λ i − K i )D i (1/2 − A 2 ). On the interval ω < C < A 1 Equation (S8.12) has a source term qualitatively similar to the Fisher-Kolmogorov equation and hence a lower bound for the minimum wave speed exists [39], v * 2 ≥ 2 √ −F (A 1 )(λ 2 + 4K g (λ i − λ g − K i + K g )) 1/2 . For all parameter regimes considered that correspond to the strong Allee effect with positive-negative F A (C) we never observe a case where the upper bound for v * 1 is higher than the lower bound for v * 2 and hence the conditions required for travelling wave solutions are not met. As expected, numerical solutions of Equation (S8.1) in these parameter regimes did not lead to travelling wave behaviour. For both the weak and the reverse Allee effect, we expect that solutions do exist as the source terms on both intervals are qualitatively equivalent to a Fisher source term.
Numerical solutions demonstrating the travelling wave behaviour of Equation (S8.1) with positive-negative F A (C) and both the reverse Allee effect and the weak Allee effect are given in Figures S31 and S32, respectively.