Population Dynamics with Fear-Induced Allee Effects: A
Mathematical Study on Prey–Predator Interaction
Rajiba Lochan Sahoo, Roll No: 24MAT207
Department of Mathematics
Institute of Chemical Technology, Mumba
Under the guidance of :
Dr. Debasmita Mukherjee
November 15 ,2025
1. Introduction
Predators not only kill prey but also create fear that changes prey behavior. Because of fear, prey animals eat
less, move less, and focus more on survival than reproduction. These anti-predator behaviors help prey avoid
being hunted but reduce their growth and reproduction rates.
This supports the Allee effect, which means that when the population is small or scared, it grows more slowly.
Earlier studies focused only on direct killing by predators. But fear itself can change the population balance
— it can lower birth rates and increase the chance of extinction. So, adding fear into mathematical models helps
us understand real-life predator–prey interactions more accurately.
2. Single-Species Model with Generalized Predator
In this section, we study how fear affects the growth of a single prey population.
Normally, the population grows following the rule:
where:
•
•
•
•
To include the fear effect (caused by predators), the model adds a fear function that reduces reproduction
1
Here:
• = predator population,
• f = level of fear (how strongly prey react to predators).
The term means that as fear or predator density increases, the growth rate of prey decreases.
The factor ,satisfies the following conditions:
(a) F(0,v) = 1: if there is no anti-predator behaviors of prey (no fear of predation) then no reduction in the prey
produc- tion.
(b) F(f,0) = 1: if there is no predator, then no reduction in prey production due to anti-predator behaviors.
(c) If anti-predator behavior is very large, then prey production reduces to zero.
(d) if predator population is very large then prey reproduction reduced to zero, due to large
anti- predator behaviors.
%% Figure 1: Population Density vs Per Capita Growth Rate (Effect of Fear)
clear; close all; clc;
% Parameters
theta = 0.2; % Allee threshold
v = 1; % predator density (constant)
f_vals = [0, 0.1, 0.5, 1]; % different fear levels
N = linspace(0, 1, 500); % prey density range (normalized)
% Function for per capita growth rate (from Model 2.1)
pgr = @(N, theta, f, v) N .* (1 - N) .* (N - theta) ./ (1 + f * v);
% Plot
figure; hold on; box on;
for i = 1:length(f_vals)
plot(N, pgr(N, theta, f_vals(i), v), 'LineWidth', 1.8);
end
% Axis labels and styling
xlabel('Prey Density (N)', 'FontName', 'Times New Roman', 'FontSize', 12);
ylabel('Per Capita Growth Rate (PGR)', 'FontName', 'Times New Roman', 'FontSize',
12);
title('Fig. 1: Population Growth Rate vs Density for Different Fear Levels', ...
'FontName', 'Times New Roman', 'FontSize', 12);
legend({'Single Allee (f = 0)','Allee with fear (f = 0.1)','Allee with fear (f =
0.5)','Allee with fear (f = 1)'}, 'Location','northwest');
2
grid on;
yline(0,'k-');
xlim([0 1]);
ylim([-0.05 0.25]);
Fig. 1. Population density versus pgr(per-capital growth Rate) plot for the Model (2.2) . We fixed the generalized
predator density as , and the other parameter is fixed as and we vary the fear effect parameter f .
The red curves are indicating the combined effect of fear Allee pgr function ( f = 0 . 1 (solid red curve), f = 0 .
5 (yellow curve) and f = 1 (violet curve)), and the blue curve is the single Allee pgr function (for f = 0 ). We see
that, the shape of the pgr curve is greatly influenced by the fear effect parameter f . (For interpretation of the
references to color in this figure legend, the reader is referred to the web version of this article.).
(e) : as the anti-predator behaviors increases, the production of prey decreases.
(f) : as the predator population size increases, the production of prey decreases.
Moreover, the Model (2.1)can be written as
3
which justify the fact that such modeling approach is also true, when we assume fear significantly reduces the
reproduction process only (does not affect the death). Further, we non-dimensionalize the above model (2.1)
and get the single species model as
where .The system (2.2) has three equilibria,namely,
Among these , is unstable, while and are stable. If the initial population , the prey
goes extinct. Otherwise, it grows and reaches its carrying capacity.
3. Prey-predator (two species) model with cost of fear
In this section, we study a prey-predator model with Allee effect in prey and the reduced reproduction of prey
due to cost of fear. For simplicity, we consider the linear functional response. Thus our two-species model
becomes,
where:
• prey density, :predator density
• prey intrinsic rate, k: carrying capacity, (with 0 ):Allee threshold
• a:predation rate , :conversion efficiency of prey into predator growth
• m:predator natural mortality
• fear level(larger stronger reduction of prey reproduction)
The factor multiplies the prey reproduction term and reduses it when either fear for predator density v
increases.
Here, the parameter a is the predation rate, is the conversion efficiency of predator by consuming prey and
m is the predator’s natural death rate. First, to non-dimensionalized the model, consider and ,
then the system (3.1) reduces to
4
In two term first term growth rate of prey .This growth is affected by two things — the Allee effect and fear
of [Link] fP value is high Decreses growth Rate. Second term So, this term represents the loss
of prey population due to predation.
In The first term (N×P) means the predator grows by eating prey(term NP), Second term (−mP)
represents the natural death of predators.
Where , is the time scale separation between two populations, .By the
following this approach, the number of parameters reduced from seven to four. Moreover, we can separate the
time scale between populations, as in nature, different species have different life span. =1, represent no time
scale difference between pop- ulations and smaller value of ( ) indicates a large time scale separation
between two populations life cycle.
By using the standard basic theory of ordinary differential equation systems, we can show that for every
initial condition , the system(3.2)has a unique solution, and the solution remains positive and
bounded.
The Model (3.2) has three boundary equilibria, namely,
which always exist. It has a unique interior equilibrium
which exists if .Below we discuss the stability conditions of all equilibria for the model (3.2).
1. The extinction equilibrium is always LAS, as the eigenvalues corresponding to this equilibrium are
2. The axial equilibrium is a source if and a saddle if , as the eigenvalues associated with
this equilibrium are and
5
3. The other axial equilibrium is a saddle if and is LAS (locally asymptotically
stable) if as the eigenvalues corresponding to this equilibrium are: and
4. The Jacobian matrix at the interior equilibrium of the model (3.2) is given by
The corresponding characteristic equation is
where
%Fig.2 :Phase Portraits for predator extinction and predation-driven
%extinction
% Parameters
eps = 1.0; % epsilon
theta = 0.2; % Allee threshold
f = 1.0; % fear
m_vals = [1.1, 0.15]; % (a) predator extinction, (b) predation-driven extinction
% RHS as anonymous functions
dN = @(N,P,eps,theta,f,m) (1/eps) * ( N.*(1-N).*(N-theta)./(1 + f.*P) - N.*P );
dP = @(N,P,m) N.*P - m.*P;
% Grid for vector field
Ngrid = linspace(0,1,30);
Pgrid = linspace(0,0.65,30);
[Nm,Pm] = meshgrid(Ngrid,Pgrid);
figure('Position',[100 100 1100 430]);
for j = 1:2
m = m_vals(j);
subplot(1,2,j); hold on; box on;
6
% Vector field (normalized for nice arrows)
U = dN(Nm,Pm,eps,theta,f,m);
V = dP(Nm,Pm,m);
speed = sqrt(U.^2 + V.^2);
speed(speed==0) = 1;
quiver(Nm, Pm, U./speed, V./speed, 0.6, 'k'); % 0.6 scale for arrow length
% Prey nullcline (explicit positive root of f P^2 + P - (1-N)(N-theta) = 0)
Nline = linspace(0.001,0.999,800);
if f == 0
% when f==0 the equation reduces to P = (1-N)*(N-theta) / 1 (linear)
P_nc = max(0, (1 - Nline).*(Nline - theta));
else
disc = 1 + 4*f.*(1 - Nline).*(Nline - theta);
disc(disc < 0) = 0; % avoid complex due to negative
roundoff
P_nc = (-1 + sqrt(disc))./(2*f); % positive branch
P_nc(P_nc < 0) = 0; % keep nonnegative
end
plot(Nline, P_nc, 'm--', 'LineWidth', 1.6); % prey nullcline (dN/dt=0)
% Predator nullclines
plot([0 1], [0 0], 'k-', 'LineWidth', 1.0); % P = 0
xline(m, '--b', 'LineWidth', 1.2); % N = m (predator nullcline)
% Sample initial conditions and trajectories
initials = [0.05 0.02; 0.3 0.02; 0.6 0.1; 0.2 0.2; 0.8 0.05; 0.4 0.25];
tspan = [0 300];
for k=1:size(initials,1)
y0 = initials(k,:);
% integrate
[t, sol] = ode45(@(t,y) [ dN(y(1), y(2), eps, theta, f, m); dP(y(1), y(2),
m) ], tspan, y0);
plot(sol(:,1), sol(:,2), 'LineWidth', 1.2);
plot(sol(1,1), sol(1,2), 'ko', 'MarkerFaceColor','k', 'MarkerSize',4); %
start point
end
% Plot equilibria E0, Eθ, E1
plot(0,0,'ks','MarkerFaceColor','c','MarkerSize',8); text(0.02, 0.02,
'E_0', 'FontSize',10);
plot(theta,0,'ks','MarkerFaceColor','g','MarkerSize',8);
text(theta+0.015, 0.02, 'E_\theta', 'FontSize',10);
plot(1,0,'ks','MarkerFaceColor','r','MarkerSize',8); text(0.92, 0.02,
'E_1', 'FontSize',10);
% Interior equilibrium E* if exists (θ < m < 1)
if (theta < m) && (m < 1)
A = (1-m)*(m-theta);
disc = 1 + 4*f*A;
7
Pstar = (-1 + sqrt(disc)) / (2*f);
if isreal(Pstar) && Pstar > 0
plot(m, Pstar, 'mo', 'MarkerFaceColor','m', 'MarkerSize',8);
text(m+0.01, Pstar+0.02, 'E_*', 'FontSize',10, 'Color','m');
end
end
% Formatting
xlim([0 1]); ylim([0 0.6]);
xlabel('Prey (N)', 'FontSize',11); ylabel('Predator (P)', 'FontSize',11);
if j==1
title(sprintf('(a) Predator extinction for m = %.2f', m));
else
title(sprintf('(b) Predation-driven extinction for m = %.2f', m));
end
grid on;
set(gca, 'FontSize', 10);
end
% Legend (manual)
subplot(1,2,1);
legend({'Vector field','Prey nullcline','P=0','N=m (predator
nullcline)','Trajectories','start','E_0','E_\theta','E_1','E_* (if present)'}, ...
'Location','northeastoutside','FontSize',8);
[Link] the two extreme outcomes: predator extinction for m =1.1 and predation-driven extinction for
m=0.15 .we fix the additional parameter vakues as =1, =0.2, .
8
3.1. Effect of time scale separation among two populations:
% Fig.3: Phase portraits for different m (eps=1, theta=0.2, f=1)
eps = 1.0; theta = 0.2; f = 1.0;
m_list = [0.65, 0.57, 0.5345, 0.50]; % (a)(b)(c)(d)
% RHS functions
dN = @(N,P,eps,theta,f,m) (1/eps) * ( N.*(1-N).*(N-theta)./(1+f.*P) - N.*P );
dP = @(N,P,eps,theta,f,m) N.*P - m.*P;
% prey nullcline: positive root of f P^2 + P - (1-N)(N-theta) = 0
prey_nc = @(N,f,theta) (-1 + sqrt(max(0, 1 + 4*f.*(1-N).*(N-theta)))) ./ (2*f);
Ngrid = linspace(0,1,26); Pgrid = linspace(0,0.45,26);
[Nm,Pm] = meshgrid(Ngrid,Pgrid);
figure('Position',[80 80 1100 800]);
for j = 1:4
m = m_list(j);
u = dN(Nm,Pm,eps,theta,f,m);
v = dP(Nm,Pm,eps,theta,f,m);
subplot(2,2,j); hold on; box on; axis square;
% vector field
q = quiver(Nm,Pm,u,v,'AutoScale','on','AutoScaleFactor',1.1,'Color',[0 0 0]);
[Link] = 0.6;
% nullclines
Nline = linspace(0,1,600);
Pnc = prey_nc(Nline,f,theta);
plot(Nline, Pnc, 'm--','LineWidth',1.6); % prey nullcline (dN/dt=0)
yline(0,'k-','LineWidth',1.0); % P=0 (predator axis)
xline(m,'Color',[0.9 0.5 0],'LineWidth',1.5); % N=m (orange-ish)
9
% sample trajectories (initial points)
inits = [0.08 0.04; 0.2 0.02; 0.35 0.1; 0.55 0.12; 0.75 0.03; 0.45 0.25];
tspan = [0 600];
colors = lines(size(inits,1));
for k=1:size(inits,1)
[t,sol] = ode45(@(t,y)[ dN(y(1), y(2), eps, theta, f, m); dP(y(1), y(2),
eps, theta, f, m) ], tspan, inits(k,:));
plot(sol(:,1),sol(:,2),'b-','LineWidth',1.6); % trajectory
plot(sol(1,1),sol(1,2),'ko','MarkerFaceColor','k'); % start point
end
% E0 = (0,0)
plot(0,0,'ks','MarkerFaceColor','k','MarkerSize',8);
text(0+0.02, 0+0.01, 'E_0','FontWeight','bold');
% E_theta = (theta,0)
plot(theta,0,'ks','MarkerFaceColor','k','MarkerSize',8);
text(theta+0.02, 0+0.01, 'E_\theta','FontWeight','bold');
% E1 = (1,0)
plot(1,0,'ks','MarkerFaceColor','k','MarkerSize',8);
text(1-0.18, 0+0.01, 'E_1','FontWeight','bold');
% interior E* if exists (theta < m < 1)
if (theta < m) && (m < 1)
disc = 1 + 4*f*(1-m)*(m-theta);
if disc >= 0
Pstar = (-1 + sqrt(disc)) / (2*f);
if Pstar > 0
plot(m,Pstar,'ro','MarkerFaceColor','r','MarkerSize',7);
text(m+0.02, Pstar+0.01, sprintf('E_* (%.3g,%.3g)', m,
Pstar),'FontWeight','bold','Color','r');
end
end
end
% formatting
xlim([0 1]); ylim([0 0.45]); xlabel('Prey (N)'); ylabel('Predator (P)'); grid
on;
titles = {sprintf('(a) m=%.4g',m_list(1)), sprintf('(b) m=%.4g',m_list(2)), ...
sprintf('(c) m=%.4g',m_list(3)), sprintf('(d) m=%.4g',m_list(4))};
title(titles{j}, 'FontSize',11);
end
subplot(2,2,1);
h_q = plot(nan,nan,'k-'); hold on;
h_nc = plot(nan,nan,'m--','LineWidth',1.6);
h_tr = plot(nan,nan,'b-','LineWidth',1.6);
h_e = plot(nan,nan,'ro','MarkerFaceColor','r');
10
Fig. 3 (reproduced). Phase planes for , , with different m. (a) : stable interior;
(b) : stable limit cycle; (c) : heteroclinic-type loop; (d) : unstable interior. Blue:
trajectories; magenta dashed: = 0;orange line: ; black line:
%% Fig. 4: Stable interior for different epsilon (eps = 1, 0.25, 0.1)
% Parameters as in the figure/caption
theta = 0.2; f = 1.0; m = 0.64; % m>(1+theta)/2=0.6 -> interior is
stable
eps_list = [1.0, 0.25, 0.1];
IC = [0.45; 0.06]; % initial state (N0, P0)
tspan = [0 200];
11
% Convenience: interior equilibrium
Pstar = (-1 + sqrt(1 + 4*f*(1-m)*(m-theta)))/(2*f); % coexistence predator level
figure('Position',[80 80 1050 320]);
for k = 1:numel(eps_list)
eps = eps_list(k);
[t,sol] = ode45(@(t,y) rhs_pp(y,eps,theta,f,m), tspan, IC);
subplot(1,3,k); hold on; box on; grid on;
plot(t, sol(:,1), 'LineWidth', 1.6); % N(t)
plot(t, sol(:,2), 'LineWidth', 1.3); % P(t)
yline(m,'--'); % N* = m
yline(Pstar,':'); % P* = Pstar
xlabel('Time'); ylabel('Population density');
title(sprintf('Stable interior for \\epsilon = %.2f (E^*=(%.2f, %.3f))', eps,
m, Pstar));
legend({'Prey N','Predator P','N^*','P^*'}, 'Location','northeast');
end
% local function (system RHS)
function dY = rhs_pp(Y,eps,theta,f,m)
N = Y(1); P = Y(2);
dN = (1/eps)*( N*(1-N)*(N-theta)/(1+f*P) - N*P );
dP = N*P - m*P;
dY = [dN; dP];
end
Fig. 4 (reproduced). Time series for with {1, 0.25, 0.1}.The interior coexistence is
stable in all cases; smaller (faster prey) gives faster convergence.
% Stability region in (theta,m) plane (reproducing Fig.4 style)
% White = E0 only (extinction equilibrium stable)
12
% Green = E0 and E* bistable (interior equilibrium exists and is LAS)
% Yellow = E0 and E1 bistable (axial E1 also stable)
% parameter grid
theta_vec = linspace(0,0.5,201); % horizontal axis (theta)
m_vec = linspace(0.5,1.5,151); % vertical axis (m)
[TH, M] = meshgrid(theta_vec, m_vec);
% Conditions from paper (non-dimensional model)
% Existence of interior equilibrium E*: theta < m < 1
existE = (TH < M) & (M < 1);
% Stability of interior equilibrium E* (LAS condition used in paper)
% Paper shows A>0 if m < (1+theta)/2 and A<0 if m > (1+theta)/2.
% But final LAS condition for E* in paper: max{theta, (1+theta)/2} < m < 1.
stableE = ( M > max(TH, (1+TH)/2) ) & (M < 1);
% Stability of E1 (axial equilibrium at N=1): E1 is LAS if m > 1
stableE1 = (M > 1);
% Regions
region_E0_only = ~existE & ~stableE1; % no interior, E1 not stable -> only
E0 stable
% however E0 is always LAS in their analysis; we mark white where other equilibria
not stable
region_E0_only = ~(stableE | stableE1); % simpler: neither E* nor E1 stable
region_E0_Estar = stableE; % green: interior exists & LAS (with
E0 also LAS)
region_E0_E1 = stableE1; % yellow: E1 LAS (with E0 also LAS)
% Plot setup
figure('Color','w','Position',[100 100 900 650]);
hold on; box on;
% Plot green region (E0 & E*) as dense dots
idxG = find(region_E0_Estar);
scatter(TH(idxG), M(idxG), 18, [0 0.5 0], 'filled', 'MarkerEdgeColor','none'); %
dark green
% Plot yellow region (E0 & E1)
idxY = find(region_E0_E1);
scatter(TH(idxY), M(idxY), 18, [0.95 0.85 0.1], 'filled',
'MarkerEdgeColor','none'); % yellow
% Plot white region by plotting faint grey dots (so empty-looking)
idxW = find(region_E0_only);
scatter(TH(idxW), M(idxW), 8, [0.85 0.85 0.85], '.'); % light grey / appears white
% Add horizontal line at m = 1 and vertical gridlines to match paper look
13
yline(1,'k-','LineWidth',1);
for t = 0:0.05:0.5
xline(t,'Color',[0.9 0.9 0.9],'LineStyle','-');
end
% Annotate typical labels as in paper
text(0.18, 1.25, 'E_0 \cup
E_1','FontSize',20,'FontWeight','bold','HorizontalAlignment','center');
text(0.2, 0.85, 'E_0 \cup E_*','FontSize',20,'FontWeight','bold','Color',[0 0.45
0],'HorizontalAlignment','center');
text(0.38, 0.55, 'E_0','FontSize',20,'HorizontalAlignment','center');
% Axes and title
xlim([0 0.5]); ylim([0.45 1.5]);
xlabel('\theta','FontSize',16); ylabel('m','FontSize',16);
title('Stability region in (theta, m) ( epsilon = 1,f =
1 )','FontSize',14,'FontWeight','normal');
set(gca,'FontSize',14,'LineWidth',1,'Box','on');
set(gca,'YTick',0.5:0.1:1.5);
axis square;
% Legend (manual small patch markers)
h1 = scatter(-1,-1,80,[0 0.5 0],'filled'); %# dummy for legend
h2 = scatter(-1,-1,80,[0.95 0.85 0.1],'filled');
h3 = scatter(-1,-1,8,[0.85 0.85 0.85],'.');
legend([h1 h2 h3],{'E_0 \cup E_*','E_0 \cup E_1','E_0
only'},'Location','southoutside','Orientation','horizontal');
legend boxoff;
% polish
set(gcf,'Renderer','painters');
hold off;
14
%% Fig. 5: Predator density at coexistence vs fear f (ε=1, θ=0.2, m=0.7)
% parameters (must satisfy θ < m < 1 for coexistence) ---
theta = 0.2;
m = 0.7; % 0.2 < 0.7 < 1 -> OK
epsi = 1; % not used in formula here, kept for caption consistency
% fear grid
fvec = linspace(0,5,400);
% (coexistence constant)
A = (1 - m) * (m - theta); % > 0 since θ < m < 1
% P*(f): closed-form with careful f=0 limit and robust sqrt ---
15
Pstar = zeros(size(fvec));
Pstar(1) = A; % limit as f -> 0 : P* -> A
disc = max(0, 1 + 4 * fvec(2:end) .* A); % protect against tiny negatives
Pstar(2:end) = (-1 + sqrt(disc)) ./ (2 * fvec(2:end));
%
figure('Position',[100 420 900 380],'Color','w');
plot(fvec, Pstar, 'LineWidth', 2); hold on; grid on; box on;
xlabel('fear level, f','FontName','Times New Roman');
ylabel('Predator density at coexistence, P^*(f)','FontName','Times New Roman');
title(sprintf('Fig. 6: P^*(f) for \\epsilon = %g, \\theta = %.1f, m = %.2f', epsi,
theta, m), ...
'FontName','Times New Roman','FontWeight','bold');
xlim([0 5]);
yl = [max(0, min(Pstar)-0.01), max(Pstar)+0.01];
ylim(yl);
% annotate stability text (optional)
text(0.15, Pstar(1), ' stable','FontName','Times New Roman');
% dashed baseline at zero (like paper style)
yline(0,'k--','LineWidth',1);
Fig.5. predator density at the unique interior equlibrium versus fear level f for = 1, The
curve decreases monotonically, showing that higher fear reduces the predator density at coexistence.
16
[Link]
From the analysis of the model, we conclude that the stability of the coexistence equilibrium does not depend
on the cost of fear parameter f. However, as the fear level increases, the predator population at equilibrium
decreases.
This happens because higher fear reduces prey reproduction, leading to less food for predators. Therefore,
the predator density gradually decreases with increasing fear, as shown in Fig. 6. The system remains stable
even when fear increases.
[Link]
[Link], X., Zanette, L., & Zou, X. (2017). Modelling the fear effect in predator–prey interactions. Journal of
Mathematical Biology, 75(5), 1231–1251.
[Link], S. K. (2018). Population dynamics with multiple Allee effects induced by fear factors – A
mathematical study on prey–predator interactions. Applied Mathematical Modelling, 64, 1–14. [Link]
10.1016/[Link].2018.07.021
[Link], J. D. (2002). Mathematical Biology I: An Introduction (3rd ed.). Springer.
[Link], N., & Adak, D. (2019). Mathematical Biology I: An Introduction. Springer.
17