# (In)Stability of equilibrium points in a vector plot (continued)?

I got a nice answer from Chris to my earlier question(though there are still some points unanswered): (In)Stability of equilibrium points depicted in a vector plot?

Now consider the following system:

\begin{align} \frac{dS}{dt} &= \nu N -\frac{\beta S I}{N} + \xi R - \nu S\\ \frac{dE}{dt} &= \frac{\beta S I}{N}- \sigma E -\nu E \\[2ex] \frac{dI}{dt} &= \sigma E -\gamma I -\nu I \\[2ex] \frac{dR}{dt} &= \gamma I -\xi R - \nu R \end{align}


where $$N = S + E + I + R$$ is the total population.

If we apply the same technique we have:

r := 1 - s - e - i;
SetModel[{
Pop[pop] -> {
Component[s] -> {Equation :> ξ + ν - ξ e - ξ i - β s i  - ξ s - ν s},
Component[e] -> {Equation :>  β s i - σ e - ν e },
Component[i] -> {Equation :> σ e - γ i - ν i }
},
Parameters :> {β > 0, γ > 0, ξ > 0, σ > 0, ν > 0}
}]
eq = SolveEcoEq[]


Which is fine. Now I want to do the following;

β = 0.7; γ = 0.3; σ = 0.3; ξ = 0.4; ν = 0.3;
eq = SolveEcoEq[]
N[EcoEigenvalues[eq[]]]
Show[PlotEcoPhasePlane[{s, 0, 2}, {e, -0.6, 1}],
RuleListPlot[eq, PlotMarkers -> EcoStableQ[eq]]]


Our system is 3 dimensional so how would I plot phase planes for say $$s$$ vs $$i$$ or $$s$$ vs $$e$$ etc? What about initial conditions for the phase portraits?

EDIT:

The paper:

https://www.mdpi.com/2227-7390/5/1/7/htm

Hopefully you can recreate the results!

My code for ODE solution:

tmax = 100;
\[Beta] = 4;
\[Sigma] = 1/3;
\[Gamma] = 1/7;
b = (1.9/100)/365;
d = (0.8/100)/365;
\[Alpha] = 1/14;
\[Nu] = 0.8;
\[Rho] = 0.5;
SExPRS = NDSolveValue[{
S'[t] == \[Alpha] +
b (1 - \[Nu]) - \[Alpha]*Ex[t] - \[Alpha]*P[t] - \[Beta]*S[t]*
P[t] - (b + \[Rho] + \[Alpha])*S[t],
Ex'[t] == \[Beta]*S[t]*P[t] - (b + \[Sigma])*Ex[t],
P'[t] == \[Sigma]*Ex[t] - (b + \[Gamma])*P[t],
S == 0.25,
Ex == 0.25,
P == 0.25},
{S, Ex, P},
{t, 0, tmax}];
{f1, f2, f3} = SExPRS;

st = Style[#, 15, Black] &;

Plot[{f1[t], f2[t], f3[t], 1 - f1[t] - f2[t] - f3[t], 1, 0}, {t, 0,
tmax}, PlotStyle -> {Blue, Orange, Red, Green,
Directive[Black, Dashed], Directive[Black, Dashed]}, Frame -> True,
FrameLabel -> st /@ {"Time", "Density"},
PlotLegends ->
Placed[LineLegend[{Blue, Orange, Red, Green, Black}, {"S(t)",
"E(t)", "I(t)", "R(t)" }, LegendFunction -> Framed], {0.85,
0.44}], ImageSize -> 500] r := 1 - s - e - i;
SetModel[{Pop[
pop] -> {Component[
s] -> {Equation :> \[Xi] + \[Nu] (1 -
c) - \[Xi] e - \[Xi] i - \[Beta] s i  - \[Xi] s - \[Nu] s \
- \[Rho] s},
Component[
e] -> {Equation :>  \[Beta] s i - \[Sigma] e - \[Nu] e },
Component[
i] -> {Equation :> \[Sigma] e - \[Gamma] i - \[Nu] i }},
Parameters :> {\[Beta] > 0, \[Gamma] > 0, \[Xi] > 0, \[Sigma] >
0, \[Nu] > 0, \[Rho] > 0, c > 0}}]
eq = SolveEcoEq[]
\[Beta] = 4; \[Gamma] = 1/7; \[Sigma] = 1/3; \[Xi] =
1/14; c = 0.8; \[Nu] = (1.9/100)/365; \[Rho] = 0.5;
eq = SolveEcoEq[]
N[EcoEigenvalues[eq[]]]
Show[PlotEcoPhasePlane[{s, 0, 1}, {i, 0, 1},
Fixed -> {e -> \[Sigma]}],
RuleListPlot[eq, PlotMarkers -> EcoStableQ[eq]]]
sol = EcoSim[{s -> 0.25, e -> 0.25, i -> 0.25}, 90];
PlotDynamics[sol] An example of $$e$$ plotted against $$i$$: • @Moo Nope, to be honest, I'm an amateur in Mathematica
– Math
Sep 14, 2021 at 11:12
• @Moo can you help?
– Math
Sep 15, 2021 at 9:36
• BTW, be sure to use :>, not -> after Equation in SetModel. Sep 19, 2021 at 14:32

Nice to see you trying my EcoEvo package. I will increase $$\beta$$ to get the disease to persist at a stable positive equilibrium.

β = 2;

eq = SolveEcoEq[]
EcoStableQ[eq]
(* {{s -> 1., e -> 0, i -> 0}, {s -> 0.6, e -> 0.233333, i -> 0.116667}} *)
(* {False, True} *)

sol = EcoSim[{s -> 0.999, e -> 0.001, i -> 0}, 100];
PlotDynamics[sol] PlotEcoPhasePlane takes the option Fixed, which you can use to fix other variables to get a 2D slice. For example,

PlotEcoPhasePlane[{s, 0, 1}, {e, 0, 1}, Fixed -> {i -> 0}] However, I'm not sure how useful this is, because the fixed variable i actually changes, which you can't see in such a 2D slice of the 3D phase space. One idea would be to assume a quasi-steady state between exposed e and infected i by setting $$di/dt=0$$, solving for $$i$$ in terms of $$e$$, then using that in your 2D phase plane:

PlotEcoPhasePlane[{s, 0, 1}, {e, 0, 1}, Fixed -> {i -> σ e/(γ + ν)}] Now we can see the stable equilibrium eq[] in a 2D slice.

I haven't built in support for 3D isoclines & streams yet, because I usually find them hard to see. But here's a start at it:

isoclines = ContourPlot3D[
Evaluate[EcoEqns[] /. ZeroLHS /. var_[t] -> var],
{s, 0, 1}, {e, 0, 1}, {i, 0, 1},
ContourStyle -> Opacity[0.8], Mesh -> None,
RegionFunction -> Function[{s, e, i}, s + e + i < 1],
AxesLabel -> {s, e, i}];

streams = StreamPlot3D[
Evaluate[EcoEqns[] /. RHS /. var_[t] -> var],
{s, 0, 1}, {e, 0, 1}, {i, 0, 1}, AxesLabel -> {s, e, i},
RegionFunction -> Function[{s, e, i}, s + e + i < 1],
StreamStyle -> {Black, Thin}, StreamColorFunction -> None,
StreamPoints -> Fine];

Show[isoclines, streams] Given your interests, I suggest including $$R$$ as a dynamic variable. You can set colors for each variable with Color in SetModel.

SetModel[{
Pop[pop] -> {
Component[s] -> {Equation :> ν (1 - c) - β s i - ν s + ξ r - ρ s, Color -> Blue},
Component[e] -> {Equation :> β s i - σ e - ν e, Color -> Orange},
Component[i] -> {Equation :> σ e - γ i - ν i, Color -> Red},
Component[r] -> {Equation :> ν c + γ i - ν r - ξ r + ρ s, Color -> Green}
}, Parameters :> {β > 0, γ > 0, ξ > 0, σ > 0, ν > 0, ρ > 0, c > 0}}]


To simulate:

β = 4; γ = 1/7; σ = 1/3; ξ = 1/14; c = 0.8; ν = (1.9/100)/365; ρ = 0.5;
sol = EcoSim[{s -> 0.25, e -> 0.25, i -> 0.25, r -> 0.25}, 100];
PlotDynamics[sol, PlotRange -> {0, 1}] To plot the dynamics of one variable vs another, use RuleListPlot.

RuleListPlot[sol, {e, i}] Here's a grid of all pairwise combinations:

GraphicsGrid[
Table[RuleListPlot[sol, {var1, var2}], {var1, {s, e, i, r}}, {var2, {s, e, i, r}}],
Spacings -> 0] • Nice answer yet again! Is there a way to plot say $s$ vs $r$? I am reading a paper(their model is a little more complex) but very similar. I can attach it in the question as an edit if you'd like to see. They have managed to produce results that I am unable to produce(albeit quite close), maybe you can try?
– Math
Sep 20, 2021 at 14:48
• Sure, attach the paper and I'll take a look when I have time. Sep 20, 2021 at 14:56
• It will give me coding errors again for some reason but I have attached it. The main interest is the solution to the ODE system and the phase plots.
– Math
Sep 20, 2021 at 15:07
• BTW you need to use  to mark LaTeX code Sep 20, 2021 at 15:32
• Based on the paper, I assume you don't want isoclines or streams, just plotting a single trajectory. Since r is only defined algebraically, the easiest thing to do is ParametricPlot[ Evaluate[{s[t], 1 - s[t] - e[t] - i[t]} /. sol], {t, 0, 100}, AspectRatio -> 1]. For variables that are listed as Components you could use RuleListPlot[sol, {s, e}], so another approach would be to explicitly model r as well. Sep 20, 2021 at 15:50