Anisotropic peridynamic model—Formulation and implementation

1.   Introduction

Fear, is defined as, "An unpleasant emotion caused by the belief that someone or something is dangerous" ^[1].

It is a complex emotion, that is critical as a safety measure, and can trigger the "fight or flight" response ^[2] - in particular it can change the way one acts, even when there is no threat present ^[3]. In predator-prey systems, this is most naturally observed among prey, due to their perceived threat of depredation ^[4]. This perception can lead to non-consumptive effects or trait-mediated interactions, which are behavioral, morphological or physiological changes in prey phenotype, due to this threat ^[4,5]. Such effects are known to strongly influence predator-prey dynamics ^[6]. From a mathematical viewpoint, the effects of fear in predator-prey systems has been intensely investigated since the seminal work of Brown et al. ^[7], where optimal foraging theory is extended to consider a game theoretic setup, played out by predator and prey, exhibiting stealth and fear, in which an animal follows a map or a "landscape of fear" (LOF), which describes its predation risk while it navigates the physical landscape. In recent work, Wang et al. ^[8], model fear of depredation, as a (predator) density dependent effect, that negatively effects the prey population. In essence, the prey's growth rate is modeled as a monotonically decreasing function of predator density. Dynamically, a key finding in ^[8] is that under the parametric restrictions of a Hopf bifurcation, an increase in the fear parameter (and prey's birth rate parameter) can alter the direction of a Hopf bifurcation from supercritical to subcritical. Thus, fear enables both supercritical and subcritical Hopf bifurcations, contrary to only the supercritical bifurcations found in classical predator-prey systems. In essence, the fear effect can change the fundamental cyclical patterns of predator-prey dynamics, leading to large scale ecological consequences ^[9]. Note, age/stage structure in the model can also have influence on the dynamics. For example, in ^[10] anti-predator sensitivity, as a consequence of fear, is introduced. Higher values of anti-predator sensitivity may cause predator extinction, thus adult prey persistence. However, lower values lead to co-existence, thus bi-stability is observed.

These results have since initiated a host of activities in diverse ecological scenarios such as when refuges are present ^[11,12] or when the prey has tendencies to avoid predators ^[13]. Predator avoidance could cause its own cascade of dynamics, such as bats avoiding cave predators could fall prey to diurnal birds ^[14], or the green anolis lizard avoiding curly predatory tailed lizards could in turn introduce refuge competition ^[15]. A predator's response could also be influenced by interference pressures, for instance, via a Beddington-DeAngelis functional response ^[16]. Various works have considered the fear effect in case of group defense by the prey ^[10,17]. It has been investigated in the context of cooperative and competitive systems within the larger predator-prey context. These include the fear effect when predators are cooperating, as earlier mentioned ^[14,18] in the hunting process, or when they are hunting for competing prey ^[19]. These effects have been investigated in the three and multi-species settings as well ^[20,21] where fear can damp population explosions ^[22]. Various authors have considered the fear effect in a stochastic setting ^[23] as well as a spatially explicit setting, in the context of taxis type movements, as well as pattern formation ^[24,25]. It can also lead to chaotic dynamics ^[26]. However, the effect of fear has been far less investigated in classical monotone systems, such as purely cooperative or competitive two species systems - that are outside the predator-prey setting.

Competition among two species, typically modeled via the Lotka–Volterra competition model and its variants have been intensely investigated in the last few decades. These models take into account growth and inter/intraspecific competition ^[27], and predict well-observed states in biology of co-existence, competitive exclusion of one competitor, and bi-stability, and find diverse applications in ecology and invasion science ^{[28,29,30,31,32]}. There are several ecological motivations for competitors being fearful of each other. This is perhaps most naturally seen to occur with intraguild predation - a widespread phenomenon in many food webs, where competitors will kill and consume each other ^[33]. Recent evidence of non-consumptive effects exerted by intraguild predator mites (Blattisocius dentriticus) on their competitor (Neoseiulus cucumeris) show this can be an important factor in determining food web dynamics in biological control ^[34,35]. However, there is strong evidence for fear in purely competitive two species systems without predatory effects. Barred owls (Strix varia) are a species of owl, native to eastern North America. They have expanded their range westward over the last century and are considered invasive in western North America. Currently, their range overlaps with the spotted owl (Strix occidentalis), which is native to the north west and western North America. This has resulted in intense competition between the two species ^[36]. Barred owls exert a strong negative influence on spotted owls, threatening their possible competitive exclusion ^[37]. Field observations report frequent barred owl attacks on spotted owls, and even on surveyors imitating spotted owl calls ^[38]. There is also evidence of barred owls aggressively chasing spotted owls out of shared habitat - but not the opposite ^[39]. There is also other very recent empirical evidence to support such investigations. In ^[15] a series of 6 year long experiments are conducted in various Caribbean islands that aim to refute the theory of adaptive predation - which suggests that predators reduce dominant competitors, thus preventing competitive exclusion and enhancing coexistence in food webs. However, non-consumptive effects such as fear of depredation can have strongly influencing effects ^{[15,40,41,42]}. ^[15] considers a series of experiments with two competing species of lizards, brown anolis (Anolis sagrei) that dwells on tree trunks, and green anolis (Anolis smaragdinus) that dwells on tree canopies. The experiments show that typically these species co-exist - due to a clear niche separation. However, the introduction of an intraguild predator, the curly tailed lizard (Leiocephalus carinatus) that dwells on the ground, causes (non-consumptive) fear driven effects. The brown anolis being fearful of possible depredation (as the lower half of tree trunks are within striking distance of the curly tailed lizard) moves upwards into the canopy, which is occupied by the green anolis. Herein interspecific competition intensifies leading to a loss of co-existence. However, what is most crucial in this study, is that the fecal analysis of the curly tail lizard shows that its diet included the brown anolis, in only 2 out of 51 samples examined. Thus the new dispersal pattern of the brown anolis and the "refuge competition" is driven strongly by a non-consumptive fear effect and not a consumptive one. Thus the brown anolis and the curly tail lizards are really competitors, as they have a strong overlap in dietary niche for several insects. However, the brown anolis is clearly fearful of the curly tailed lizard ^[15]. There is further evidence of non consumptive effects such as fear among competing aphid species, as well as competitors that feed on aphids ^[43,44]. Such interplay between competition and predation has been investigated ^[45], where it is proposed that in many ecological processes, competition and predation are interlinked, and depending on niche overlap, one of them will dominate to drive the underlying dynamics. Such evidence clearly motivates considering fear dynamics into competitive two species model where one of the competitors is fearful of the other. Furthermore, one could look at more complex settings as well, consisting of several interacting predators and prey, across several trophic levels ^[46].

There are also several socio-economic-political settings, where pure competitors may be fearful of each other. Small/new businesses (as weaker competitors) may be fearful of large businesses (that are stronger/dominant competitors), due to their already large market share ^[47]. But large business may also be fearful of small local businesses, due to their familiarity with local nuances, that may yield competitive advantage at a small local scale ^[48] - Thus there could be strong spatial heterogeneity associated to the way in which competitors are fearful of each other. Fear is also conceivable among two competing political parties, where the weaker party on a national scale, may have a stronger voter bank at a regional scale ^[49] - This again is a clear example of spatial heterogeneity, in that a weaker competitor overall is dominant or exerts fear on a stronger competitor at a local/regional scale. Or perhaps two warring drug cartels, where the weaker cartel has certain local/territorial strongholds ^[50,51] - within which they might be able to induce fear among the stronger cartel ^[50]. Such phenomenon becomes even more interesting in the spatially explicit case where this fear could be heterogeneous in the spatial domain of interest. This connects back to the LOF concept, where the fear function is essentially the map that describes how the fear levels change as a species disperses over a physical landscape.

2.   Model formulation

2.1.   Recap of basic competition model

Consider the classical two species Lotka-Volterra ODE competition model,

where $u$ and $v$ are the population densities of two competing species, $a_1$ and $a_2$ are the intrinsic (per capita) growth rates, $b_1$ and $b_2$ are the intraspecific competition rates, $c_1$ and $c_2$ are the interspecific competition rates. All parameters considered are positive. The dynamics of this system are well studied ^[52]. We recap these briefly,

● $E_0 = (0, 0)$ is always unstable.

● $E_u = (\frac{a_1}{b_1}, 0)$ is globally asymptotically stable if $\dfrac{a_{1}}{a_{2}} > \max\left\lbrace\dfrac{b_{1}}{c_{2}}, \dfrac{c_{1}}{b_{2}}\right\rbrace$. Herein $u$ is said to competitively exclude $v$.

● $E_v = (0, \frac{a_2}{b_2})$ is globally asymptotically stable if $\dfrac{a_{1}}{a_{2}} < \min\left\lbrace\dfrac{b_{1}}{c_{2}}, \dfrac{c_{1}}{b_{2}}\right\rbrace$. Herein $v$ is said to competitively exclude $u$.

● $E^* = \Big(\frac{a_1b_2-a_2c_1}{b_1b_2-c_1c_2}, \frac{a_2b_1-a_1c_2}{b_1b_2-c_1c_2}\Big)$ exists when $b_1b_2-c_1c_2 \neq 0$. The positivity of the equilibrium holds if $\frac{c_2}{b_1} < \frac{a_2}{a_1} < \frac{b_2}{c_1}$ and is globally asymptotically stable if $b_1b_2-c_1c_2 > 0$. This is said to be the case of weak competition.

● If $b_1b_2-c_1c_2 < 0$, then $E^* = \Big(\frac{a_1b_2-a_2c_1}{b_1b_2-c_1c_2}, \frac{a_2b_1-a_1c_2}{b_1b_2-c_1c_2}\Big)$ is unstable as a saddle. In this setting, one has initial condition dependent attraction to either $E_u(\frac{a_1}{b_1}, 0)$ or $E_v(0, \frac{a_2}{b_2})$. This is the case of strong competition.

We proceed by considering the effects of fear on the classical model (2.1), when one of the competitors is fearful of the other.

2.2.   The fear model for two competitors

Species diffusion is ubiquitous in spatial ecology ^[53]. Species disperse to find mates, food and shelter ^[54]. Such movement is modeled often via reaction diffusion systems ^[53]. The spatially explicit Lotka-Volterra model, particularly in the case of heterogeneity in spatial resources has been intensely investigated, ^{[27,55,56,57,58,59,60,61,62,63,64]}. Herein, we consider the spatially explicit version of the Lotka-Volterra competition model for two competing species $u, v$, where the competitor $v$ is fearful of $u$. We model the fear effect as in ^[8], where the growth rate of the fearful competitor $v$, is not constant but rather density dependent. Essentially, the growth rate is decreased by a factor $\approx \frac{1}{1+k u}$, where $k \geq 0$ is a fear coefficient. Thus a higher density of the competitor $u$ increases the fear in $v$. When $k = 0$, the assumption is there is no fear and one recovers the classical spatially explicit model. If fear is present, we obtain the following spatially explicit model,

Here $u(x, t), v(x, t)$ are the densities of two competing species, where $v$ is also fearful of $u$. The species diffuse in a bounded domain $\Omega \subset \mathbb{R}^{n}$, with dispersal speeds $d_{1}$ and $d_{2}$, respectively. We prescribe Neumann boundary conditions, modeling the effect that the species do not immigrate or emigrate from the domain $\Omega$. Also the initial data is chosen to be positive. We now proceed to study the dynamics of the above model, when various forms of fear are considered. We divide the dynamical analysis into two components, $d_{1} = d_{2} = 0$, the ODE case discussed next, after which we consider the PDE case $d_{1} > 0, d_{2} > 0$, where we allow a spatially heterogenous fear function as well, so $k = k(x)$.

3.   The ODE case

3.1.   The case of $v$ fearing $u$

We consider the dynamics of (2.2), when $d_{1} = d_{2} = 0$, and we are in a bonafide ODE setting, where the competitor $v$ is fearful of $u$,

3.1.1.   Existence

The nullclines associated with the system (3.1) are

Hence, the boundary equilibrium points are obtained by substituting $u = 0$ and $v = 0$ in the above equations of the nullclines, respectively. Denote the boundary equilibrium points as $\widehat{E}_1 = (0, 0)$, $\widehat{E}_2 = (\frac{a_1}{b_1}, 0)$ and $\widehat{E}_3 = (0, \frac{a_2}{b_2})$.

For the interior equilibrium, substitute $u^* = \frac{a_1}{b_1} - \frac{c_1}{b_1} v^*$ in the second nullcline equation, i.e.,

On simplification, we have that $v^*$ solves a quadratic equation of the form $A_1(v^*)^2+B_1v^*+C_1 = 0$, where

Let

be the two roots of above quadratic equation. Define discriminant as $\Delta_1 = B_1^2-4A_1C_1$. WLOG assume $v_1^* < v_2^*.$ Moreover, consider the following parametric restriction

We can prove the existence of a positive equilibrium point $\widehat{E}_4$ with the choice of specific parameters. Let us use Descartes's rule of sign to establish some sufficient conditions for the existence of one or two positive equilibrium points.

Two positive equilibrium points: Under the assumption $A_1 > 0, B_1 < 0, C_1 > 0$, i.e., $b_2b_1 > 2c_1c_2 > c_1c_2, k < \dfrac{1}{a_1^2c_2} \Big(b_1^2a_2-a_1c_2b_1\Big)$ and $\Delta_1 > 0$, we have two positive roots. In order to claim that these two roots correspond to two positive interior equilibria, we need some extra assumption given by:

Hence, if $b_1(a_1 c_2-a_2 b_1)+a_1^2 c_2 k < [(b_2 b_1-c_1 c_2)(2 a_1k-b_1)-a_1kc_1 c_2] \frac{a_1}{c_1},$ we have two positive interior equilibrium points $\widehat{E}_4 = (u_i^*, v_i^*)$ for $i = 1, 2$.

Remark 1. The existence of two positive internal equilibrium points can also be observed under the assumption $A_1 > 0, B_1 < 0, C_1 > 0$, and $\Delta_1 > 0$, which ensures that

and there must be two positive internal equilibria. Similarly, under the assumption $A_1 < 0, B_1 > 0, C_1 < 0$, and $\Delta_1 > 0$, we have the existence of two positive internal equilibria.

One positive equilibrium point: Under the assumption $A_1 > 0, B_1 < 0$, i.e., $b_2b_1 > 2c_1c_2 > c_1c_2$, we have at least one positive root of the quadratic equation. If $C_1 < 0$, which is $k > \dfrac{1}{a_1^2c_2} \Big(b_1^2a_2-a_1c_2b_1\Big)$ along with $(3.4)$ gives existence of one positive equilibrium point $\widehat{E}_4 = (u^*, v^*)$. Moreover, if $C_1 > 0$, which is $k < \dfrac{1}{a_1^2c_2} \Big(b_1^2a_2-a_1c_2b_1\Big)$, $\Delta_1 > 0$, and $v_1^* < \frac{a_1}{c_1}$ and $v_2^* > \frac{a_1}{c_1}$, then we have existence of one positive equilibrium point $\widehat{E}_4 = (u^*, v^*)$.

We formulate all these restrictions as an existence theorem:

Theorem 3.1. For the given ODE system (3.1), we always have three boundary equilibrium points, namely $\widehat{E}_1: = (0, 0), \widehat{E}_2: = (\frac{a_1}{b_1}, 0)$ and $\widehat{E}_3: = (0, \frac{a_2}{b_2})$. For the case of two positive interior equilibrium points $\widehat{E}_4 = \Big(u_i^*, v_i^*\Big)_{i = 1, 2}$, we have the following parametric restrictions:

Lastly, for the case of one positive interior equilibrium point $\widehat{E}_4 = (u^*, v^*)$, we have either one of the following choices of parametric restrictions:

1) $b_2b_1 > 2c_1c_2 > c_1c_2, \quad k > \dfrac{1}{a_1^2c_2} \Big(b_1^2a_2-a_1c_2b_1\Big)$ and $(3.4)$ holds true.

2) $b_2b_1 > 2c_1c_2 > c_1c_2, \quad k < \dfrac{1}{a_1^2c_2} \Big(b_1^2a_2-a_1c_2b_1\Big)$, $\Delta_1 > 0$ and $v_1^* < \frac{a_1}{c_1}$ and $v_2^* > \frac{a_1}{c_1}$.

3) $b_2b_1 < c_1c_2 < 2c_1c_2, \quad k > \dfrac{1}{a_1^2c_2} \Big(b_1^2a_2-a_1c_2b_1\Big)$, $\Delta_1 > 0$ and $v_1^* < \frac{a_1}{c_1}$ and $v_2^* > \frac{a_1}{c_1}$,

where $v_i^*$ are the roots of the quadratic equation defined as $(3.3).$

We now provide several lemmas, so that we can compare the effect of fear to the classical competition case.

Lemma 3.2. Consider the given ODE system (3.1), with $k = 0$, such that the system is in the weak competition setting with $b_2b_1 > 2c_1c_2$. Then for a fear coefficient $k$ such that $k > k_{c} = \dfrac{1}{a_1^2c_2} \Big(b_1^2a_2-a_1c_2b_1\Big)$, $(\frac{a_1}{b_1}, 0)$ is globally asymptotically stable. That is, $u$ will competitively exclude $v$.

Remark 2. For the existence of one positive interior equilibrium point for the given system (3.1), we have either one of the following choices of parametric restrictions: For weak competition,

and for strong competition,

3.1.2.   Linear stability analysis

We next perform stability analysis on the equilibrium points of system (3.1). The Jacobian matrix of system (3.1) is given by

We state the following lemmas:

Lemma 3.3. $\widehat{E}_1$ is unstable.

Proof. On evaluating (3.7) at $\widehat{E}_1$, we have

Being a triangular matrix, we know that the above matrix has two positive eigenvalues $a_1$ and $a_2$. Hence, the equilibrium point $\widehat{E}_1$ is unstable.

Lemma 3.4. $\widehat{E}_2$ is locally stable iff $k > \dfrac{a_2 b_1^2-c_2a_1b_1}{a_1^2c_2}$.

Proof. We again evaluate (3.7) at $\widehat{E}_2$ and obtain

Being a triangular matrix, the above matrix has two eigenvalues, $\lambda_1 = -a_1$ and $\lambda_2 = \frac{a_2b_1}{b_1 + k a_1} -\frac{c_2 a_1}{b_1}$. As $\lambda_1$ is always negative, if we can show that $\lambda_2$ is negative, we are done. We make the assumption that,

Therefore, the boundary equilibrium point $\widehat{E}_2$ is locally stable.

Local stability of $\widehat{E}_2$ actually implies global stability, we can see this via a simple geometric argument.

Lemma 3.5. $\widehat{E}_2$ is globally stable if $k > k_{c} = \dfrac{a_2 b_1^2-c_2a_1b_1}{a_1^2c_2}$.

Proof. Consider the nullclines of $u$ and $v$, where,

In order to establish the global stability of $\widehat{E}_2$, via the geometry of the nullclines, it suffices to show that

when $u = \frac{a_1}{b_{1}}$, i.e., when $k > \frac{a_2 b_1^2-c_2a_1b_1}{a_1^2c_2}$. Herein, the $v$-nullcline lies completely below the $u$-nullcline, and via the convexity of the $v$-nullcline, it lies completely below the straight line connecting its $v$ and $u$ intercepts - which lies completely below the $u$-nullcline. Now, via the standard theory of competition and a comparison argument, where $v$ is compared to the $\tilde{v}$ that is a solution to the straight line nullcline connecting the $v$ and $u$ intercepts of the $v$-nullcline, we have the global stability of $\widehat{E}_2$.

Lemma 3.6. The equilibrium point $\widehat{E}_3$ is locally stable iff $a_1b_2 < c_1a_2$.

Proof. Similar evaluation of (3.7) at $\widehat{E}_3$ yields

Being a triangular matrix, the above matrix has two eigenvalues, $\lambda_1 = -a_1$ and $\lambda_2 = a_1- \dfrac{c_1a_2}{b_2}$. Under the assumed parametric restriction,

Hence, the equilibrium point $\widehat{E}_3$ is locally stable.

Remark 3. This is essentially the parametric restriction for the case when $v$ excludes $u$, see Figure 4. Thus under this parametric restriction, $\widehat{E}_3$ is locally stable for all $k$, but globally stable only for certain $k$ values.

Lemma 3.7. The interior equilibrium $\widehat{E}_4$ exists and is locally stable if $0 < k < \frac{1}{a_{2}} \left(\frac{b_{1}b_{2}}{c_{1}} - c_{2}\right).$

Proof. On evaluating (3.7) again at $\widehat{E}_4$, we have

For the local stability of $\widehat{E}_4$, it is enough to show that $\text{Tr}(J(\widehat{E}_4)) < 0$ and $\det(J(\widehat{E}_4)) > 0.$ Simple computations yield

and

Note that

Therefore, if $k$ is chosen such that, $0 < k < \frac{1}{a_{2}} \left(\frac{b_{1}b_{2}}{c_{1}} - c_{2}\right)$, then,

and the result follows.

Lemma 3.8. The interior equilibrium $\widehat{E}_4$ exists and is a saddle if

Proof. In order to claim that the interior equilibrium $\widehat{E}_4$ is a saddle, it is enough to show that $\text{Tr}(J(\widehat{E}_4)) < 0$ and $\det(J(\widehat{E}_4)) < 0.$ We have that

and

Under the assumption and density of reals, we have

Hence, $\widehat{E}_4$ is a saddle.

Lemma 3.9. Consider the system (3.1). For any given $k > 0$ such that, there exist two positive interior equilibria, a sink and a saddle.

Proof. Theorem 3.1 provides conditions under which two positive interior equilibria exists. Since stability has to alternate by standard theory ^[65] for planar systems, one of the equilibrium is stable while the other is unstable. Instability as a source is impossible due to the lack of periodic dynamics in the system via Lemma 3.10. Thus the unstable equilibrium must be a saddle. This proves the lemma.

Lemma 3.10. Consider the ODE system $(3.1)$. There do not exist any periodic orbits for the system, for any values of the fear parameter $k$.

Proof. Consider the function $\phi(u, v) = \frac{1}{u v}$ where $u$ and $v$ are both non-zero. Let,

Then we have

The result follows by application of the Dulac criterion ^[65].

Corollary 1. Consider the ODE system $(3.1)$, with the fear term modeled as $\dfrac{a_2 v}{1+kF(u)}$, where $F(u) \in C^{0, 1}$. There do not exist any periodic orbits for the system, for any values of the fear parameter $k$.

Proof. The result follows via direct application of the Dulac function from Lemma 3.10.

Remark 4. Some observations:

$1)$ From Figure 4, we observe that when there is no fear, then $(0, v^*)$ is globally stable. For a sufficiently large level of fear in species $v$, a bi-stability situation is created. That is, for a certain initial data, species $u$ is completely excluded by $v$ and initial data is attracted to the co-existence state (See Figure 4).

$2)$ If $(0, v^*)$ is globally attracting, a much higher level of fear in species $v$ ($\approx 200$) can change the dynamics to a strong competition type case. However most data in this setting is attracted to $(u^{*}, v^{*})$. For initial data $(u_{0}, v_{0})$ to be attracted to $(0, v^*)$, we would need $v_{0} \quad \gg \quad 1, u_{0} \quad \ll \quad 1$.

$3)$ Lemma 3.10 shows the lack of periodic orbits in (3.1), which could be interpreted as the inability of the fear factor to destabilize a steady state through a Hopf bifurcation to cyclical dynamics. This result could be extended to a wider class of fear models if fear in the second competitor was modeled as $\frac{a_{2}v}{1+kF(u)}$, where $F(u)$ was a general density dependent fear effect, not necessarily linear. See Corollary 1 which shows non existence of periodic orbits. Thus the result can be extended to a more general class of models than (3.1).

$4)$ One could consider the case of $u$ fearing $v$, but since the fear model is symmetric, most of the results will be similar. Thus we do not consider this case in detail, rather relegate a few technicalities to the appendix Section B, that rigorously show that the case of $u$ fearing $v$ and $v$ fearing $u$, are indeed models that are symmetric.

3.2.   Bifurcation Analysis

A bifurcation is said to occur in a dynamical system when the behavior of solutions changes when a parameter is varied. Bifurcation analysis is useful in understanding and measuring these qualitative changes as the system switches from stable to unstable and vice-versa.

3.2.1.   Saddle-node bifurcation

The following theorem is connected to the existence of a saddle-node bifurcation for the growth rate $a_1$ when there is no fear effect in competitor $u$.

Theorem 3.11. The ODE system $(3.1)$ undergoes a saddle-node bifurcation around $\widehat{E_4}$ at $a_1 = a_1^* = \frac{2 (k u^*+1)^2 \left(b_1 u^* \left(2 b_2 v^*+c_2 u^*\right)+b_2 c_1 v^*{^2}\right)-a_2 \left(2 b_1 u^* (k u^*+1)+c_1 v^* (2 k u^*+1)\right)}{(k u^*+1) \left((k u^*+1) \left(2 b_2 v^*+c_2 u^*\right)-a_2\right)}$ and when $c_1 > \frac{a_2\left(2b_2v^*+c_2u^* \right)}{u^* v^* \left(ka_2+\left(1+ku^* \right)^2c_2 \right)}$.

Proof. We shall use the Sotomayor's Theorem ^[65] to show the occurrence of a saddle-node bifurcation at $a_1 = a_1^*$. At $a_1 = a_1^*$, we can have $\det(J) = 0$ and for $c_1 > \frac{a_2\left(2b_2v^*+c_2u^* \right)}{u^* v^* \left(ka_2+\left(1+ku^* \right)^2c_2 \right)}$, $\text{Tr}(J) < 0$. This implies that, $\det(J)$ admits a zero eigenvalue. Now let $P = (p_1, p_2)^T$ and $Q = (q_1, q_2)^T$ be the eigenvectors of $J$ and $J^{T}$ corresponding to the zero eigenvalue respectively.

We have that, $P = \left(\dfrac{-B}{A}, 1 \right)^T$ and $Q = \left(\dfrac{-A}{F}, 1 \right)^T$

where $A = a_1-2b_1 u^*-c_1 v^*$, $B = -c_1u^*$, $F = -\dfrac{a_2 k v^*}{(1+ku^*)^2}-c_2 v^*$.

Now, let $X = (X_1, X_2)^T$ where

Furthermore,

and

Therefore by the Sotomayor's theorem system (3.1) undergoes a saddle-node bifurcation at $a_1 = a_1^*$ around $\widehat{E_4}$.

Remark 5. Consider the case of two interior equilibria, such as in Figure 4. Then decreasing the fear coefficient $k$ results in a saddle-node bifurcation, where the interior equilibrium (the saddle $\widehat{E_{4}}$ and the node $\widehat{E_{5}}$) collide and disappear, resulting in the boundary equilibrium $\widehat{E_{2}} = (u^{*}, 0)$ becoming globally asymptotically stable. Thus in this setting, a certain critical level of "fear" can maintain a co-existence state, but fear less than this level takes the system back to a competitive exclusion type scenario. This transition occurs via a saddle-node bifurcation. This can be rigorously proven by adopting the methods of Theorem 3.11, to the parameter $k$ instead of using $a_{1}$.

3.2.2.   Transcritical bifurcation

Theorem 3.12. The model $(3.1)$ undergoes a transcritical bifurcation around $\widehat{E}_3$ when $a_1 = a_1^* = \dfrac{c_1a_2}{b_2}$ and $c_1 \neq \dfrac{b_1b_2}{ka_2+c_2}$.

Proof. The Jacobian matrix for system (3.1) evaluated at $\widehat{E}_3$ with $a_1^* = \dfrac{c_1a_2}{b_2}$ is given as

The corresponding eigenvalues to the Jacobian of (3.1) in Eq (3.8) are $\lambda_1 = 0$ and $\lambda_2 = -a_2.$ Clearly, there is a zero eigenvalue at $a_1 = a_1^* = \dfrac{c_1a_2}{b_2}$. Next, we let $W = (w_1, w_2)^T$ and $Z = (z_1, z_2)^T$ represent the eigenvectors related to the zero eigenvalue of the matrices $J$ and $J^{T}$ respectively.

We obtain $W = \left(-\dfrac{b_2}{ka_2+c_2}, 1 \right)^T$ and $Z = \left(1, 0 \right)^T$. Now, let $X = (X_1, X_2)^T$ as already defined. Presently, we validate the transversality conditions using the Sotomayor's theorem ^[65]. Now,

Also,

and

Therefore by the Sotomayor's theorem system (3.1) experiences a transcritical bifurcation at $a_1 = a_1^* = \dfrac{c_1a_2}{b_2}$ around $\widehat{E}_3$.

Remark 6. Consider the case of one interior equilibrium, such as in Figure 2-that is, we are in the weak competition case, when there is no fear or $k = 0$. Now increasing the fear coefficient $k$, results in a transcritical bifurcation, where the interior equilibrium (the node $\widehat{E}_{4}$) and the boundary saddle equilibrium $\widehat{E}_{2} = (u^{*}, 0)$ collide, exchange stability, after which $\widehat{E}_{4}$ now moves to the 4$^{th}$ quadrant, while $\widehat{E}_{2}$ becomes globally asymptotically stable. Thus in this setting, a certain critical level of "fear" can move the system from a weak competition setting of coexistence to a competitive exclusion type scenario. This transition occurs via a transcritical bifurcation. This is rigorously proven in Theorem 3.12.

4.   The PDE Case

In this section we consider the dynamics of (2.2), when $d_{1}, d_{2} > 0$, and we are in a bonafide PDE setting, where the competitor $v$ is fearful of $u$.

4.1.   Notations and preliminary observations

To prove global existence of solutions to (2.2), it suffices to derive uniform estimate on the $\mathbb{L}^{p}$ norms of the R.H.S. of (2.2), for some $p > \frac{n}{2}$. Classical theory will then yield global existence ^[66]. The usual norms in spaces $\mathbb{L}^{p}(\Omega)$, $\mathbb{L}^{\infty }(\Omega)$ and $\mathbb{C}\left(\overline{\Omega }\right)$ are respectively denoted by

To this end, we use standard techniques ^[67]. We first recall classical results guaranteeing non-negativity of solutions, local and global existence ^[67,68]:

Lemma 4.1. Let us consider the following $m\times m$ - reaction diffusion system: for all $i = 1, ..., m,$

where $d_i \in(0, +\infty)$, $f = (f_1, ..., f_m):\mathbb{R}^m \rightarrow \mathbb{R}^m$ is $C^1(\Omega)$ and $u_{i0}\in L^{\infty}(\Omega)$. Then there exists a $T > 0$ and a unique classical solution of (4.2) on $[0, T).$ If $T^*$ denotes the greatest of these $T's$, then

If the nonlinearity $(f_i)_{1\leq i\leq m}$ is moreover quasi-positive, which means

then

Lemma 4.2. Using the same notations and hypotheses as in Lemma 4.1, suppose moreover that $f$ has at most polynomial growth and that there exists $\mathbf{b}\in \mathbb{R}^m$ and a lower triangular invertible matrix $P$ with nonnegative entries such that

Then, for $u_0 \in L^{\infty}(\Omega, \mathbb{R}_+^m),$ the system (4.2) has a strong global solution.

Under these assumptions, the following local existence result is well known, see D. Henry ^[66].

Theorem 4.3. The system (4.4) admits a unique, classical solution $(u, v)$ on $[0, T_{\max }]\times \Omega$. If $T_{\max } < \infty$ then

where $T_{\max }$ denotes the eventual blow-up time in $\mathbb{L}^{\infty }(\Omega).$

The next result follows from the application of standard theory ^[69].

Theorem 4.4. Consider the reaction diffusion system (2.2). For spatially homogenous initial data $u_{0} \equiv c, v_{0} \equiv d$, with $c, d > 0$, then the dynamics of (4.4) and its resulting kinetic (ODE) system, when $d_{1} = d_{2} = 0$ in (2.2), are equivalent.

4.2.   Spatially heterogeneous fear

Our objective now is to consider the case of a fear function that may be heterogeneous in space. A motivation for this comes from several ecological and sociological settings. For example it is very common for prey to be highly fearful closer to a predators lair, but less fearful in a region of refuge ^[12], or in regions of high density due to group defense ^[70]. Furthermore, a conceivably weaker drug cartel, could have certain localized strongholds, within which they would be more feared by stronger groups. To these ends, it is conceivable that the fear coefficient $k$ is not a constant, but actually varies in the spatial domain $\Omega$, so $k = k(x)$, which could take different forms depending on the application at hand. This is also in line with the LOF concept ^[7]. Thus we consider the following spatially explicit version of (3.1), with heterogeneous fear function $k(x)$, resulting in the following reaction diffusion system,

where $\Omega \subset \mathbb{R}^{n}$. We assume Neumann boundary conditions. Also, we prescribe spatially homogeneous (flat) initial conditions $u(x, 0) = u_0 (x) \equiv c > 0, \quad v(x, 0) = v_0(x) \equiv d > 0.$ Furthermore, we impose the following restrictions on the fear function $k(x)$,

Remark 7. If $k(x) \equiv 0$ on $\Omega_{1} \subset \Omega$, with $|\Omega_{1}| > \delta > 0$, or $k(x) \equiv 0$ on $\cup^{n}_{i = 1}\Omega_{i} \subset \Omega$, with $\Sigma^{n}_{i = 1}|\Omega_{i}| > \delta > 0$, that is, on non-trivial parts of the domain, the analysis is notoriously difficult, as one now is dealing with a degenerate problem. See ^[71,72] for results on this problem. This case is not in the scope of the current manuscript.

Since the nonlinear right hand side of (4.4) is continuously differentiable on $\mathbb{R}^{+}\times$ $\mathbb{R}^{+}$, then for any initial data in $\mathbb{C}\left(\overline{\Omega }\right)$ or $\mathbb{L} ^{p}(\Omega), \; p\in \left(1, +\infty \right)$, it is standard to estimate the $\mathbb{L}^{p}-$norms of the solutions and thus deduce global existence. Standard theory will apply even in the case of a bonafide fear function $k(x)$, because due to our assumptions on the form of $k$, standard comparison arguments will apply ^[73]. Thus applying the classical methods above, via Theorem 4.3, and Lemmas 4.1 and 4.2, we can state the following lemmas:

Lemma 4.5. Consider the reaction diffusion system (4.4), for $k(x)$ such that the assumtions via (4.5) hold. Then solutions to (4.4) are non-negative, as long as they initiate from positive initial conditions.

Lemma 4.6. Consider the reaction diffusion system (4.4). For $k(x)$ such that the assumtions via (4.5) hold, the solutions to (4.4) are classical, that is for $(u_{0}, v_{0}) \in \mathbb{L}^{\infty }(\Omega)$, $(u, v) \in C^{1}(0, T; C^{2}(\Omega))$, $\forall T$.

Our goal in this section is to investigate the dynamics of (4.4). Herein we will use the comparison technique, and compare to the ODE cases of classical competition, or the constant fear function case, where the dynamics are well known.

Remark 8. The analysis in this section are primarily focused on the choice of spatially homogenous (flat) initial data.

We begin by defining the following systems of PDEs,

where

We assume Neumann boundary conditions for all of the reaction diffusion systems (4.6)–(4.9). Also in each of the systems we prescribe spatially homogenous (flat) initial conditions $u(x, 0) = u_0 (x) \equiv c > \; 0, \quad v(x, 0) = v_0(x) \equiv d > 0.$

We now state the following lemma:

Lemma 4.7. Consider the reaction diffusion system (4.4), as well as the reaction diffusion systems (4.6)–(4.9). Then the following point wise comparison holds,

Proof. Note via positivity of solutions to (4.4), (4.6)–(4.9), the definitions via (4.10), and the upper bound on $u$, which is a solution to (4.4) of $\frac{a_{1}}{b_{1}}$, (derived via comparison to the logistic equation), we have

Thus the result follows via standard comparison theory ^[73].

4.3.   The competitive exclusion case

Theorem 4.8. Consider the reaction diffusion system (4.4), for a fear function $k(x)$, such that the assumptions via (4.5) are met, and

Then the solution $(u, v)$ to (4.4) converges uniformly to the spatially homogenous state $(\frac{a_{1}}{b_{1}}, 0)$ as $t \to \infty$.

Proof. From the classical theory of competition ^[52], we know the dynamics for (4.6), that is in the competitive exclusion case, when

we have

Moreover, under the assumption

and making use of Lemma 3.5, along with the use of Theorem 4.4, we have $(\widetilde{u}, \widetilde{v}) \to \Big(\dfrac{a_1}{b_1}, 0\Big)$. Now using Lemma 4.7 we have,

which entails,

subsequently,

Now using a squeezing argument, in the limit that $t \rightarrow \infty$, we have uniform convergence of solutions of (4.4), i.e.,

as $t \rightarrow \infty$.

Using the positivity of solutions the requirement on $k(x)$ and so in turn on $\mathbf{\widetilde{k} }$, can be weakened to derive a stronger result.

Theorem 4.9. Consider the reaction diffusion system (4.4), for a fear function $k(x)$, such that the assumptions via (4.5) are met, and $\dfrac{a_1}{a_2} > \max \Big\{ \dfrac{b_1}{c_2}, \dfrac{c_1}{b_2} \Big\}$. Then the solution $(u, v)$ to (4.4) converges uniformly to the spatially homogeneous state $(\frac{a_{1}}{b_{1}}, 0)$ as $t \to \infty$.

Proof. From the classical theory of competition ^[52], we know the dynamics for (4.6), that is in the competitive exclusion case, when

we have

Now using Lemma 4.7, we have

Using the non negativity of solutions to (4.4) via Lemma 4.5, entails,

subsequently,

Theorem 4.10. Consider the reaction diffusion system (4.4), for a fear function $k(x)$, such that the assumptions via (4.5) are met, and

Then the solution $(u, v)$ to (4.4) converges uniformly to the spatially homogeneous state $(0, \frac{a_{2}}{b_{2}})$ as $t \to \infty$.

Proof. From the classical theory of competition ^[52], we know the dynamics for (4.6), that is in the competitive exclusion case, when

we have

Consider the nullclines of (4.8). If

it follows via the geometry of the nullclines, that the $v$-nullcline remains above the $u$-nullcline, and thus $(\widetilde{u}, \widetilde{v}) \to \Big(0, \dfrac{a_2}{b_2}\Big)$. On simplification, and using the upper bound estimate for $u$, we have

Now using Lemma 4.7 we have,

which entails,

subsequently,

Now using a squeezing argument, in the limit that $t \rightarrow \infty$, we have uniform convergence of solutions of (4.4), i.e.,

as $t \rightarrow \infty$.

4.4.   The weak competition case

We state the following result:

Lemma 4.11. Consider the reaction diffusion system (4.4), for a fear function $k(x)$, such that the assumptions via (4.5) are met, with

Then the solution $(u, v)$ to (4.4) converges uniformly to the spatially homogeneous state $(u^*, 0)$ as $t \to \infty$.

Proof. Via Lemma 3.2, and the parametric restrictions assumed we have,

and the result follows via similar analysis as in Theorem 4.8.

Lemma 4.12. Consider the reaction diffusion system (4.4), such that system is in the weak competition case when $k(x) \equiv 0$. Then given $0 < \epsilon \quad \ll \quad 1$, there exists a fear function $k_{\epsilon}(x)$, for which the assumptions via (4.5) are met, such that the solution $(u, v)$ to (4.4) with the fear function $k_{\epsilon}(x)$, converges uniformly to a spatially homogeneous state $(u^*, v^{*})$ as $t \to \infty$.

Proof. Given $0 < \epsilon \quad \ll \quad 1$, we can always construct a $k_{\epsilon}$, s.t $k_{c} - \epsilon \leq \widehat{k}_{\epsilon}$, whereas $\widetilde{k}_{\epsilon} \leq k_{c} + \epsilon$. Thus via Lemma 4.7 we have,

Lemma 3.7 ensures we have $(\widetilde{u}, \widetilde{v}) \to (u^*, v^{*})$ and $(\widehat{u}, \widehat{v}) \to (u^{**}, v^{**})$, where the spatially homogeneous solutions may be different. Hence, via squeezing argument, we can take $\epsilon \rightarrow 0$, to yield the uniform convergence of solutions, i.e.,

This proves the lemma.

The numerical simulations above motivate the following conjecture:

Conjecture 1. Consider the reaction diffusion system (4.4), for a fear function $k(x)$, such that the assumptions via (4.5) are met, with

and the parameters follow $(3.5)$ and Theorem $3.1$, then the solution $(u, v)$ to (4.4) converges uniformly to the spatially homogenous state $(u^*, v^{*})$ as $t \to \infty$.

Remark 9. We see via Lemma 3.5, and Lemma 3.7, that in the ODE case, if we are in the weak competition setting without fear, then a critical amount of fear $k_{c}$ is both sufficient and necessarily required to change the system's dynamics to a competitive exclusion type scenario. In the PDE case, where the fear function $k(x)$ can be spatially heterogeneous, this requirement is certainly sufficient, as seen via Theorem 4.13, but not necessary, in a point wise sense.

This result is stated next:

Theorem 4.13. Consider the reaction diffusion system (4.4), with initial data $u_{0}(x) \geq 1$, s.t $u \nearrow u^{*}$($u$ increases to $u^{*})$, such that system is in the weak competition setting when $k(x) \equiv 0$. For a fear function $k(x)$, s.t the assumptions via (4.5) are met, the solution $(u, v)$ converges uniformly to the spatially homogeneous state $(\frac{a_{1}}{b_{1}}, 0)$ as $t \to \infty$, if the following condition holds,

where $k_{c}$ is as defined in Lemma 3.5.

Proof. Consider (4.4), integrating the equation over $\Omega$ yields,

This follows via the mean value theorem for integrals. We can now compare,

where $\tilde{v}$ solves $\frac{d \tilde{v}}{dt} = \frac{a_{2}}{1+k_{c}\frac{a_{1}}{b_{1}}}\tilde{v} - b_2 \tilde{v}^2 - c_2 u\tilde{v}.$ Thus using Theorem 4.4, Lemma 3.2 and the positivity of solutions, this entails

The bounds on $v$, standard Lebesgue convergence theorems, and squeezing argument entail,

which implies the uniform convergence,

Remark 10. The $C_{1}$ in Theorem 4.13 is a pure constant, that will depend on the size of the domain $\Omega$, and the other problem parameters, but not on the initial data or the spatial variable $x$.

Remark 11. Clearly $k(x)$ could be chosen such that it lies below $k_{c}$ for a portion of the domain, and above $k_{c}$ on some portion of the domain - thus $\min_{x\in \Omega} k(x)$ does not lie uniformly above $k_{c}$. Yet via Theorem 4.13, we see that, one can change the system's dynamics and bring it to a competitive exclusion type scenario, from a coexistence situation.

The next result gives a lower estimate for the fear function $k$.

Lemma 4.14. Consider the fear function $k(x)$ in Theorem 4.13, then we have,

Proof. We have that for any fear function $k(x)$ satisfying (4.5),

Now using Theorem 4.13 the result follows.

4.5.   The strong competition case

Theorem 4.15. Consider the reaction diffusion system (4.4), for a fear function $k(x)$, such that the assumptions via (4.5) are met, with $b_{1}b_{2} < c_{1}c_{2}$. Then there exists sufficiently large positive initial data $[u_0, v_0]$ for which the solution $(u, v)$ to (4.4) converges uniformly to the spatially homogenous state $(0, \frac{a_{2}}{b_{2}})$ as $t \to \infty$, while there also exists sufficiently small positive initial data $[u_1, v_1]$ for which solution $(u, v)$ to (4.4) converges uniformly to the spatially homogenous state $(\frac{a_{1}}{b_{1}}, 0)$ as $t \to \infty$.

Proof. Consider the system (4.6). From the classical strong competition parametric restrictions, $b_{1}b_{2} < c_{1}c_{2}$ that are assumed, as well as Lemma 4.7, we can make use of standard competition theory and use the stable manifold theorem, i.e., $\exists W_{s}(E_{4}) \in C^{1}$ separatrix, such that for initial data $(\overline{u}_0, \overline{v}_0)$ chosen above $W_{s}(E_{4})$ the solution $(\overline{u}, \overline{v}) \to (0, v^*)$ and for initial data chosen below $W_{s}(E_{4})$, $(\overline{u}, \overline{v}) \to (u^*, 0)$. Here $E_{4}$ is the interior saddle equilibrium to the kinetic (ODE) system for (4.6). Moreover, since $\dfrac{a_{2}}{1 + \mathbf{\widetilde{k}} \frac{a_{1}}{b_{1}}} \le a_{2}$, and $b_{1}b_{2} < c_{1}c_{2}$, we have that for (4.10) we still remain in the strong competition case, and via standard theory again, $W^{1}_{s}(E^{*}_{4}) \in C^{1}$ separatrix, such that for initial data $(\widehat{u}_0, \widehat{v}_0)$ chosen above $W^{1}_{s}(E^{*}_{4})$, the solution $(\tilde{u}, \tilde{v}) \to (0, v^*)$ and for initial data chosen below $W^{1}_{s}(E_{4})$, $(\tilde{u}, \tilde{v})\to (u^*, 0)$. Here $E^{*}_{4}$ is the interior saddle equilibrium to the kinetic (ODE) system for (4.9). Now since $\dfrac{a_{2}}{1 + \mathbf{\widetilde{k}} \frac{a_{1}}{b_{1}}} \le a_{2}$, the $v$ component of $E^{*}_{4}$ is higher than the $v$ component of $E_{4}$. Now using the $C^{1}$ property of the separatricies $W^{1}_{s}(E^{*}_{4}), W_{s}(E_{4})$, we have the existence of a wedge $V$ emanating from $E_{4}$, s.t within $V$ we have $W^{1}_{s}(E^{*}_{4}) \geq W_{s}(E_{4})$. Note via Lemma 4.7, we have $\tilde{v} \leq v \leq \overline{v}$. Let us consider positive initial data $(u_0, v_0)$ chosen large enough, within $V$ s.t. $(u_0, v_0) > W^{1}_{s}(E^{*}_{4}) > W_{s}(E_{4})$, we will have

On the other hand, for positive initial data $(u_1, v_1)$ chosen small enough via an analogous construction we will have

This proves the theorem.

From the observations in the above numerical simulations, we can state a conjecture concerning positive interior equilibrium.

Conjecture 2. Consider the system $(4.4)$. For any non-negative fear function $k(x)$ such that

and parametric restrictions given by (4.5) and Theorem 3.1 hold true, for $\mathbf{k} = \mathbf{\widetilde{k}}, \mathbf{\hat{k}}$. There exists some data $[u_0, v_0]$, for which the solution $(u, v) \to (u^*, v^*)$, and for some choice of data $[u_1, v_1]$, the solution converges to the boundary equilibrium $(0, v^*)$.

5.   Numerical simulation

All of our PDE simulations for the case of a spatially heterogeneous fear function were performed in $\verb|MATLAB|$ R2021b, using the $\verb|MATLAB|$ inbuilt function pdepe, which is used in solving $1-$D parabolic and elliptic PDEs. The simulations were run on an $8-$core CPU, Apple $M1$ pro-chip-based workstation. In this configuration, a typical simulation takes between 5–7 seconds when the spatial domain is taken to be $[0, 1]$, and is partitioned into $1000$ sub-intervals. The simulations' parameters are chosen to meet the parametric constraints of the theorems. In the spatially explicit setting, we have used the comparison theory to determine point-wise restrictions on the fear functions such that competitive exclusion or strong competition-type dynamics abound. These are shown via Theorems 4.8, 4.10 and 4.15, whereas Figures 8, 9 and 12 can infer the numerical validation of this result. We also have numerically validated that in the spatially explicit setting, fear can change a situation of weak competition to a competitive exclusion type scenario for fear functions with certain $\mathbb{L}^{1}$ restrictions. This is proved in Theorem 4.13 and Lemma 4.14. Moreover, we have performed a time series analysis over a long period to numerically validate the result, see Figures 14–16. Various heterogeneous fear functions are constructed to demonstrate these results numerically; see Figure 17b. All the details about the parameters and their choice are explicitly mentioned in the captions of the figures of numerical simulation. Essentially in all of the PDE and ODE simulation, parametric values are chosen that meet the parametric restrictions of the Lemmas and Theorems that are being numerically demonstrated. Also we assume fairly competitive populations, and do not want one population exceeding more than three times the other. This (under scaling) restricts the choice of the parameters to lie in the range $(0, 3]$. This also assists in easier visualization of the dynamics.

6.   Discussion

In this section we discuss several aspects of the lemmas and theorems described in the current manuscript. We first briefly recap our key results.

● Sufficiently large fear can change a situation of competitive exclusion, to a strong competition type scenario, where there is bi-stability between boundary equilibrium. Dynamically, this occurs via a transcritical bifurcation. This is shown via Lemma 3.12, see Figure 6.

● Fear in a certain parametric regime can change a situation of competitive exclusion to bi-stability between boundary equilibrium and interior equilibrium, see Figure 4. Dynamically, this occurs via a saddle-node bifurcation. This is shown via Lemma 3.11, see Figure 5. This is in sharp contrast with classical competition theory, where bi-stability occurs only between boundary equilibria.

● Sufficiently large fear can change a situation of weak competition to a competitive exclusion type scenario. This is shown via Lemma 3.5, see Figure 2.

● Fear cannot qualitatively change a strong competition type scenario. This is shown via Lemma 3.8, see Figure 3. Also, fear cannot produce periodic orbits. This is demonstrated via Lemma 3.10 and Corollary 1.

● In the spatially explicit setting, comparison theory is used to determine point-wise restrictions on the fear functions such that competitive exclusion or strong competition type dynamics abounds. These are shown via Theorems 4.8, 4.10 and 4.15, see Figures 8, 9 and 12.

● In the spatially explicit setting, fear can change a situation of weak competition to a competitive exclusion type scenario, for fear functions with certain $\mathbb{L}^{1}$ restrictions. This is shown via Theorem 4.13 and Lemma 4.14, see Figures ^[14,15,16]. In particular, the fear functions need not lie uniformly above the critical fear levels derived in the ODE case via Lemma 3.5.

● Various heterogeneous fear functions are constructed to demonstrate these results numerically, see Figure 17b. Applications of these to ecological as well as socio-political settings are discussed in Section 6.

We see that in the competitive exclusion case, where $(u^{*}, 0)$ is globally stable, a small amount of fear in $u$, can bring about a strong competition type situation - via a transcritical bifurcation. Where after bifurcation there appears an interior equilibrium, which is a saddle, and one has initial condition dependent attraction to $(u^{*}, 0)$ or $(0, v^{*})$. See Figures 4 and 6]. However, a slightly larger level of fear, can bring down the separatrix drastically see Figure 17a. Furthermore, here (depending on parametric restrictions) only a very small quantity of fear can create two interior equilibriums, see Figure 6, where there are both weak and strong competition dynamics at play, where most initial conditions are attracted to an interior equilibrium, but certain initial conditions $(u_{0}, v_{0})$ (where $u_{0} \ll 1, v_{0} \gg 1$), are attracted to $(0, v^{*})$. It is interesting to think about this in an invasion setting. If $u$ were an invasive species, and $v$ a resident species, then in the absence of fear in $u$, it can invade and thus exclude $v$. However, if the resident $v$ can instill just a small amount of fear in the invader $u$, coexistence for most initial conditions is possible. If one plays devil's advocate, and switches the role of $u$ and $v$, ($u$ resident, $v$ invader), then we see in the absence of fear in the resident, it will exclude the invader, but if it can instill some fear, then coexistence is possible, as earlier alluded to. However, having too much fear in either case, yields only one interior equilibrium, which is a saddle and we are in a strong competition type setting, with initial condition dependent attraction to either $(u^{*}, 0)$ or $(0, v^{*})$. Thus if coexistence (for most initial conditions) is sought after, it is advantageous to induce a little fear - but not too much fear. Dynamically, this is seen because in the former case we have a saddle-node bifurcation occurring first, followed by a transcritical bifurcation, see Figure 6b, as opposed to the later case, where only a transcritical bifurcation occurs, see Figure 6a. Also, a 2 parameter bifurcation in the $k-c_{1}$ parameter space reveals a rich higher co-dimensional bifurcation structure. A Bogdanov Takens (BT) bifurcation, as well as a saddle node transcritical (SNTC) bifurcation is observed. Future work could entail further theoretical investigations of the SNTC ^[74], herein numerical investigations are also apt ^[75,76]. In this work we assume fear comes from the physical presence of a competitor. However even past presence such as by scent or other chemical cues can cause fear ^[77]. It is worthwhile considering modeling such a scenario in future work.

This phenomenon becomes even more interesting in sociological or political settings, particularly when we are in the spatially explicit scenario. Consider the case of law enforcement trying to control crime in high crime areas, where traditionally control efforts have failed, indicative of the competitive exclusion scenario, or a coexistence scenario, where the criminal groups have not been managed so coexist at high levels, and have not been eradicated. The question becomes, what is the optimal design of the fear function $k(x)$, that could now yield a co-existence scenario–or one that can yield a competitive exclusion scenario. We see via Theorem 4.13, that this function could be very small in some areas of the spatial domain, and large enough in others, so that it would change dynamics - this is also seen numerically, where "hotspot" type fear functions have been utilized to generate a competitive exclusion type scenario, see Figures 14–16. Similar ideas for control and policing activities have been explored in ^[78,79]. Essentially, the fear function can be modeled through various functional types, that also preserve spatial features of the underlying domain. For instance, consider competition between two warring drug cartels, where the weaker cartel has certain territorial strongholds. This situation can be modeled through oscillatory functions like $\sin^2(nx)$, where $n\in \mathbb{N}$. Moreover, for complex ecological models, the result of these theorems also holds, where we need to model the fear responses by a sequence of functions or in the form of piece-wise functions. We display some of these functions in Figure 17b. Many of the relevant theorems have been tested by choosing such functions (See Figures 8, 10, 12 and 13). However, we have not incorporated non-smooth fear function or non-smooth density dependence herein ^{[31,80,81,82]}, and this remains an apt direction for future work. In essence the fear model considered can be placed in a much broader context than ecological. Corollary 1 enables us to provide a mechanistic description of the fear function in much broader contexts. Also, in sociopolitical domains, the fear effect does not have to be thought of as consumptive vs non-consumptive. For example a new small business that is fearful of a larger established business might not invest in advertising or expanding clientele in established market, rather search for niche markets - which might hinder its growth. Interesting future work could consider competition models with a fear effect where the outcome of the competition whether contest or scramble was explicit.

In the weak competition case without fear, once there is a critical level of fear, in either $u$ or $v$, a competitive exclusion type phenomenon will occur, see Lemma 3.5. This occurs again via a transcritical bifurcation, see Theorem 3.12 and Figure 6. This again is interesting in the spatially explicit setting, as via Theorem 4.8, where the fear function need not be above $k_{c}$ so as to induce competitive exclusion. We have explored certain types of fear functions herein, see Figures 8, 14–16, but the effect of functions such as in Figure 17b, remain to be explored. Here one could look at the effect of fear in fragmented domains vs domains that are not, via deriving or enforcing conditions relating the fear function $k(x)$ to the resource function $m(x)$, see similar ideas explored in ^[83]. The results are symmetric even if $u$ is assumed to be fearful of $v$, and one considers the spatially explicit case, i.e., when a competitive exclusion situation, can change to a bi-stability type situation, where the bi-stability is between the $(u^{*}, 0)$ and a coexistence equilibrium. A further increase of the fear function causes a strong stability situation to occur. Similar dynamics are seen in changing a weak competition case to a competitive exclusion case. We are also motivated to studying the fear effect in competitive systems modeled differently - such as bringing in this effect in the competitive coefficients, rather than the birth rate. In addition, for competitive systems, one might intuitively think that the "dominant" population causes fear in the inferior population, which further weakens the inferior population. However, simulations show this is incorrect. Particularly if one considers fear in both competitors. The interplay of the heterogenous fear functions $k(x), f(x)$, can actually lead to a co-existence scenario, even when one competitor is dominant with and without the fear effect. This counterintuitive situation warrants much further investigation.

An interesting question is the dynamical consequence of Allee effects on the fear models considered. Note, Allee effects in competition systems have been considered ^[84,85], and results show the possibility of alternate stable states, and even tri-stability situations possible, where the occurrence of two interior saddle equilibrium, cause the phase to be partitioned such that initial data can be attracted to $(u^{*}, 0)$, $(0, v^{*})$ or (0, 0). However, whether a weak or strong Allee effect is considered, periodic orbits have not been shown to occur ^[84,85,86]. The possibility of two interior saddles does not preclude the existence of heteroclinic orbits, say via a saddle-saddle collision–although such dynamics have not been rigorously proved to the best of our knowledge. The interplay of the fear effect with an Allee effect could well change this situation. For example with a simple strong Allee effect the Dulac function used in Lemma 3.10 would not be sufficient to claim non-existence of limit cycle dynamics. Thus such investigations would make for interesting future studies. It would be interesting to explore applications of these results to the case of competing political parties, where a nationally weaker party with smaller voter bank, could have several local strongholds. What levels of fear need to be instilled in those strongholds so as to enforce competitive advantage, is an apt question. Herein estimates on the $\mathbb{L}^{1}(\Omega)$ norm of $k(x)$, would be interesting to derive. Currently Lemma 4.14, only gives lower estimates - but sharp upper estimates or even just upper estimates are unknown, and would make for interesting future work. The effect of fear on the strong competition setting is perhaps the least interesting dynamically. Herein, having fear in $u$ or $v$, only shifts the interior saddle equilibrium - but qualitatively the dynamics remain the same, as in a strong competition type scenario persists, no matter what level of fear ($f$ or $k$) is chosen. Also it would be of interest to rigorously prove Conjectures 1 and 2. Our current numerical evidence clearly motivates their validity. It would also be of interest to consider fractional derivative extensions of the fear model (3.1). Such problems have been much investigated recently ^{[87,88,89,90]}.

Also, in ^[8], it is found that various monotonically decreasing fear functions essentially yield the same dynamics. This has not been tested in the current work, in the case of competitive systems. Neither have we investigated rigorously the case of both competitors being fearful of each other with possible different fear functions, indicative of different LOF for each competitor. Here again, upper estimates on the fear functions would be useful and we could allude to methods and techniques explored in ^[91] and in ^[92]. Here, one could consider a spatially dependent growth function as well, such as in ^[92], and attempt to derive conditions relating the function describing the resources to the fear function. Furthermore, it would make for interesting future work if certain choices of (density) dependence or fear functions which lead to degenerate dynamics, can cause periodic orbits. In the current scenario Lemma 3.10 and Corollary 1, do not allow limit cycle dynamics. All in all, we hope these questions lead to future investigations of the fear effect in competitive systems, as a host of rich applications exist.

Acknowledgement

VS and RP acknowledge valuable summer support from the National Science Foundation via DMS 1715044.

Conflict of interest

The author declare that they have no conflict of interest.

Appendix

A.1.   The case of $u$ fearing $v$

Here we shall consider the case when the competitor $u$ is being fearful of $v$. Modeling as earlier yields,

A.2.   Existence

The system (A.1) can have one or two positive interior equilibria, but not three. This is established via the following lemma:

Lemma A.2.1. Consider the ODE system (A.1), with $f = 0$. If we are in the competitive exclusion, weak-competition, or strong-competition setting, then for any $f > 0$, under the parametric restrictions in (A.2), there do not exist three positive interior equilibria.

Proof. On substituting $u^* = \frac{1}{b_1}(\frac{a_1}{1+fv^*}-c_1 v^*)$ in the second nullcline equation, we get the following third order polynomial in $v^*$:

where

We will require the following conditions to satisfy Descartes's rule of signs, so as to obtain three positive roots to the cubic equation. These are:

Let us prove this is an impossible claim by contradiction: First, assume the parameter set satisfies both competitive exclusion-state parametric restriction and (A.2). We know that competitive exclusion-state is asymptotically stable if $\frac{a_1}{a_2} > \max \Big\{ \frac{b_1}{c_2}, \frac{c_1}{b_2} \Big \}$. On using these parametric restriction, $-a_1c_2 + a_2b_1 < 0$, which is a contradiction to the last inequality in (A.2).

If the parameter set satisfies both the strong competition state parametric restriction and (A.2), then we have a contradiction because of the first inequality in (A.2), as under the strong competition state parametric restriction that inequality should be positive.

For the weak competition, recall the parametric restrictions:

Let us re-write the third inequality in (A.2),

On adding the second and third inequality in (A.2), we have

which is a contradiction as the added inequality should be negative by the parametric restrictions of weak-competition. This proves the lemma.

Also note:

Remark 12. The cubic equation obtained in (A.1) can be transformed to a quadratic polynomial by reducing $(b_1 f v^*+b_1)$ on both sides, i.e.,

where

Thus in essence models (3.1) and (A.1) are symmetric, and it suffices to consider (3.1).