Supplementary notebook 

Heterogeneity in behaviour and movement can influence the stability of predator-prey periodic travelling waves - Renato Andrade*¹ and Christina A. Cobbold²

Bulletin of Mathematical Biology

¹ School of Mathematics and Statistics,  University of Glasgow, Glasgow, G12 8QQ, United Kingdom

² Boyd Orr Centre for Population and Ecosystem Health, University of Glasgow, Glasgow, G12 8QW, United Kingdom

*Corresponding author: Renato.AntunesCostadeAndrade@glasgow.ac.uk

(*

Computation and evaluation of ω1 according to the procedure described by Sherrat et al. (2003) The Effects of the Size and Shape of Landscape Features on the Formation of Traveling Waves in Cyclic Populations. 

\partial_T h= \partial^2_X h+(1-h)h- \Big[\frac{C h}{1+C h}+ \frac{\beta C\eta h}{1+C\eta h}\Big]p
    \partial_T p= \delta\partial^2_X p- \frac{p}{E_1B_1}+ \Big[\frac{C h}{B_1(1+C h)}+ \frac{\beta C\eta h}{B_1(1+C\eta h)}\Big]p

*)

Section 0: Define the reaction functions and compute the steady states

ClearAll[β,c,η,B1,B2,E1,F,G,p,h,w, fullsol, hs, ps]
$Assumptions=η >0 && c >0 &&B1 > 0 &&B2 > 0 && A > 0 &&β>0;

F[h_,p_]:= c p h/(B1(1+ c h))+c β η p h/(B1(1+c η h)) - p/(E1 B1) 
G[h_,p_]:= - c p h/(1+ c h)-β c η p h/(1+ c η h)  + h(1 - h)
w[h_]:= (1-h)/(c/(1+c h) + β c η/(1 + c η h))
fullsol = FullSimplify[Solve[F[h,p] == 0 && G[h,p] == 0, {h,p}]];
(*coexistence steady state*)
hs = fullsol[[3,1,2]]
ps = fullsol[[3,2,2]]

Section 1: Compute the jacobian and C_crit

ClearAll[β,J,c,η,B1,B2,E1,μ,trace,det,r,i,Λ0,Ω0,Ω1,p,h, ccrit1,ccrit2,ξ1,ξ2]

J ={{D[F[h,p],p],D[F[h,p],h]},{D[G[h,p],p],D[G[h,p],h]}};
trace = Tr[J] /.  {h -> hs,p-> ps};
det = Det[J]/.  {h -> hs,p-> ps};
ccrit1 = Solve[trace==0,c][[2,1,2]];

(*equivalent way to solve Tr[J] = 0, but more efficient - See stability analysis in Cobbold et al (2022) for details*)
ccrit2 =FullSimplify[Solve[w'[hs]==0,c][[1,1,2]]]
Simplify[ccrit2/. η -> 1] 
FortranForm[ccrit2 /. {η-> eta, β -> beta}];

ClearAll[β,J,c,η,B1,B2,E1,r,i,Λ0,Ω0,Ω1,p,h, ccrit]
ccrit = ccrit2;

(*Relevant parameters for the normal form*)
r = trace/2;
i =Sqrt[det-trace*trace/4];
Λ0 = D[r,c] /. c -> ccrit;
Ω0 = i /. c -> ccrit;
Ω1 = D[i,c] /. c -> ccrit;

Set the fixed parameters

ClearAll[β,E1,B1,η ]
(*Fixed paramaters*)
"β = "
β =1.0
"η = "
η=7.0

(* Tests and constraints for existence of the Hopf-bifurcation - See Cobbold et al. (2022)*)
"E1 must be bigger than"
1/(1 + β)

Assuming[E1>0,Simplify[ccrit]]
"must be bigger than"
(β η + 1)/(β η^2 +1)

Section 2: Normal form transformation

ClearAll[c,CurlyF,CurlyG,p,h,S,μ,J]

J ={{D[F[h,p],p],D[F[h,p],h]},{D[G[h,p],p],D[G[h,p],h]}} /.{h -> hs,p-> ps};

(*Construct the transformation matrix*)
S = {{i/J[[2,1]],-(J[[2,2]]-r)/J[[2,1]]},{0,1}};
c=ccrit + μ;
$Assumptions = μ>0;
(*To obtain the non-linear part of the transformed equations, subtract the linear part with μ = 0*)
{p,h} = S.{P,H}+{ps,hs};
(*normal form non-linear part*)
{CurlyF,CurlyG} = Inverse[S].{F[h,p],G[h,p]} -{-Ω0 H,Ω0 P} /. μ -> 0;

(*testing - all first derivatives should be 0 at (H,P) = (0,0). Use random values for E1,B1 because it takes large amounts of time to calculate in the general model*)

Chop[{CurlyF,CurlyG} /. {H-> 0,P -> 0, E1 -> RandomReal[], B1 -> RandomReal[]}]
Chop[{D[CurlyF,H],D[CurlyF,P]}  /. {H-> 0,P -> 0, E1 -> RandomReal[], B1 -> RandomReal[]}]
Chop[{D[CurlyG,H],D[CurlyG,P]}  /. {H-> 0,P -> 0, E1 -> RandomReal[], B1 -> RandomReal[]}]

Section 4: Evaluate ω_1 and generate the values for the E_1xB_1 plot (Figure 5)

ClearAll[E1,B1,listω1,stabilitytest,Nx,Ny,Maxy,Maxx,Miny,Minx,p,h,ω1,HS,PS,Listccrit,listomega0]
(*grid size*)
Nx = 20;
Ny =20;

(*parameter range*)
Maxy = 10;
Miny =1.00;
Maxx = 5;
Minx = 0.55;
μ = 0;
"ccrit,hs,ps ="

AbsoluteTiming[

(*testing if hs, ps and ccrit are within the correct range*)
HS=Table[If[hs < 0 || hs > 1,-1, 1],{B1,Miny,Maxy,(Maxy-Miny)/(Ny-1)},{E1,Minx,Maxx,(Maxx-Minx)/(Nx-1)}];
PS=Table[If[ps < 0 ,-1, 1],{B1,Miny,Maxy,(Maxy-Miny)/(Ny-1)},{E1,Minx,Maxx,(Maxx-Minx)/(Nx-1)}];
Listccrit = Table[If[ccrit < (β η + 1)/(β η^2 +1) ,-1, 1],{B1,Miny,Maxy,(Maxy-Miny)/(Ny-1)},{E1,Minx,Maxx,(Maxx-Minx)/(Nx-1)}];
ArrayPlot[HS,AxesLabel->Automatic,FrameTicks->{{Miny,Maxy},{Minx,Maxx}},FrameLabel->{B1,E1},ColorRules->{1 -> White,-1-> Black},DataReversed->True]
ArrayPlot[PS,AxesLabel->Automatic,FrameTicks->{{Miny,Maxy},{Minx,Maxx}},FrameLabel->{B1,E1},ColorRules->{1 -> White,-1-> Black},DataReversed->True]
ArrayPlot[Listccrit,AxesLabel->Automatic,FrameTicks->{{Miny,Maxy},{Minx,Maxx}},FrameLabel->{B1,E1},ColorRules->{1 -> White,-1-> Black},DataReversed->True]

(*Computing Λ1 and Ω2 calculating each derivative separately for efficiency*)
ClearAll[dfppp,dfphh,dgpph,dfpp,dgpp,dfhh,dghh,dfph,dgph,dgppp,dfpph,dfhhh,listΩ0,Lambda1,Omega2,dummy]
dummy = D[CurlyF,{P,3}] /. {H-> 0,P -> 0};
dfppp = Table[dummy,{B1,Miny,Maxy,(Maxy-Miny)/(Ny-1)},{E1,Minx,Maxx,(Maxx-Minx)/(Nx-1)}];
dummy = D[CurlyF,P,H,H] /. {H-> 0,P -> 0};
dfphh = Table[dummy,{B1,Miny,Maxy,(Maxy-Miny)/(Ny-1)},{E1,Minx,Maxx,(Maxx-Minx)/(Nx-1)}];
dummy =D[CurlyG,P,P,H]/. {H-> 0,P -> 0};
dgpph = Table[dummy,{B1,Miny,Maxy,(Maxy-Miny)/(Ny-1)},{E1,Minx,Maxx,(Maxx-Minx)/(Nx-1)}];
dummy =D[CurlyG,{H,3}]/. {H-> 0,P -> 0};
dghhh = Table[dummy,{B1,Miny,Maxy,(Maxy-Miny)/(Ny-1)},{E1,Minx,Maxx,(Maxx-Minx)/(Nx-1)}];
dummy =D[CurlyF,{P,2}]/. {H-> 0,P -> 0};
dfpp = Table[dummy,{B1,Miny,Maxy,(Maxy-Miny)/(Ny-1)},{E1,Minx,Maxx,(Maxx-Minx)/(Nx-1)}];
dummy =D[CurlyG,{P,2}]/. {H-> 0,P -> 0};
dgpp = Table[dummy,{B1,Miny,Maxy,(Maxy-Miny)/(Ny-1)},{E1,Minx,Maxx,(Maxx-Minx)/(Nx-1)}];
dummy =D[CurlyF,{H,2}]/. {H-> 0,P -> 0};
dfhh = Table[dummy,{B1,Miny,Maxy,(Maxy-Miny)/(Ny-1)},{E1,Minx,Maxx,(Maxx-Minx)/(Nx-1)}];
dummy =D[CurlyG,{H,2}]/. {H-> 0,P -> 0};
dghh = Table[dummy,{B1,Miny,Maxy,(Maxy-Miny)/(Ny-1)},{E1,Minx,Maxx,(Maxx-Minx)/(Nx-1)}];
dummy =D[CurlyF,P,H]/. {H-> 0,P -> 0};
dfph = Table[dummy,{B1,Miny,Maxy,(Maxy-Miny)/(Ny-1)},{E1,Minx,Maxx,(Maxx-Minx)/(Nx-1)}];
dummy =D[CurlyG,P,H]/. {H-> 0,P -> 0};
dgph =Table[dummy,{B1,Miny,Maxy,(Maxy-Miny)/(Ny-1)},{E1,Minx,Maxx,(Maxx-Minx)/(Nx-1)}];
dummy= D[CurlyG,{P,3}] /. {H-> 0,P -> 0};
dgppp = Table[dummy,{B1,Miny,Maxy,(Maxy-Miny)/(Ny-1)},{E1,Minx,Maxx,(Maxx-Minx)/(Nx-1)}];
dummy =D[CurlyG,P,H,H]/. {H-> 0,P -> 0};
dgphh = Table[dummy,{B1,Miny,Maxy,(Maxy-Miny)/(Ny-1)},{E1,Minx,Maxx,(Maxx-Minx)/(Nx-1)}];
dummy =D[CurlyF,P,P,H]/. {H-> 0,P -> 0};
dfpph = Table[dummy,{B1,Miny,Maxy,(Maxy-Miny)/(Ny-1)},{E1,Minx,Maxx,(Maxx-Minx)/(Nx-1)}];
dummy =D[CurlyF,{H,3}]/. {H-> 0,P -> 0};
dfhhh = Table[dummy,{B1,Miny,Maxy,(Maxy-Miny)/(Ny-1)},{E1,Minx,Maxx,(Maxx-Minx)/(Nx-1)}];

dummy = Ω0/. {H-> 0,P -> 0};
listΩ0 = Table[dummy,{B1,Miny,Maxy,(Maxy-Miny)/(Ny-1)},{E1,Minx,Maxx,(Maxx-Minx)/(Nx-1)}]; 



Lambda1 = -(dfppp+dfphh+dgpph + dghhh)/16+(dfpp dgpp - dfhh dghh - dfph(dfpp + dfhh)+ dgph(dghh + dgpp))/(16 listΩ0);
Omega2 = (dgppp + dgphh - dfpph - dfhhh)/16 +(dfpp(dgph - dfhh) + dghh(dfph - dgpp) - 3dfhh(dgph+dfhh)-3dgpp(dfph+dgpp) - 2(dfpp+dfhh - dgph)^2-2(dgpp+dghh - dfph)^2)/(48 listΩ0);


(*ω1 = -Ω2/Λ1*)
ω1 = - Omega2/Lambda1;

(*stability criterion: |ω1| < 1.110468 - 0: stable, 1: unstable*)
stabilitytest = Table[If[Im[ω1[[i,j]]] ==0,If[ω1[[i,j]]^2<1.233,0,1], -1],{i,1,Ny},{j,1,Nx}];

(*plot*)
ArrayPlot[stabilitytest,AxesLabel->Automatic,ColorRules->{1 -> Blue, 0 -> Red,-1-> Black},FrameTicks->{{Miny,Maxy},{Minx,Maxx}},DataReversed-> True,FrameLabel->{B1,E1}]
]

(*save data to generate the E1xB1 plot: E1B1Plots_Minx=*,Maxx=*,Miny=*,Maxy=*,beta=*,eta=*.txt*)
Export["E1B1Plots_Minx="<>ToString[Minx]<>",Maxx="<>ToString[Maxx]<>",Miny="<>ToString[Miny]<>",Maxy="<>ToString[Maxy]<>",beta="<>ToString[β]<>",eta="<>ToString[η]<>".txt" ,stabilitytest]