Twocomponent model with p = const
The hostvirus dynamics may be considered within the framework of the Volterratype models
\begin{array}{l}\mathit{x}`=\mathit{x}(1\mathit{\text{ax}}\mathit{bz}\left(1\mathit{p}\right))\equiv \mathit{P}(\mathit{x},\mathit{z}),\\ \mathit{z}`=\mathit{z}(\mathit{d}\mathit{\text{bxp}}+\mathit{\text{bMx}}\left(1\mathit{p}\right))\equiv \mathit{Q}(\mathit{x},\mathit{z})\\ \phantom{\rule{.3em}{0ex}}\mathit{x}\ge 0,\phantom{\rule{0.36em}{0ex}}\mathit{z}\ge 0;\mathit{a}\ge 0,\mathit{b},\mathit{d},\mathit{M}>0\end{array}
(M1)
If the parameter p is constant and a = 0 then model (M1) is Hamiltonian (conservative) with the Hamiltonian G(x, z) = ln z + d ln x  b(M(1  p)  p)x  b(1  p)z. If 0<\mathit{p}<\frac{\mathit{M}}{\mathit{M}+1} then the system has a saddle in the origin O and a center in the equilibrium
\mathit{B}\left(\overline{\mathit{x}}=\frac{\mathit{d}}{\mathit{b}\left(\mathit{M}\left(1\mathit{p}\right)\mathit{p}\right)},\overline{\mathit{z}}=\frac{1\mathit{a}\overline{\mathit{x}}}{\mathit{b}\left(1\mathit{p}\right)}\right).
If a > 0 (logistic case), then for constant 0<\mathit{p}<\frac{\mathit{M}}{\mathit{M}+1} system (M1) has a saddle in the origin, equilibria A(1/a, 0) and \mathit{B}\left(\mathit{x}=\frac{\mathit{d}}{\mathit{b}\left(\mathit{M}\left(1\mathit{p}\right)\mathit{p}\right)},\mathit{z}=\frac{\mathit{b}\left(\mathit{M}\left(1\mathit{p}\right)\mathit{p}\right)\mathit{\text{ad}}}{{\mathit{b}}^{2}\left(1\mathit{p}\right)\left(\mathit{M}\left(1\mathit{p}\right)\mathit{p}\right)}\right); equilibrium B belongs to the positive quadrant (x, z) if b(M(1  p)  p)  ad > 0. In this case, B is a stable node/spiral and A is a saddle , A is a stable node if B is not positive (see, e.g., [54]).
Twocomponent model with p = p(z)
Malthusian version of the model (M1) (a = 0).
This system has equilibrium in the origin (0,0), which is a saddle; z coordinate of any other equilibrium \left(\overline{\mathit{x}},\overline{\mathit{z}}\right) has to be a root of the equation 1  bz(1  p(z)) = 0 and \overline{\mathit{x}}=\frac{\mathit{d}}{\mathit{b}\left(\mathit{M}\mathit{p}\left(\overline{\mathit{z}}\right)\left(1+\mathit{M}\right)\right)}. The point \left(\overline{\mathit{x}},\overline{\mathit{z}}\right) is positive only if \mathit{p}\left(\overline{\mathit{z}}\right)<\mathit{M}/\left(1+\mathit{M}\right).
Proposition 1. If p_{
z
}(z) < 0, then the Malthusian model has only one nontrivial equilibrium \mathit{B}\left(\overline{\mathit{x}},\overline{\mathit{z}}\right),which is an unstable node/spiral for every parameter values.
Proof. The Jacobian of the system in the equilibrium \left(\overline{\mathit{x}},\overline{\mathit{z}}\right) is
\mathit{J}\left(\overline{\mathit{x}},\overline{\mathit{z}}\right)=\left(\begin{array}{cc}0& \mathit{b}\overline{\mathit{x}}(1+\mathit{p}\left(\overline{\mathit{z}}\right)+\overline{\mathit{z}}{\mathit{p}}_{\mathit{z}}(\overline{\mathit{z}}))\\ \mathit{b}\overline{\mathit{z}}(\mathit{M}+\left(1+\mathit{M}\right)\mathit{p}(\overline{\mathit{z}}))& \mathit{b}\left(1+\mathit{M}\right)\overline{\mathit{x}}\overline{\mathit{z}}{\mathit{p}}_{\mathit{z}}(\overline{\mathit{z}})\end{array}\right),
and its determinant and trace are:
\begin{array}{l}\mathit{\text{Det}}(\mathit{J}\left(\overline{\mathit{x}},\overline{\mathit{z}}\right))={\mathit{b}}^{2}\overline{\mathit{x}}\overline{\mathit{z}}(\mathit{M}(1+\mathit{M})\mathit{p}(\overline{\mathit{z}}))(1\mathit{p}(\overline{\mathit{z}})\overline{\mathit{z}}{\mathit{p}}_{\mathit{z}}(\overline{\mathit{z}})),\\ \begin{array}{l}\mathit{\text{Trace}}(\mathit{J}\left(\overline{\mathit{x}},\overline{\mathit{z}}\right))=\mathit{b}\left(1+\mathit{M}\right)\overline{\mathit{x}}\overline{\mathit{z}}{\mathit{p}}_{\mathit{z}}(\overline{\mathit{z}})\end{array}.\end{array}
(M2)
If the function p(z) decreases monotonically, i.e. p_{
z
}(z) < 0, then p(z) and monotonically increasing function \mathit{h}\left(\mathit{z}\right)=1\frac{1}{\mathit{\text{bz}}} can intersect only once. Thus, equation 1  bz(1  p(z)) = 0 has only one root \overline{\mathit{z}}, and the Malthusian model has only one nontrivial equilibrium \mathit{B}\left(\overline{\mathit{x}},\overline{\mathit{z}}\right). Next, \mathit{\text{Det}}\left(\mathit{J}\left(\overline{\mathit{x}},\overline{\mathit{z}}\right)\right)>0,\mathit{\text{Tr}}\phantom{\rule{0.24em}{0ex}}\left(\mathit{J}\left(\overline{\mathit{x}},\overline{\mathit{z}}\right)\right)>0 for positive \overline{\mathit{x}},\overline{\mathit{z}} if p_{
z
}(z) < 0. Thus \mathit{B}\left(\overline{\mathit{x}},\overline{\mathit{z}}\right) is unstable node or unstable spiral.
Logistic version of model (M1) (a = 1).
Equilibria of system (M1) with a = 1 and p = p(z) defined by (3) are the points (0,0), A(1,0), and \mathit{B}\left(\overline{\mathit{x}},\overline{\mathit{z}}\right) where coordinates \overline{\mathit{x}},\overline{\mathit{z}} solve the system
\begin{array}{l}\overline{\mathit{x}}=\frac{\mathit{d}}{\mathit{b}\left(\mathit{M}\mathit{p}\left(\overline{\mathit{z}}\right)\left(1+\mathit{M}\right)\right)},1\overline{\mathit{x}}\mathit{b}\overline{\mathit{z}}\left(1\mathit{p}\left(\overline{\mathit{z}}\right)\right)=0,\\ \mathit{p}\left(\overline{\mathit{z}}\right)<\mathit{M}/\left(\mathit{M}+1\right).\end{array}
(M3)
Denote \mathit{h}\left(\mathit{z}\right)\equiv \frac{\mathit{d}\mathit{\text{bM}}+{\mathit{b}}^{2}\mathit{\text{Mz}}}{\mathit{b}\left(1\mathit{M}+\mathit{\text{bMz}}\right)}, then \overline{\mathit{z}} is a root of the equation p(z) = h(z). Solutions of this equation are the points of intersection of the curves p(z) and h(z); up to two equilibrium points {\mathit{B}}_{1}\left({\overline{\mathit{x}}}_{1},{\overline{\mathit{z}}}_{1}\right),{\mathit{B}}_{2}\left({\overline{\mathit{x}}}_{2},{\overline{\mathit{z}}}_{2}\right) can appear in the model with parameter variations. Denote J(x, z) the Jacobian of system (M1), (3) with a = 1:
\mathit{J}\left(\mathit{x},\mathit{z}\right)=\left(\begin{array}{cc}12\mathit{x}\mathit{bz}(1\mathit{p}\left(\mathit{z}\right))& \mathit{bx}(1\mathit{p}\left(\mathit{z}\right)\mathit{z}{\mathit{p}}_{\mathit{z}}(\mathit{z}))\\ \mathit{bz}(\mathit{M}+\left(1+\mathit{M}\right)\mathit{p}(\mathit{z}))& \mathit{d}+\mathit{bMx}(1\mathit{p}\left(\mathit{z}\right))\mathit{bxp}(\mathit{z})\mathit{b}(1+\mathit{M})\mathit{xz}{\mathit{p}}_{\mathit{z}}(\mathit{z})\end{array}\right),
\begin{array}{l}\mathit{J}\left(0,0\right)=\left(\begin{array}{cc}1& 0\\ 0& \mathit{d}\end{array}\right),\\ \mathit{J}\left(1,0\right)=\left(\begin{array}{cc}1& \mathit{b}(1\mathit{p}\left(0\right))\\ 0& \mathit{d}+\mathit{\text{bM}}(1\mathit{p}\left(0\right))\mathit{\text{bp}}(0)\end{array}\right)=\left(\begin{array}{cc}1& 0\\ 0& \mathit{d}\mathit{b}\end{array}\right).\end{array}
Thus O(0,0) is a saddle and A(1,0) is a stable node for all parameter values.
Now let us consider the system consisting of equation (M3) supplemented with the equation \mathit{\text{Det}}\left(\mathit{J}\left(\overline{\mathit{x}},\overline{\mathit{z}}\right)\right)=0. In the case {\mathit{\lambda}}_{1}=\mathit{\text{Trace}}\left(\mathit{J}\left(\overline{\mathit{x}},\overline{\mathit{z}}\right)\right)\ne 0 this system defines coordinates of saddlenode point \mathit{B}\left(\overline{\mathit{x}},\overline{\mathit{z}}\right) whose second eigenvalue λ_{2} = 0. The last system supplemented with the equation \mathit{\text{Trace}}\left(\mathit{J}\left(\overline{\mathit{x}},\overline{\mathit{z}}\right)\right)=0 defines an additional degeneration in the system in the point \mathit{B}\left(\overline{\overline{\mathit{x}}},\overline{\overline{\mathit{z}}}\right) where λ_{1} = λ_{2} = 0.
Lemma 1. For a wide range of fixed parameters b, d, s there exist a point \left(\overline{\mathit{M}},\overline{\mathit{k}}\right)in the (M, k) parameter plane such that the BogdanovTakens bifurcation of codimension 2 is realized in the model (M1), a = 1 under variations of M and k close to \left(\overline{\mathit{M}},\overline{\mathit{k}}\right).The values \left(\overline{\mathit{M}},\overline{\mathit{k}}\right) and coordinates of the equilibrium \mathit{B}\left(\overline{\overline{\mathit{x}}},\overline{\overline{\mathit{z}}}\right)where \phantom{\rule{0.24em}{0ex}}\overline{\overline{\mathit{x}}}=\mathit{x}\left(\overline{\mathit{M}},\overline{\mathit{k}}\right),\overline{\overline{\mathit{z}}}=\mathit{z}\left(\overline{\mathit{M}},\overline{\mathit{k}}\right)are defined by the system consisting of equations (M3) and equations \mathit{\text{Det}}\left(\mathit{J}\left(\overline{\mathit{x}},\overline{\mathit{z}}\right)\right)=0, \mathit{Trace}\left(\mathit{J}\left(\overline{\mathit{x}},\overline{\mathit{z}}\right)\right)=0(see, e.g. , [41]).
We have found this bifurcation for some reasonable fixed values of the parameters d, b, s with the help of computer package LOCBIF [55]. Using the Lemma and Proposition 1 we prove the following statement.
Proposition 2. (1) System ( M1 ), (3) has equilibria: the saddle O(0,0) and stablenode A(1, 0) for all positive values of the parameters b, d, M, k = 1, s;
(2) there exists positive M* such that

a)
for M < M* the system has only the equilibria O and stable equilibrium A;

b)
for M > M* the system has two more equilibria, a saddle {\mathit{B}}_{1}\left({\overline{\mathit{x}}}_{1},{\overline{\mathit{z}}}_{1}\right) and a stable topological node/spiral
\phantom{\rule{0.12em}{0ex}}{\mathit{B}}_{2}\left({\overline{\mathit{x}}}_{2},{\overline{\mathit{z}}}_{2}\right);

c)
there exist M ^{* *} > M* such that for M > M ^{* *} the spiral \phantom{\rule{0.12em}{0ex}}{\mathit{B}}_{2}\left({\overline{\mathit{x}}}_{2},{\overline{\mathit{z}}}_{2}\right) is placed inside an unstable limit cycle.
The bifurcation diagram of the system under variation of the parameter M is shown in Additional file 1.
Threecomponent model: proof of Statement 1
Let us formulate Statement 1 in more details:

(1)
Equilibrium O(x = y = z = 0) has eigenvalues λ _{1}(O) = 1, λ _{2}(O) =  d, λ _{3}(O) = 1  l; it is unstable for all parameter values;

(2)
If a > 0 then equilibrium A(0,1/a,0) has eigenvalues λ _{1}(A) =  l, λ _{2}(A) =  1, \phantom{\rule{0.12em}{0ex}}{\mathrm{\lambda}}_{3}\left(\mathrm{A}\right)=\frac{\mathit{da}+\mathit{b}\left(\mathit{M}\left(1\mathit{s}\right)\mathit{s}\right)}{\mathit{a}}; it is stable for \mathit{M}<\frac{\mathit{da}+\mathit{bs}}{\mathit{b}\left(1\mathit{s}\right)} and unstable for \mathit{M}>\frac{\mathit{da}+\mathit{bs}}{\mathit{b}\left(1\mathit{s}\right)}.
Proof. Jacobian of system (1) is of the form: J(x, y, z) = (a_{i,j})i, j = 1, 2, 3, where
a_{1,1} = 1  l  ax  a(x + y)  bz(1  p(z)), a_{1,2} =  ax + esz, a_{1,3} = esy  bx(1  p(z)) + bxzp_{
z
}(z); a_{2,1} = l  ay, a_{2,2} = 1  ay  a(x + y)  b(1  s)z  esz, a_{2,3} =  b(1  s)y  esy,
{\mathit{a}}_{3,1}=\mathit{bMz}\left(1\mathit{p}\left(\mathit{z}\right)\right)\mathit{bzp}\left(\mathit{z}\right),{\mathit{a}}_{3,2}=\mathit{bM}\left(1\mathit{s}\right)\mathit{z}\mathit{bsz},
{\mathit{a}}_{3,3}=\mathit{d}+\mathit{b}\left(\mathit{M}\left(\mathit{x}+\mathit{y}\right)\left(\mathit{M}+1\right)\left(\mathit{p}\left(\mathit{z}\right)\mathit{x}+\mathit{sy}\right)\right)\mathit{bxz}{\mathit{p}}_{\mathit{z}}\left(\mathit{z}\right)\left(1+\mathit{M}\right).
p_{
z
}(z) = 0, if p = const and p_{
z
}(z) = k(s  p(z)), if p(z) = (1  s)e^{ kz} + s.
\begin{array}{l}\mathit{J}\left(0,0,0\right)=\left(\begin{array}{ccc}1\mathit{l}& 0& 0\\ \mathit{l}& 1& 0\\ 0& 0& \mathit{d}\end{array}\right),\\ \mathit{J}\left(0,1/\mathit{a},0\right)=\left(\begin{array}{ccc}\mathit{l}& 0& \mathit{bes}/\mathit{a}\\ 1+\mathit{l}& 1& \mathit{b}\left(1+\mathit{s}\mathit{es}\right)/\mathit{a}\\ 0& 0& (\mathit{ad}\mathit{b}\left(\mathit{M}\left(1\mathit{s}\right)\mathit{s}\right)/\mathit{a}\end{array}\right).\end{array}
So, the equilibrium O(0, 0, 0) is unstable for both Malthusian and logistic models, the equilibrium \mathit{A}\left(0,\frac{1}{\mathit{a}},0\right) of logistic model is stable for \mathit{M}<\frac{\mathit{ad}+\mathit{bs}}{\mathit{b}\left(1\mathit{s}\right)} and unstable for \mathit{M}>\frac{\mathit{ad}+\mathit{bs}}{\mathit{b}\left(1\mathit{s}\right)}.The Statement is proven.
Threecomponent Malthusian model with constant immunity p, p ≥ s > 0
Let us analyze the dynamics of model (1) depending on the parameters l, e. The system has trivial unstable equilibrium O(0,0,0). Let firstly p = s. Then system (1) by changing of variables u = x + y, z = z is reduced to 2component Malthusian system (M1) with respect to u, z  variables. For \mathit{s}=\mathit{p}<\frac{\mathit{M}}{\mathit{M}+1} the orbits {u, z} of this system belongs to the positive quadrant, the nontrivial equilibrium \overline{\mathit{u}}=\frac{\mathit{d}}{\mathit{b}\left(\mathit{M}\left(1\mathit{s}\right)\mathit{s}\right)},\overline{\mathit{z}}=\frac{1}{\mathit{b}\left(1\mathit{s}\right)} is a center, i.e., it is located inside closed orbits (similar to model (M1)) and trajectories demonstrate periodic oscillations. In variables x(t), y(t), z(t) model (1) also demonstrates periodic oscillations; the model has equilibrium B(x_{
e
}, y_{
e
}, z_{
e
}) where
\begin{array}{l}{\mathit{x}}_{\mathit{e}}=\frac{\mathit{des}}{\mathit{b}\left(\mathit{M}\left(1\mathit{s}\right)\mathit{s}\right)\left(\mathit{bl}\left(1\mathit{s}\right)+\mathit{es}\right)},\\ {\mathit{y}}_{\mathit{e}}=\frac{\mathit{dl}\left(1\mathit{s}\right)}{\left(\mathit{M}\left(1\mathit{s}\right)\mathit{s}\right)\left(\mathit{bl}\left(1\mathit{s}\right)+\mathit{es}\right)},\\ {\mathit{z}}_{\mathit{e}}=\frac{1}{\mathit{b}\left(1\mathit{s}\right)}\end{array}
(M4)
If p > s then any nontrivial equilibrium (x_{
e
}, y_{
e
}, z_{
e
}) of model (1) is such that the coordinate z_{
e
} solves the quadratic equation
\left(1\mathit{l}\right)+\mathit{b}\left(2+\mathit{l}+\mathit{p}+\mathit{s}\mathit{es}\left(1\mathit{l}\right)\right)\mathit{z}+{\mathit{b}}^{2}\left(1\mathit{p}\right)\left(1\mathit{s}+\mathit{es}\right){\mathit{z}}^{2}=0,
(M5)
and coordinates x_{
e
}, y_{
e
} can be expressed via z = z_{
e
} as follows:
{\mathit{x}}_{\mathit{e}}=\frac{\mathit{d}\left(1\mathit{b}\left(1\mathit{s}\right)\mathit{z}\right)}{\mathit{b}\left(\mathit{p}\mathit{s}\right)\left(1+\mathit{M}\mathit{bz}\right)},\phantom{\rule{0.48em}{0ex}}{\mathit{y}}_{\mathit{e}}=\frac{\mathit{d}\left(1\mathit{b}\left(1\mathit{p}\right)\mathit{z}\right)}{\mathit{b}\left(\mathit{p}\mathit{s}\right)\left(1+\mathit{M}\mathit{bz}\right)}.
(M6)
Let \mathit{l}\left(\mathit{e}\right)=\left(\mathit{M}+\mathit{p}+\mathit{Mp}\right)\left(1+\frac{\left.\mathit{s}(1+\mathit{M}\right)}{\left.\mathit{M}\mathit{s}(1+\mathit{M}\right)}\mathit{e}\right). This equation defines a boundary line Γ = (e, l(e)) in the parametric plane (e, l).
Proposition 3.

1)
System (1) with p = const and a = 0 has at most one positive equilibrium B(x _{
e
}, y _{
e
}, z _{
e
}) where x _{
e
}, y _{
e
}, z _{
e
} are defined by formulas (M5), (M6) for s < p and by formula (M4) for s = p;

2)
the system has a single positive equilibrium B (x _{
e
}, y _{
e
}, z _{
e
}) if and only if one of the following conditions holds:

a)
\mathit{s}<\mathit{p}<\frac{\mathit{M}}{\mathit{M}+1};

b)
\mathit{s}<\frac{\mathit{M}}{\mathit{M}+1}<\mathit{p} and l > l(e);

3)
in the case a) the equilibrium B is asymptotically stable;

4)
for s = p the system demonstrates periodic oscillations of variables x ( t ), y(t) while z(t) tends to a stable value for a wide domain of initial values (x _{
0,
} y _{0,} z _{0}) close to the equilibrium B.
Proof.
Taking the solutions of quadratic equation (M5) in the form
\begin{array}{l}\mathit{z}={\mathit{z}}_{1}\left(\mathit{l},\mathit{e}\right)=\frac{2\mathit{l}\left(\mathit{p}+\mathit{s}\right)+\mathit{es}(1\mathit{l})+\sqrt{\mathit{D}}}{2\mathit{b}\left(1\mathit{p}\right)\left(1\mathit{s}+\mathit{es}\right)},\\ \mathit{z}={\mathit{z}}_{2}\left(\mathit{l},\mathit{e}\right)=\frac{2\mathit{l}\left(\mathit{p}+\mathit{s}\right)+\mathit{es}(1\mathit{l})+\sqrt{\mathit{D}}}{2\mathit{b}\left(1\mathit{p}\right)\left(1\mathit{s}+\mathit{es}\right)}\end{array}
where D = l^{2}(1  s)^{2} + (p  s + es)^{2}  2l(p(1  s + 2es)  s(1 + e  s + es) we can easily verify that both "branches" z_{1}(l, e), z_{2}(l, e) are real for any positive (e, l) because the expression under the radical is nonnegative. The branches z_{1}(l, e), z_{2}(l, e) are positive both if l < 1 and only z_{1}(l, e) is positive if l > 1.
Analysis of formulas (M5), (M6) shows that only the branch z_{1}(l, e) can define positive coordinates of the equilibrium x_{
e
} = x(z_{1}), y_{
e
} = y_{
e
}(z_{1}). Substituting z_{1}(l, e) into formulas (M6) we obtain that x_{
e
}(e, l), y_{
e
}(e, l) are positive if the point (e, l) in the parametric plane is placed above the boundary line Γ given by equation
\begin{array}{l}\left(1+\mathit{p}\right)\left(1+\left(1+\mathit{e}\right)\mathit{s}\right)(\mathit{l}\left(\mathit{M}\left(1+\mathit{s}\right)+\mathit{s}\right)\\ \phantom{\rule{1.75em}{0ex}}+\left(\mathit{M}\left(1+\mathit{p}\right)+\mathit{p}\right)\left(\left(1+\mathit{e}\right)\mathit{s}+\mathit{M}\left(1+\left(1+\mathit{e}\right)\mathit{s}\right)\right)=0.\end{array}
Statements 1 and 2 of the Proposition are proven.
Let us analyze a stability of equilibrium B (x_{
e
}, y_{
e
}, z_{
e
}) of the system. For p = s characteristic polynomial of the system in the point B, whose coordinates are given by (M4), is of the form \mathit{E}\left({\mathit{\mu}}_{0}\right)=\mathit{Det}(\left(\mathit{J}\left(\mathit{B}\right){\mathit{\mu}}_{0}\mathit{I}\right)=\frac{\left(\mathit{d}+{\mathit{\mu}}_{0}^{2}\right)\left(\mathit{b}\left(\mathit{l}+{\mathit{\mu}}_{0}\right)\left(1\mathit{s}\right)+\mathit{es}\right)}{\mathit{b}\left(1\mathit{s}\right)}. Thus, two eigenvalues of the point are imaginary, \phantom{\rule{0.25em}{0ex}}{{\mathit{\mu}}_{0}}^{1,2}=\pm \mathit{i}\sqrt{\mathit{d}},\phantom{\rule{0.25em}{0ex}} and the third is negative, \phantom{\rule{0.25em}{0ex}}{{\mathit{\mu}}_{0}}^{3}=\frac{\mathit{es}+\mathit{bl}\left(1\mathit{s}\right)}{\mathit{b}\left(1\mathit{s}\right)}.Thus, statement 4 is proven.
We show now that for p > s the point B is a sink, i.e., its eigenvalues have negative real parts (more exactly, one eigenvalue is real negative, and two others are complex with negative real part).
Introduce the parameter α = p  s and write the right hands of (1), a = 0 in the form:
\begin{array}{l}\mathit{P}\left(\mathit{x},\mathit{y},\mathit{z}\right)=\mathit{x}\mathit{lx}\mathit{b}(1\mathit{\alpha}\mathit{s})\mathit{xz}+\mathit{esyz},\\ \mathit{Q}\left(\mathit{x},\mathit{y},\mathit{z}\right)=\mathit{lx}+\mathit{y}\mathit{b}(1\mathit{s})\mathit{yz}\mathit{esyz},\\ \mathit{R}\left(\mathit{x},\mathit{y},\mathit{z}\right)=\mathit{z}(\mathit{d}+\mathit{bM}((1\mathit{\alpha}\mathit{s})\mathit{x}+(1\mathit{s})\mathit{y})\\ \phantom{\rule{6.5em}{0ex}}\mathit{b}((\mathit{\alpha}+\mathit{s})\mathit{x}+\mathit{sy})).\end{array}
From the condition R(x, y, z) = 0, Q(x, y, z) = 0 we can express x = x_{
e
}, y = y_{
e
} via z = z_{
e
}:
\begin{array}{l}{\mathit{x}}_{\mathit{e}}=\frac{(\mathit{d}\left(1\mathit{b}\left(1\mathit{s}\right)\mathit{z}+\mathit{esz}\right)}{\mathit{b}(\mathit{\alpha}\left(1+\mathit{M}\right)\left(1\mathit{b}\left(1\mathit{s}\right)\mathit{z}+\mathit{esz}\right)+\left(\mathit{M}\left(1\mathit{s}\right)\mathit{s}\right)\left(1+\mathit{l}+\mathit{esz}+\mathit{bz}\left(1\mathit{s}\right)\right)}.\\ {\mathit{y}}_{\mathit{e}}=\frac{\mathit{dl}}{\mathit{b}(\mathit{\alpha}\left(1+\mathit{M}\right)\left(1\mathit{b}\left(1\mathit{s}\right)\mathit{z}+\mathit{esz}\right)+\left(\mathit{M}\left(1\mathit{s}\right)+\mathit{s}\right)\left(1+\mathit{l}+\mathit{esz}+\mathit{bz}\left(1\mathit{s}\right)\right)}\end{array}
Substituting x = x_{
e
}, y = y_{
e
} to P(x, y, z) we obtain:
\mathit{P}\left({\mathit{x}}_{\mathit{e}},{\mathit{y}}_{\mathit{e}},\mathit{z}\right)\equiv \mathit{H}\left(\mathit{z}\right)=\frac{\mathit{d}\left(\mathit{l}\left(1\mathit{b}\left(1\mathit{s}\right)\mathit{z}\right)\left(1\mathit{b}\left(1\mathit{s}\right)\mathit{z}\mathit{esz}\right)\left(1\mathit{b}\left(1\mathit{s}\mathit{\alpha}\right)\mathit{z}\right)\right)}{\mathit{b}(\mathit{\alpha}\left(1+\mathit{M}\right)\left(1\mathit{b}\left(1\mathit{s}\right)\mathit{z}+\mathit{esz}\right)+\left(\mathit{M}\left(1\mathit{s}\right)+\mathit{s}\right)\left(1+\mathit{l}+\mathit{esz}+\mathit{bz}\left(1\mathit{s}\right)\right)}.
Solving the equation H(z) = 0 we get two roots z_{e1}, z_{e2}; supposing z = z_{0} + αh_{
z
} and expanding H(z) in series by α we get the two roots up to o(α) in the form:
\begin{array}{ll}{\mathit{z}}_{\mathit{e}1}={\mathit{z}}_{01}+\mathit{\alpha}{\mathit{h}}_{1\mathit{z}},& {\mathit{z}}_{\mathit{e}2}={\mathit{z}}_{02}+\mathit{\alpha}{\mathit{h}}_{2\mathit{z}},\\ {\mathit{z}}_{01}=\frac{1}{\mathit{b}\left(1\mathit{s}\right)},& {\mathit{h}}_{1\mathit{z}}=\frac{\mathit{es}}{\mathit{b}{\left(1\mathit{s}\right)}^{2}\left(\mathit{bl}\left(1\mathit{s}\right)+\mathit{es}\right)};\\ {\mathit{z}}_{02}=\frac{1\mathit{l}}{\mathit{b}\left(1\mathit{s}\right)+\mathit{es}},& {\mathit{h}}_{2\mathit{z}}=\frac{\mathit{bl}\left(1\mathit{l}\right)}{\left(\mathit{b}\left(1\mathit{s}\right)+\mathit{es}\left)\right(\mathit{bl}\left(1\mathit{s}\right)+\mathit{es}\right)}.\end{array}
(M7)
Notice, that the root z_{e1} is positive, whereas z_{e2} is positive only for 0 < l < 1. The coordinates x = x_{
e
}(α), y = y_{
e
}(α) of B are both positive only for z = z_{e1}.
Letting x_{e1} = x_{01} + αh_{1x}, y_{e1} = y_{01} + αh_{1y}, we get
\begin{array}{l}{\mathit{x}}_{01}=\frac{\mathit{des}}{\mathit{b}\left(\mathit{M}\left(1\mathit{s}\right)\mathit{s}\right)\left(\mathit{bl}\left(1\mathit{s}\right)+\mathit{es}\right)},\\ \phantom{\rule{0.12em}{0ex}}{\mathit{h}}_{1\mathit{x}}=\frac{\mathit{des}\left({\mathit{e}}^{2}{\mathit{s}}^{2}\left(1+\mathit{M}\right)+{\mathit{b}}^{2}\mathit{el}\left(1\mathit{s}\right)\left(\mathit{M}\left(1\mathit{s}\right)\mathit{s}\right)+\mathit{bel}\left(1+2\mathit{M}\left(1\mathit{s}\right)\mathit{s}2{\mathit{s}}^{2}\right)\right)}{\mathit{b}{\left(\mathit{M}\left(1\mathit{s}\right)\mathit{s}\right)}^{2}{\left(\mathit{bl}\left(1\mathit{s}\right)+\mathit{es}\right)}^{3}};\\ {\mathit{y}}_{01}=\frac{\mathit{dl}\left(1\mathit{s}\right)}{\mathit{b}\left(\mathit{M}\left(1\mathit{s}\right)\mathit{s}\right)\left(\mathit{bl}\left(1\mathit{s}\right)+\mathit{es}\right))},\\ \phantom{\rule{0.12em}{0ex}}{\mathit{h}}_{1\mathit{y}}=\frac{\mathit{dels}\left(\mathit{es}\mathit{b}\left(1\mathit{s}\right)\left(\mathit{M}\mathit{el}\left(1+\mathit{M}\right)\left(1\mathit{s}\right)+\mathit{s}+\mathit{Ms}\right)\right)}{\mathit{b}{\left(\mathit{M}\left(1\mathit{s}\right)\mathit{s}\right)}^{2}{\left(\mathit{bl}\left(1\mathit{s}\right)+\mathit{es}\right)}^{3}}.\end{array}
(M8)
Thus, coordinates of the equilibrium B_{
e
} can be presented in the form (x_{
e
}, y_{
e
}, z_{
e
}) = (x_{1e} + αh_{1x}, y_{1e} + αh_{1y}, z_{1e} + αh_{1z}), see formulas (M7), (M8).
We apply the method of small parameter expansion to verify stability of the equilibrium B_{
e
}. Characteristic polynomial of the system can be presented in the form
\begin{array}{ll}\phantom{\rule{5.5pt}{0ex}}\mathit{E}\left(\mathit{\mu}\right)=\mathit{Det}(\left(\mathit{J}\left({\mathit{B}}_{\mathit{e}}\right)\mathit{\mu I}\right)& =\frac{\left(\mathit{d}+{\mathit{\mu}}_{0}^{2}\right)\left(\mathit{b}\left(\mathit{l}+{\mathit{\mu}}_{0}\right)\left(1\mathit{s}\right)+\mathit{es}\right)}{\mathit{b}\left(1\mathit{s}\right)}\\ \phantom{\rule{1em}{0ex}}\frac{\mathit{\alpha}}{\mathit{b}{\left(1\mathit{s}\right)}^{2}\left(\mathit{M}\left(1\mathit{s}\right)+\mathit{s}\right)\mathit{bl}\left(1\mathit{s}\right)+\mathit{es})}\mathit{Z}\left(\mathit{\mu}\right)\end{array}
where μ = μ_{0} + αh_{
μ
} and
\begin{array}{ll}\phantom{\rule{6pt}{0ex}}\mathit{z}\left(\mathit{\mu}\right)& =(({\mathit{b}}^{2}\mathit{l}(\mathit{d}(1+{\mathit{h}}_{\mathit{\mu}}(1\mathit{s}))+{\mathit{\mu}}_{0}\left({\mathit{\mu}}_{0}\left(13{\mathit{h}}_{\mathit{\mu}}\left(1\mathit{s}\right)\right)\right.\\ \phantom{\rule{1.5em}{0ex}}\left.+2\mathit{l}{\mathit{h}}_{\mathit{\mu}}\left(1\mathit{s}\right)\right)(((1\mathit{s}){}^{2}\mathit{M}(1\mathit{s})\mathit{s})\\ \phantom{\rule{1.5em}{0ex}}\mathit{e}{\mathit{s}}^{2}(\mathit{d}+{\mathit{\mu}}_{0}({\mathit{\mu}}_{0}+2{\mathit{h}}_{\mathit{\mu}}\left(1\mathit{s}\right)))((1\mathit{s}){}^{2}(\mathit{M}(1\mathit{s})\mathit{s})\\ \phantom{\rule{1.5em}{0ex}}+\mathit{be}(1\mathit{s})\mathit{s}\left({\mathit{\mu}}_{0}\right({\mathit{\mu}}_{0}(13{\mathit{h}}_{\mathit{\mu}}\left(1\mathit{s}\right))\\ \phantom{\rule{1.5em}{0ex}}4\mathit{l}{\mathit{h}}_{\mathit{\mu}}(1\mathit{s}))(\mathit{M}(1+\mathit{s})+\mathit{s})\\ \phantom{\rule{1.5em}{0ex}}+\mathit{d}\left(\left(1{\mathit{h}}_{\mathit{\mu}}\left(1\mathit{s}\right)\right)\right(\mathit{M}(1+\mathit{s})+\mathit{s})+\mathit{l}(\mathit{M}\left(1+\mathit{s}\right)\\ \phantom{\rule{1.5em}{0ex}}{\mathit{\mu}}_{0}(1+\mathit{M})(1+\mathit{s})+\mathit{s}\left)\right)\left)\right)).\end{array}
Let now {\mathit{\mu}}_{0}^{2}=\mathit{d}. Substituting this value to Z(μ) and solving equation Z(μ) = 0, we find
{\mathit{h}}_{\mathit{\mu}}=\frac{\left(\mathit{bdels}\left(\mathit{M}+{\mathit{\mu}}_{0}\left(1+\mathit{M}\right)\left(1+\mathit{s}\right)\mathit{s}\mathit{Ms}\right)\right)}{\left(2\left(\mathit{M}\left(1+\mathit{s}\right)+\mathit{s}\right)\left(\mathit{bl}\left(1+\mathit{s}\right)\mathit{es}\right)\left(\mathit{b}\left(\mathit{d}\mathit{l}{\mathit{\mu}}_{0}\right)\left(1+\mathit{s}\right)+\mathit{es}{\mathit{\mu}}_{0}\right)\right)}
and Re({\mathit{h}}_{\mathit{\mu}})=\frac{\mathit{e}\mathrm{el}\mathit{s}}{2{\left(1\mathit{s}\right)}^{3}}\left(\frac{\mathit{es}\left(\mathit{M}\left(1+\mathit{s}\right)\mathit{s}\right)+\mathit{b}\left(1\mathit{s}\right)\left(\mathit{l}\left(\mathit{M}\left(1+\mathit{s}\right)\mathit{s}\right)\mathit{d}\left(1+\mathit{M}\right)\left(1\mathit{s}\right)\right)}{{\mathit{b}}^{2}\mathit{d}\left(\mathit{M}\left(1\mathit{s}\right)\mathit{s}\right)}\frac{{\left(1\mathit{s}\right)}^{2}}{\mathit{es}+\mathit{bl}\left(1\mathit{s}\right)}\right)<0for reasonable parameter values. Thus, equilibrium B_{
e
} is asymptotically stable for s < p.
The Proposition is proven.
Threecomponent logistic model with constant immunity p
The logistic version of model (1) with p = s is reduced to 2component logistic system (M1) with respect to variables u = x + y, z. For \mathit{s}=\mathit{p}<\frac{\mathit{M}}{\mathit{M}+1} it has two equilibria, A(0, 1/a) and \mathit{B}\left(\mathit{u}=\frac{\mathit{d}}{\mathit{b}\left(\mathit{M}\left(1\mathit{s}\right)\mathit{s}\right)},\mathit{z}=\frac{\mathit{b}\left(\mathit{M}\left(1\mathit{s}\right)\mathit{s}\right)\mathit{ad}}{{\mathit{b}}^{2}\left(1\mathit{s}\right)\left(\mathit{M}\left(1\mathit{s}\right)\mathit{s}\right)}\right). According to Proposition 2, these equilibria are stable in different parameter domains.
In variables x, y, z the equilibrium B(x_{
e
}, y_{
e
}, z_{
e
}) of system (1) has coordinates
\begin{array}{l}{\mathit{x}}_{\mathit{e}}=\frac{\mathit{des}\left(\mathit{b}\left(\mathit{M}\left(1\mathit{s}\right)\mathit{s}\right)\mathit{ad}\right)}{\mathit{b}\left(\mathit{M}\left(1\mathit{s}\right)\mathit{s}\right)\left({\mathit{b}}^{2}\mathit{l}\left(1\mathit{s}\right)\left(\mathit{M}\left(1\mathit{s}\right)\mathit{s}\right)+\mathit{es}\left(\mathit{b}\left(\mathit{M}\left(1\mathit{s}\right)\mathit{s}\right)\mathit{ad}\right)\right)},\\ {\mathit{y}}_{\mathit{e}}=\frac{\mathit{bdl}\left(1\mathit{s}\right)}{\mathit{b}\left(\mathit{M}\left(1\mathit{s}\right)\mathit{s}\right)\left({\mathit{b}}^{2}\mathit{l}\left(1\mathit{s}\right)\left(\mathit{M}\left(1\mathit{s}\right)\mathit{s}\right)+\mathit{es}\left(\mathit{b}\left(\mathit{M}\left(1\mathit{s}\right)\mathit{s}\right)\mathit{ad}\right)\right)},\\ {\mathit{z}}_{\mathit{e}}=\frac{\mathit{b}\left(\mathit{M}\left(1\mathit{s}\right)\mathit{s}\right)\mathit{ad}}{{\mathit{b}}^{2}\left(1\mathit{s}\right)\left(\mathit{M}\left(1\mathit{s}\right)\mathit{s}\right)}.\end{array}
(M9)
Let p be a constant, p > s > 0 . According to Statement 1, the system has trivial equilibrium O(x = 0, y = 0, z = 0), which is unstable for all parameter values, and the equilibrium \phantom{\rule{0.24em}{0ex}}\mathit{A}\left(\mathit{x}=0,\mathit{y}=\frac{1}{\mathit{a}},\mathit{z}=0\right), which is stable if \mathit{M}<\frac{\mathit{ad}+\mathit{bs}}{\mathit{b}\left(1\mathit{s}\right)} and unstable if \mathit{M}>\frac{\mathit{ad}+\mathit{bs}}{\mathit{b}\left(1\mathit{s}\right)}. The system has also a nontrivial equilibrium B whose x, y  coordinates are expressed via z coordinate:
\begin{array}{l}\mathit{x}=\frac{\mathit{d}+\mathit{b}{\mathit{f}}^{\pm}\left(\mathit{z}\right)\left(\mathit{M}\left(1\mathit{s}\right)\mathit{s}\right)}{\mathit{b}\left(1+\mathit{M}\right)\left(\mathit{p}\mathit{s}\right)},\\ \mathit{y}=\frac{\mathit{d}\mathit{b}{\mathit{f}}^{\pm}\left(\mathit{z}\right)\left(\mathit{M}\left(1\mathit{p}\right)\mathit{p}\right)}{\mathit{b}\left(1+\mathit{M}\right)\left(\mathit{p}\mathit{s}\right)},\text{where}\\ {\mathit{f}}^{\pm}\left(\mathit{z}\right)=\frac{1+\mathit{M}\mathit{bz}\pm \sqrt{{\left(1+\mathit{M}\mathit{bz}\right)}^{2}4\mathit{ad}\left(1+\mathit{M}\right)\mathit{z}}}{2\mathit{a}\left(1+\mathit{M}\right)}\end{array}
(M10)
and z coordinate solves the equation:
\begin{array}{l}\left(1+\mathit{l}+\mathit{a}{\mathit{f}}^{\pm}\right)(\mathit{d}\mathit{b}{\mathit{f}}^{\pm}(\mathit{M}(1\mathit{s})\mathit{s}))+(\mathit{des}{\mathit{b}}^{2}{\mathit{f}}^{\pm}(1\mathit{p})(\mathit{M}(1\mathit{s})\mathit{s})\\ \phantom{\rule{1.78em}{0ex}}+\mathit{b}(\mathit{d}\left(1\mathit{p}\right)\mathit{e}{\mathit{f}}^{\pm}(\mathit{M}(1\mathit{p})\mathit{p})\mathit{s}))\mathit{z}=0.\end{array}
(M11)
Denote
\begin{array}{l}{\mathit{h}}^{\pm}\left(\mathit{z}\right)\equiv \frac{\left(\mathit{d}\left(\mathit{l}1+\mathit{a}{\mathit{f}}^{\pm}\right)+\mathit{b}{\mathit{f}}^{\pm}\left(\mathit{M}\left(1\mathit{a}{\mathit{f}}^{\pm}\mathit{l}\left(1\mathit{s}\right)\right)+\mathit{ls}\right)\right)+\left(\mathit{d}\mathit{b}{\mathit{f}}^{\pm}\mathit{M}\right)\left(\mathit{b}\left(1\mathit{s}\right)+\mathit{es}\right)\mathit{z}}{\mathit{b}{\mathit{f}}^{\pm}\left(1+\mathit{M}\right)\left(1\mathit{a}{\mathit{f}}^{\pm}\left(\mathit{b}\left(1\mathit{s}\right)+\mathit{es}\right)\mathit{z}\right)}\end{array}
then \overline{\mathit{z}} is a root of the equations p = h^{±}(z). Solutions of the latter equation are the points of intersection of the line 0 < p ≤ 1 and the curves h^{±}(z). Two cases with different values of the parameter M are presented in Additional file 2. It demonstrates that at most two positive values of \overline{\mathit{z}} are possible and hence the system may have at most two positive equilibria, whose x, y ‒ coordinates are expressed via \overline{\mathit{z}} by the equations (M9).
Let us define the critical value of the parameter M, {\mathit{M}}^{\mathit{tc}}=\frac{\mathit{ad}+\mathit{bs}}{\mathit{b}\left(1\mathit{s}\right)}.
Proposition 4. For \mathit{s}\le \mathit{p}<\frac{\mathit{M}}{\mathit{M}+1}and fixed positive values of d, b < 1, l, e system (1) with a = 1 has a single stable equilibrium \phantom{\rule{0.24em}{0ex}}\mathit{A}\left(\mathit{x}=0,\mathit{y}=\frac{1}{\mathit{a}},\mathit{z}=0\right)if M < M^{tc}and a single stable equilibrium B(x_{0}, y_{0}, z_{0}) if M > M^{tc}; here the coordinates (x_{0}, y_{0}, z_{0}) are defined by formulas (M10, M11) if s < p and by (M9) if s = p. Transcritical bifurcation "changing of stability" of the points A and B happens as M = M^{tc}.
For p = s system (1) has nontrivial point B, whose coordinates are expressed by (M9); it is easily to verify that {\mathit{u}}_{\mathit{e}}={\mathit{x}}_{\mathit{e}}+{\mathit{y}}_{\mathit{e}}=\frac{\mathit{d}}{\mathit{b}\left(\mathit{M}\left(1\mathit{s}\right)\mathit{s}\right)},{\mathit{z}}_{\mathit{e}}=\frac{\mathit{b}\left(\mathit{M}\left(1\mathit{s}\right)\mathit{s}\right)\mathit{ad}}{{\mathit{b}}^{2}\left(1\mathit{s}\right)\left(\mathit{M}\left(1\mathit{s}\right)\mathit{s}\right)}, and 1  a(x_{
e
} + y_{
e
})  b(1  s)z_{
e
} = 0. Hence, z_{
e
} > 0 if > M^{tc}; at this condition the equilibrium A looses stability according to Statement 1.
Characteristic equation at the equilibrium B for p = s is of the form:
\phantom{\rule{0.25em}{0ex}}\mathit{Det}(\left(\mathit{J}\left(\mathit{B}\right)\mathit{\mu I}\right)=\left(1+\mathit{l}+\mathit{\mu}+\mathit{a}\left({\mathit{x}}_{\mathit{e}}+{\mathit{y}}_{e}\right)+\mathit{b}\left(1\mathit{s}\right){\mathit{z}}_{e}+{\mathit{sz}}_{e}\right)*
\left({\mathit{\mu}}^{2}+{\mathit{b}}^{2}\left(1\mathit{s}\right)\left(\mathit{M}\left(1\mathit{s}\right)\mathit{s}\right)\left({\mathit{x}}_{\mathit{e}}+{\mathit{y}}_{\mathit{e}}\right){\mathit{z}}_{\mathit{e}}+\mathit{\mu}\left(1+2\mathit{a}\left({\mathit{x}}_{\mathit{e}}+{\mathit{y}}_{\mathit{e}}\right)+{\mathit{bz}}_{\mathit{e}}\left(1\mathit{s}\right)\right)\right)=0.\phantom{\rule{0.25em}{0ex}}
Accounting the equality 1  a(x_{
e
} + y_{
e
})  b(1  s)z_{
e
} = 0, we obtain from the first term of the last equation: μ_{1} = 1  l  au_{
e
}  b(1  s)z_{
e
}  esz_{
e
} =  l  sz_{
e
} < 0, and from the second term {\mathit{\mu}}_{2,3}=\frac{1}{2}\left(\mathit{au}\pm \sqrt{4\left(\mathit{b}\left(\mathit{M}\left(1\mathit{s}\right)\mathit{s}\right)\mathit{ad}\right){\mathit{u}}_{\mathit{e}}+{\left(\mathit{au}\right)}^{2}}\right). For = M^{tc}μ_{2} = 0, μ_{3} =  au < 0. Hence, the transcritical (pitchfork) bifurcation happens with points A, B [41], ch.7. When M > M^{tc} the eigenvalues μ_{2,3} have negative real parts. So, the equilibrium point B is stable for p = s and M > M^{tc}. Due to continuity arguments the point B is also stable for p > s if p is close to s. A general case p > s was verified by computation of eigenvalues of complete characteristic equation in B.
Threecomponent Malthusian model with the immunity p = p(z); proof of Statement 2
The coordinates of nontrivial equilibria are given by formulas (M5) and (M6). Solving equation (M5) with respect to p = p(z), we obtain the equation for zcoordinate:
\begin{array}{ll}\phantom{\rule{5pt}{0ex}}\mathit{p}\left(\mathit{z}\right)\equiv \left(1\mathit{s}\right){\mathit{e}}^{\mathit{kz}}+\mathit{s}& =\frac{\left(1\mathit{l}\right)+\left(\mathit{b}\left(2\mathit{l}\left(1\mathit{s}\right)\mathit{s}\right)+\mathit{es}\right)\mathit{z}\mathit{b}\left(\mathit{b}\left(1\mathit{s}\right)+\mathit{es}\right){\mathit{z}}^{2}}{\mathit{bz}\left(1\left(\mathit{b}\left(1\mathit{s}\right)+\mathit{es}\right)\mathit{z}\right)}\\ \phantom{\rule{1em}{0ex}}\equiv \mathit{h}\left(\mathit{z}\right)\end{array}
(M12)
Evidently, the zcoordinate of possible equilibrium does not depend on the parameter M. Next, h(z) → 1, p(z) → s < 1 as z → ∞, so that two general cases of mutual placing of p(z) and h(z) for l > 1 and l < 1 are possible, which are shown in Additional file 3. Equation (M12) has up to two ("small") positive roots z = z_{
e
} if l < 1 (see Additional file 3a) and one ("large") positive root z_{
e
} if l > 1 (see Additional file 3b). Equilibrium values x_{
e
}(z), y_{
e
}(z) defined by the formulas (M6) have different signs for small values of z_{
e
} and are both positive for large z_{
e
}. Accordingly, the system can have only one positive equilibrium. Eigenvalues of this equilibrium can be computed from the equation Det(J(x_{
e
}, y_{
e
}, z_{
e
})  μI) = 0. Using the LOCBIF software [55], we show that one eigenvalue is real and negative but two other are complex with a positive real part. Thus, this equilibrium is unstable. Statement 2 is proven.
Threecomponent logistic model with the immunity p = p(z)
Proof of Statement 3
Let B_{
M
}(x, y, z) be a nontrivial stable equilibrium of model (1), (3) whose coordinates (x(M), y(M), z(M)) depend on the parameter M and satisfy system (2); let P(μ) ≡ Det(J(B)  μI) = 0 be the characteristic polynomial of the system around B_{
M
}. The supercritical Hopf bifurcation, corresponding to changing of stability of B_{
M
} accompanied by appearance a stable limit cycle happens in the system when for some M a pair of eigenvalues becomes imaginary and certain conditions of nondegeneracy are fulfilled (see, for example, [41]). The Hopf bifurcation is supercritical if the first Lyapunov value becomes negative. The existence of this bifurcation in logistic version of model (1), (3) was verified using LOCBIF [55]. The program numerically finds coordinates of equilibrium with imaginary eigenvalues under variation of parameter M and one more parameter of the model (e.g., e, l, or s) for fixed values of other parameters, checks the sign of the first Lyapunov value and verifies the conditions of nondegeneracy formulated in the Hopf theorem. Applying this software we have found the parameter curves of Hopf bifurcation e_{
H
}(M), l_{
H
}(M); s_{
H
}(M) for fixed values of other parameters of the model (see Additional file 4); it proves the assertions of the Statement.