1
$\begingroup$

I am working on the paper "A short study of an SIR model with inclusion of an alert class, two explicit nonlinear incidence rates and saturated treatment rate" by Kumar, Nilam and Kishor and I am trying to plot the bifurcation diagram $I$ vs $R_0$ as they do in their paper. The plot is this enter image description here

And the system of ODES in the paper is this: enter image description here

I am trying to plot the graph, but I have no idea. I tried to mimic the method explained in this link by writing I* in terms of R_0, but it was not very helpful, since I obtained a polynomial of degree 3 in the variable I*.

Can you help me with this plotting? Thanks in advance. The value of R_0 given in the paper is enter image description here

Improve this question
$\endgroup$
5
  • 3
    $\begingroup$ Please make clear just what you are trying to do and show the code you have tried. $\endgroup$
    – bbgodfrey
    May 20, 2021 at 0:44
  • $\begingroup$ @bbgodfrey I am trying to plot the bifurcation diagram $R_0$ vs. $I$, but I have no idea about how to do it. The code I tried is useless now, since I couldn't put $I$ in terms of $R_0$ $\endgroup$
    – rowcol
    May 20, 2021 at 1:31
  • $\begingroup$ If making progress depends on the contents of the paper you mention at the beginning of your question, then I recommend that you provide access to that paper. Incidentally, in my experience, bifurcation diagrams look like the second figure in 96407, rather than the first figure in your question. $\endgroup$
    – bbgodfrey
    May 20, 2021 at 1:52
  • $\begingroup$ @bbgodfrey Oh I see. I thought that the name of that kind of graph were "bifurcation diagram" (since I googled it before haha). This is the paper I am working on: drive.google.com/file/d/1OE-ubS2KW092-fs8e8GskgkBQmuYRQEC/… $\endgroup$
    – rowcol
    May 20, 2021 at 2:05
  • $\begingroup$ @bbgodfrey I'd call OP's first figure a bifurcation diagram -- it's just a pretty boring one, with only a transcritical bifurcation at $R_0=1$. The paper seems to have another issue though, in that $R_0$ is not a parameter of the model per se but a combination of parameters. In the figure legend, all the model parameters are defined, so what is exactly being varied along the x-axis? Unclear. But that's more a problem with the paper than OP's question. $\endgroup$
    – Chris K
    May 20, 2021 at 2:51

1 Answer 1

Reset to default
2
$\begingroup$

As I mentioned in my comment, the paper is not clear what is exactly being varied to in Fig. 2, since all parameters have been given numerical values. I'll just go with pi as the hidden parameter being varied.

α = 0.5;
γ = 0.0009;
μ = 0.007;
d = 0.05;
δ = 0.002;
b = 0.2;
a = 0.2;
θ = 0.002;
β = 0.0012;

Solve for equilibria (ugly):

eq = Solve[{
   0 == pi - δ s - μ s - β s i/(1 + α i),
   0 == δ s - μ A - γ A i/(1 + α i),
   0 == β s i/(1 + α i) + γ A i/(1 + α i) - (μ + d + θ) i - a i/(1 + b i)
   }, {s, A, i}];
Length[eq]
(* 4 *)

Turns out the fourth equilibrium is the relevant one:

Plot[Chop[i /. eq[[4]]], {pi, 0, 4}]

enter image description here Now to get this as a function of R0, we define R0 then use ParametricPlot (with Max[0, ...] thrown in to avoid negative values of i):

R0 := pi (μ β + γ δ)/(μ (μ + δ) (d + μ + θ + a))

ParametricPlot[{R0, Max[0, Chop[i /. eq[[4]]]]}, {pi, 0, 4}, 
 PlotRange -> {{0, 2}, {-0.1, 4}}, AspectRatio -> 0.6, 
 AxesLabel -> {"R0", "I"}]

enter image description here

Evidently this isn't exactly what's plotted in the paper's figure 2, but the paper doesn't seem to give enough information to make exact recreation easy.

Improve this answer
$\endgroup$
1
  • 1
    $\begingroup$ Thank you very much; I started to make simulations and definitely it works. $\endgroup$
    – rowcol
    May 20, 2021 at 17:29

Not the answer you're looking for? Browse other questions tagged or ask your own question.