Trying to model red and grey squirrels

  • Context: Mathematica 
  • Thread starter Thread starter LETS_GO
  • Start date Start date
  • Tags Tags
    Model population sir
Click For Summary
SUMMARY

This discussion focuses on modeling the population dynamics of red and grey squirrels using a system of differential equations in Mathematica. The model incorporates variables such as natural mortality rates, virus transmission rates, and reproductive rates for both species. Key parameters include a natural mortality rate of 0.4, a virus transmission rate of 0.7, and specific carrying capacities for reds (60) and greys (80). The conversation highlights the need to define the maximum reproductive rate of greys as a function of time to accurately simulate population changes.

PREREQUISITES
  • Understanding of differential equations and population dynamics
  • Familiarity with Mathematica and its NDSolve function
  • Knowledge of SIR (Susceptible, Infected, Recovered) models
  • Basic concepts of ecological modeling and competition between species
NEXT STEPS
  • Implement time-dependent functions in Mathematica for population modeling
  • Explore advanced features of NDSolve for dynamic systems
  • Research ecological modeling techniques for species competition
  • Study the effects of varying reproductive rates on population stability
USEFUL FOR

Ecologists, mathematicians, and researchers interested in population dynamics and mathematical modeling of species interactions, particularly those focusing on wildlife management and disease impact studies.

LETS_GO
Messages
4
Reaction score
0
TL;DR
Hello, we are trying to change the max reproductive rate of grey squirrels. with a value starting at 1.2 and want it to go down to 0.3 with a step of 0.05.
Code: 
ClearAll["Global`*"]
(*R = reds, G = greys*)
(*S = susceptible, I = infected, R = recovered*)

tseries = {t, 0, 3};
vars = {HG[t], HR[t], SG[t], IG[t], RG[t], SR[t], IR[t], aG[t], qG[t]};

b = 0.4;    (*natural mortality rate, both species*)
\[Beta] = 0.7;    (*rate of virus transmission*)

aR = 1;       (*Reds max. reproductive rate*)
\[Alpha] = 26;       (*Reds mortaility rate due to virus*)
cR = 0.61;(*Reds competative effect on greys*)
KR = 60;   (*Reds carrying capacity, 5x5 km*)

\[Gamma] = 13;      (*Greys recovery rate from virus*)
cG = 1.65;(*Greys competative effect on reds*)
KG = 80;   (*Greys carrying capacity, 5x5 km*)

qR = (aR - b)/KR;  (*Reds crowding susceptibility*)
  (*Greys crowding susceptibility*)

(*initial conditions*)
SGinit = 10;
IGinit = 2;
RGinit = 0;
SRinit = 60;         
IRinit = 0;
HGinit = SGinit + IGinit + RGinit;
HRinit = SRinit + IRinit;

 eqns =
    (*total populations*)
  {HG[t] == SG[t] + IG[t] + RG[t],
   HR[t] == SR[t] + IR[t],
 
   aG[t] == aG[t - 1] - 0.05,(*aG = Greys max. reproductive rate -
   possible birth control??*)
   qG[t] == (aG[t] - b)/KG,  (*Greys crowding susceptibility*)
   (*SIR growth rates*)
 
   SG'[t] == ((aG[t] - (qG[t]*(HG[t] + (cR*HR[t]))))*HG[t]) - (b*
       SG[t]) - (\[Beta]*SG[t]*(IG[t] + IR[t])),
   IG'[t] == (\[Beta] *SG[t]*(IG[t] + IR[t])) - (b*IG[t]) - (\[Gamma]*
       IG[t]),
   RG'[t] == (\[Gamma]*IG[t]) - (b*RG[t]),
   SR'[t] == ((aR - (qR*(HR[t] + (cG*HG[t]))))*HR[t]) - (b*
       SR[t]) - ((\[Beta]*SR[t])*(IR[t] + IG[t])),
   IR'[t] == (\[Beta]*SR[t]*(IG[t] + IR[t])) - ((\[Alpha] + b)*IR[t]),
 
 
   (*call initial conditions*)
   HG[0] == HGinit, HR[0] == HRinit,
   aG[0] == 1.2, qG[0] == (1.2 - b)/KG,
   SG[0] == SGinit, IG[0] == IGinit, RG[0] == RGinit,
   SR[0] == SRinit,
   IR[0] ==
    IRinit                                                            \
        };

sol = NDSolve[eqns, vars, tseries];

{Plot[Evaluate[{SG[t], IG[t], RG[t], SR[t], IR[t]} /. sol], tseries,
  ImageSize -> Large, PlotLegends -> {"SG", "IG", "RG", "SR", "IR"}],
 Plot[Evaluate[{HG[t], HR[t]} /. sol], tseries, ImageSize -> Large,
  PlotLegends -> {"HG", "HR"}]}
<Moderator's note: Please use CODE tags when posting code.>
 
Last edited by a moderator:
Physics news on Phys.org
What is your question?
 
how do we change aG with time
 
You need to make it an actual function of time
Code: 
aG[t] == a0 + a1 * t
with appropriate constants a0 and a1.
 
LETS_GO said:
how do we change aG with time
[CODE title="It looks like it already is changing with each time step on line 37:"]aG[t] == aG[t - 1] - 0.05,(*aG = Greys max. reproductive rate -
possible birth control??*)[/CODE]
That looks like the change that you want. I'm not familiar with NDsolve and don't see anything wrong with your code.
 
FactChecker said:
[CODE title="It looks like it already is changing with each time step on line 37:"]aG[t] == aG[t - 1] - 0.05,(*aG = Greys max. reproductive rate -
possible birth control??*)[/CODE]
That looks like the change that you want. I'm not familiar with NDsolve and don't see anything wrong with your code.
That doesn't work, because NDSolve does not iterate the solution like that. It needs to be a proper function of time. (And note that t is a real variable, not an index, so t - 1 means "one unit of time earlier".)
 
Reply
  • Like
Likes   Reactions: FactChecker
DrClaude said:
You need to make it an actual function of time
Code: 
aG[t] == a0 + a1 * t
with appropriate constants a0 and a1.

So I want my starting value for aG to be 1.2 and I want it to go down to at least 0.35 with a decrease of 0.5 for each time step
 
LETS_GO said:
So I want my starting value for aG to be 1.2 and I want it to go down to at least 0.35 with a decrease of 0.5 for each time step
But there is no time step! The problem is not fundamentally discrete, but continuous. You need to figure out what aG is as a function of time.
 
hmmm I have tried the last couple of days, and im really not sure how to do this