To make it easy, I'll use my [EcoEvo package][1]. First time, you'll need to install it: ``` PacletInstall["EcoEvo", "Site" -> "http://raw.githubusercontent.com/cklausme/EcoEvo/master"] ``` Then, load the package to get started: ``` <<EcoEvo` ``` Both models are a little bit funny, in that total population size is conserved. Thus, model (i) is effectively one-dimensional and model (ii) is effectively two-dimensional. Model (i): SI -------- ``` SetModel[{ Pop[pop] -> { Component[s] -> {Equation :> -β s i/n + γ i}, Component[i] -> {Equation :> β s i/n - γ i} }, Parameters :> {β > 0, γ > 0, n > 0} }] ``` Let's go straight to the phase planes using [PlotEcoPhasePlane][2], manually adding the total population constraint as a pink straight line. The no-disease case (disease-free equilibrium `eq[[1]]` is stable, the other equilibrium `eq[[2]]` is biologically meaningless since `s > n` and `i < 0`): ``` β = 0.95; γ = 1; n = 1; Show[ PlotEcoPhasePlane[{s, 0.9, 1.1}, {i, -0.1, 0.1}], Plot[n - s, {s, 0.9, 1.1}, PlotStyle -> Pink], RuleListPlot[{{s -> n, i -> 0}, {s -> (n γ)/β, i -> n - (n γ)/β}}, PlotMarkers -> {True, False}] ] ``` [![enter image description here][3]][3] and the endemic case (endemic equilibrium `eq[[2]]` is stable): ``` β = 4; γ = 1; n = 1; Show[ PlotEcoPhasePlane[{s, 0, 1}, {i, 0, 1}], Plot[n - s, {s, 0, n}, PlotStyle -> Pink], RuleListPlot[{{s -> n, i -> 0}, {s -> (n γ)/β, i -> n - (n γ)/β}}, PlotMarkers -> {False, True}] ] ``` [![enter image description here][4]][4] Again, because of the total-population constraint, this is effectively a one-dimensional system with either two or one feasible (non-negative) equilibria. Note that the two different isoclines ($S$ and $I$) overlap completely because of this and just look gold. Model (ii): SIR ------------- Here we can get rid of the degeneracy by defining `r := n - s - i`, then work in the SI phase-plane. ``` r := n - s - i; SetModel[{ Pop[pop] -> { Component[s] -> {Equation :> -β s i/n + ξ r}, Component[i] -> {Equation :> β s i/n - γ i} }, Parameters :> {β > 0, γ > 0, ξ > 0, n > 0} }] ``` To verify your analytical results: ``` eq = SolveEcoEq[] ``` [![enter image description here][5]][5] ``` EcoEigenvalues[eq[[1]]] ``` [![enter image description here][6]][6] ``` EcoEigenvalues[eq[[2]]] ``` [![enter image description here][7]][7] The eigenvalues of the non-trivial equilibrium are ugly, but we can check stability using Routh-Hurwitz criteria in [EcoStableQ][8]: ``` Simplify[EcoStableQ[eq[[2]]]] ``` [![enter image description here][9]][9] On to the phase planes. The $S$-isocline is blue, the $I$-isocline is gold. No-disease case (`eq[[1]]` is stable, `eq[[2]]` is a biologically meaningless saddle point): ``` β = 0.95; γ = 1; ξ = 1; n = 1; eq = SolveEcoEq[] N[EcoEigenvalues[eq[[2]]]] Show[ PlotEcoPhasePlane[{s, 0.9, 1.1}, {i, -0.1, 0.1}], RuleListPlot[eq, PlotMarkers -> EcoStableQ[eq]] ] (* {{s -> 1., i -> 0}, {s -> 1.05263, i -> -0.0263158}} *) (* {-1.02384, 0.0488359} *) ``` [![enter image description here][10]][10] Endemic case 1 (`eq[[2]]` is a stable focus, due to complex eigenvalues): ``` β = 4; γ = 1; ξ = 1; n = 1; eq = SolveEcoEq[] N[EcoEigenvalues[eq[[2]]]] Show[ PlotEcoPhasePlane[{s, 0, 1.1}, {i, 0, 1.1}], RuleListPlot[eq, PlotMarkers -> EcoStableQ[eq]] ] (* {{s -> 1, i -> 0}, {s -> 1/4, i -> 3/8}} *) (* {-1.25 + 1.19896 I, -1.25 - 1.19896 I} *) ``` [![enter image description here][11]][11] Endemic case 2 (`eq[[2]]` is a stable node, due to negative real eigenvalues): ``` β = 4; γ = 1; ξ = 10; n = 1; eq = SolveEcoEq[] N[EcoEigenvalues[eq[[2]]]] Show[ PlotEcoPhasePlane[{s, 0, 1.1}, {i, 0, 1.1}], RuleListPlot[eq, PlotMarkers -> EcoStableQ[eq]] ] (* {{s -> 1, i -> 0}, {s -> 1/4, i -> 15/22}} *) (* {-9.60337, -3.1239} *) ``` [![enter image description here][12]][12] [1]: https://github.com/cklausme/EcoEvo [2]: https://www.wolframcloud.com/obj/EcoEvo/docs/ref/PlotEcoPhasePlane.nb [3]: https://i.sstatic.net/kecl7m.png [4]: https://i.sstatic.net/kCuVjm.png [5]: https://i.sstatic.net/eGvLy.png [6]: https://i.sstatic.net/YoFfb.png [7]: https://i.sstatic.net/6VcNf.png [8]: https://www.wolframcloud.com/obj/EcoEvo/docs/ref/EcoStableQ.nb [9]: https://i.sstatic.net/l3S5d.png [10]: https://i.sstatic.net/riA6lm.png [11]: https://i.sstatic.net/wqDxhm.png [12]: https://i.sstatic.net/j3s5bm.png