- #1

- 4

- 0

- TL;DR Summary
- 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"}]}
```

Last edited by a moderator: