Analysis of yellow fever prevention strategy from the perspective of mathematical model and cost-effectiveness analysis

1.   Introduction

The traditional predator–prey model was built to explain the survival of marine fish, which has been extensively researched and improved to more accurately describe the development of populations. However, there is no upper limit to the growth rate of predators, which is not suitable for most populations in nature. Based on this fact, Leslie ^[1,2] considered that the growth rates of both predator and prey have upper limits by assuming that the "carrying capacity" of the predator's environment is proportional to the number of prey and proposed the following well-known Leslie–Gower model:

where $\frac{y}{nx}$ denotes the Leslie–Gower term and $nx$ is the "carrying capacity" of the predator. Korobeinikov ^[3] successively verified that the fixed point of the Leslie–Gower model is globally stable. System (1.1) was modified by introducing different types of functional response, and the stability and bifurcations were analyzed ^{[4,5,6,7,8,9,10]}.

In the above system, the average growth rate $\frac{\dot{x}}{x}$ of the traditional logistic model is a linear function of the prey. However, in 1963, Smith ^[11] found that the average growth rate of large daphnia populations does not follow linear functions and demonstrated that food is needed for the maintenance and growth of the population during the growth stage but is only required for maintenance when the population reaches saturation. Then proposed the "limited food" model, sometimes called the "limited resources" model. When introducing Smith growth into system (1.1), the system takes the following form:

where $rx\big(\frac{K-x}{K+ax}\big)$ is the Smith growth function. With the influence of environmental toxicants, Kumar ^[12] discussed the spatiotemporal dynamics and delay-induced instability of a Leslie-Gower prey-predator model with Smith growth rate. Li and He ^[13] studied the Hopf bifurcation of a delayed food-limited model with feedback control; Feng et al. ^[14] considered the stability and Hopf bifurcation of a modified Leslie-Gower predator–prey model with Smith growth rate and Beddington-DeAngelis functional response; Bai et al. ^[15] proposed a predator–prey model with a Smith growth function and the additive predation in prey and studied the stability and Hopf bifurcation of the system. For more results on dynamical behaviors of systems with Smith growth function, please see^[16,17].

On the other hand, sustainable harvesting plays an important role in maintaining ecological balance, and its influence on population development has been extensively studied. Huang et al. ^[18] investigated the bifurcation of predator-prey models with constant-yield prey harvesting. Xu et al. ^[19] discussed the degenerate Bogdanov-Takens bifurcation of a Holling-Tanner predator-prey model with constant-yield prey harvesting. Wu et al. ^[20] studied the bifurcations such as the degenerate Bogdanov-Takens bifurcation and degenerate Hopf bifurcation of a Holling-Tanner model with generalist predator and constant-yield prey harvesting. Zhu and Lan ^[21] modified system (1.1) by introducing constant-yield prey harvesting

where $H > 0$ is the constant-yield prey harvesting. The authors discussed the saddle node bifurcation, supercritical, and subcritical Hopf bifurcations of the system.

Motivated by systems (1.2) and (1.3), in this paper we propose the following system:

where $x(t)$ and $y(t)$ are the densities of the prey and predator at time $t$, respectively; $r > 0$ and $\delta > 0$ are the intrinsic growth rates of the prey and predator, respectively; $K > 0$ is the carrying capacity of the prey; $\beta > 0$ is the maximum rate of predation; $H > 0$ is the constant-yield prey harvesting; $\frac{r}{a}$ with $a > 0$ is the mass substitution rate of the population at $K$; $n > 0$ is the quality of prey provided to the predator.

For simplicity, we make the following transformations:

and drop the bars, then system (1.4) can be rewritten as

where

We can easily verify that all the solutions of system (1.5) are positive and bounded with positive initial conditions.

The rest of the paper is organized as follows: In Section 2, we will investigate the existence and stability of the boundary and positive equilibria of system (1.5). In Section 3, we will prove the Bogdanov–Takens bifurcation of codimension 2 and the degenerate Hopf bifurcation of codimension 3. In Section 4, we will present some bifurcation diagrams and phase portraits to support our theoretical results. We give a brief discussion in the last section.

2.   Existence and stability of equilibria

In this section, we discuss types of the boundary and positive equilibria of system (1.5) in a positive invariant and bounded region

For any nonnegative equilibrium $E(x, y)$, the Jacobian matrix of system (1.5) is

and

2.1.   Existence and stability of boundary equilibria

Substituting $y = 0$ into the first equation of (1.5), we have

whose discriminant is

Denote

Then we have the following theorem.

Theorem 2.1. For the boundary equilibria of system (1.5), we have:

$\rm(1)$ If $h < \frac{1}{4}$ and $0 < a < \frac{1-2\sqrt{h}}{h}$, then system (1.5) has two boundary equilibria: an unstable node $E_{10}(x_{10}, 0)$ and a saddle $E_{20}(x_{20}, 0)$;

$\rm(2)$ if $h < \frac{1}{4}$ and $a = \frac{1-2\sqrt{h}}{h}$, then system (1.5) has a unique boundary equilibrium $\bar{E}_{*}(\bar{x}_{*}, 0)$, which is a saddle node, including an unstable parabolic sector in the left.

Proof. If $ah < 1$ and $\Delta_{f_0} > 0$, that is $h < \frac{1}{4}$ and $0 < a < \frac{1-2\sqrt{h}}{h}$, $f_0(x)$ has two positive roots in the interval $(0, 1)$, which implies that system (1.5) has two boundary equilibria $E_{10}(x_{10}, 0)$ and $E_{20}(x_{20}, 0)$.

If $ah < 1$ and $\Delta_{f_0} = 0$, that is $h < \frac{1}{4}$ and $a = \frac{1-2\sqrt{h}}{h}$, $f_0(x)$ has a unique positive root in the interval $(0, 1)$, which implies that system (1.5) has a unique boundary equilibrium $\bar{E}_{*}(\bar{x}_{*}, 0)$.

The Jacobian matrix of system (1.5) at $E_0(x, 0)$ is

We can easily verify that

Hence, $E_{10}$ is an unstable node and $E_{20}$ is a saddle.

When $a = \frac{1-2\sqrt{h}}{h}$, make the following transformations successively

then system (1.5) becomes

By Theorem 7.1 in Chapter 2 of Zhang et al. ^[22], $\bar{E}_{*}$ is a saddle node, which includes an unstable parabolic sector in the left if $h < \frac{1}{4}$ and $a = \frac{1-2\sqrt{h}}{h}$. The proof is completed.

Choose $(h, b, s) = (0.08, 1.2, 0.1)$, and we present the relationship between the boundary equilibria of system (1.5) and $a$, as well as the phase portrait when $a = 0.1$, shown in Figure 1. Figure 1(a) shows that when $0 < a < 5.428932188$, system (1.5) has two boundary equilibria; when $a\doteq5.428932188$, system (1.5) has a unique boundary equilibrium; and when $a > 5.428932188$, system (1.5) has no boundary equilibrium. Figure 1(b) shows that system (1.5) admits two boundary equilibria, where $E_{10}$ is an unstable node and $E_{20}$ is a saddle. Further, when the initial values lie inside the region to the right of the stable manifold of the saddle (the positive equilibrium), both populations will coexist in fixed sizes; otherwise, both populations will not coexist.

2.2.   Existence and stability of positive equilibria

Note that any positive equilibrium of system (1.5) satisfies

From Eq (2.2), we denote

then

whose discriminant

Obviously, in the interval $(0, 1)$, $f'(x) = 0$ has a unique positive root

For any positive equilibrium $E(x, y)$, from $f(x) = 0$, we obtain

Substituting this into ${\rm Det}(J_{E})$ and $f'(x)$, there is

Theorem 2.2. For the positive equilibria of system (1.5), the following conclusions hold:

$\rm(1)$ If $ah < 1$ and $f(x_*) < 0$, then system (1.5) has two positive equilibria: a hyperbolic saddle $E_{11}(x_{11}, x_{11})$, and an elementary and antisaddle equilibrium $E_{12}(x_{12}, x_{12})$, with $x_{11} < x_* < x_{12}$;

$\rm(2)$ if $ah < 1$ and $f(x_*) = 0$, then system (1.5) has a unique positive degenerate equilibrium $E_{*}(x_{*}, x_{*})$;

$\rm(3)$ if $ah\geq1$ or $f(x_*) > 0$, then system (1.5) has no positive equilibrium.

Proof. From (2.3) and (2.4), it is obvious that ${\rm Det}(J_{E_{11}}) < 0$, ${\rm Det}(J_{E_{12}}) > 0$ and ${\rm Det}(J_{E_{*}}) = 0$. Thus $E_{11}$ is a hyperbolic saddle, $E_{12}$ is an elementary and antisaddle equilibrium, and $E_{*}$ is a degenerate equilibrium. The proof is completed.

Next, we consider case (2) of Theorem 2.2. From $f(x_*) = f'(x_*) = 0$, we can express

furthermore, we let

Theorem 2.3. Assume that $b = b_*$, $h = h_*$ and $a < \frac{1-2x_*}{x_*^2}$ with $x_* < \frac{1}{2}$, then system (1.5) admits a degenerate positive equilibrium $E_*(x_*, x_*)$. Moreover,

$\rm(1)$ If $s > s_* \; ($or $0 < s < s_*)$, then $E_*$ is a saddle node, which includes a stable $($or an unstable$)$ parabolic sector in the left;

$\rm(2)$ If $s = s_*$, then $E_*$ is a cusp of codimension $2$.

Proof. (1) When $s\neq s_*$, make the following transformation

then system (1.5) can be written as

In order to translate the linear part of this system to Jordan form, we make the following transformation:

Then system (2.5) becomes

where

with

Note that when $a < \frac{1-2x_*}{x_*^2}$ with $x_* < \frac{1}{2}$, $A_1 < 0$, $A_2 < 0$ and $A_3 > 0$, thus $a_{20} < 0$. By Theorem 7.1 in Chapter 2 of Zhang et al. ^[22], $E_*$ is a saddle node, which includes a stable (or an unstable) parabolic sector in the left if $s > s_* \; ($ or $0 < s < s_*)$.

(2) When $s = s_*$, by the following transformations successively

system (1.5) can be rewritten as

By Lemma 3.1 of ^[18], system (2.7) near the origin is equivalent to

where

It is obvious that when $a < \frac{1-2x_*}{x_*^2}$ with $x_* < \frac{1}{2}$, we have $D > 0$ and $\check{E} < 0$, thus $E_*$ is a cusp of codimension 2. The proof is completed.

3.   Bifurcations

In this section, we investigate the Bogdanov-Takens bifurcation of codimension 2 and the degenerate Hopf bifurcation of codimension 3.

3.1.   Bogdanov-Takens bifurcation of codimension 2

According to Theorem 2.3 (2), system (1.5) may have a Bogdanov-Takens bifurcation of codimension 2 around $E_*$. Now we choose $a$ and $h$ as bifurcation parameters and obtain the following unfolding system

where $\lambda = (\lambda_1, \lambda_2)$ is a parameter vector in a small neighborhood of $(0, 0)$.

Theorem 3.1. If the conditions of Theorem $2.3 \; (2)$ hold, then system (1.5) undergoes a Bogdanov–Takens bifurcation of codimension $2$ around $E_*$.

Proof. In order to move $E_{*}$ to the origin, we make a transformation

then system (3.1) becomes

where $P_1(x_1, y_1, \lambda)$ and $P_2(x_1, y_1, \lambda)$ are $C^{\infty}$ functions of at least third order with respect to $(x_1, y_1)$, and

By the following transformation

system (3.2) can be rewritten as

where $P_3(x_2, y_2, \lambda)$ is a $C^{\infty}$ function of at least third order with respect to $(x_2, y_2)$, and

Make a transformation

still denoting $\tau$ as $t$, then system (3.3) becomes

where $P_4(x_3, y_3, \lambda)$ is a $C^{\infty}$ function of at least third order with respect to $(x_3, y_3)$, and

Note that when $\lambda_1$ and $\lambda_2$ are small enough,

where $A_1$ and $A_2$ are given in (2.6).

Let

still denoting $\tau$ as $t$, then we can rewrite system (3.4) as

where $P_5(x_4, y_4, \lambda)$ is a $C^{\infty}$ function of at least third order with respect to $(x_4, y_4)$.

It is obvious that when $\lambda_1$ and $\lambda_2$ are small enough,

Making the following transformations successively

we obtain the following versal unfolding of system (3.1)

where $P_6(x_6, y_6, \lambda)$ is a $C^{\infty}$ function of at least third order with respect to $(x_6, y_6)$, and

By Maple software, we can obtain

By the results of ^[23], system (3.1) is the versal unfolding of Bogdanov–Takens sigularity (cusp case) of codimension 2, when $\lambda = (\lambda_1, \lambda_2)$ is a parameter vector in a small neighborhood of $(0, 0)$. Hence, system (1.5) undergoes a Bogdanov-Takens bifurcation of codimension 2. The proof is completed.

3.2.   Degenerate Hopf bifurcation of codimension 3

The proof of Theorem 2.2 shows that ${\rm Det}(J_{E_{12}}) > 0$, that is, system (1.5) may undergo a Hopf bifurcation around $E_{12}$ when ${\rm Tr }(J_{E_{12}}) = 0$. For simplicity, we denote $x_{12}$ by $z$. From $f(z) = 0$ and ${\rm Tr }(J_{E_{z}}) = 0$, we obtain

By $b > 0$ and ${\rm Det}(J_{E_{z}}) > 0$, we denote

To obtain the focal values around $E_{12}(z, z+k)$, we make the following transformations successively

where

Still denoting $\tau$ by $t$, system (1.5) can be rewritten as

where

The first three Lyapunov coefficients ^[22] at $E_{12}$, respectively, are

where

and $f_2$, $f_3$ are polynomials with respect to $a, s$ and $z$, whose expressions are omitted.

For the polynomials $f, g$ and $h_i \; (i = 1, 2, \cdots, n)$, let $V(h_1, h_2, \cdots, h_n)$ be the set of common zeros of $h_i \; (i = 1, 2, \cdots, n)$, ${\rm res}(f, g, x)$ be the Sylvester resultant of $f$ and $g$ with respect to $x$, and ${\rm lcoeff}(f, x)$ be the leading coefficient of $f$ with respect to $x$. By Maple software, we can obtain

where

$g_i\ (1 = 1, 2, 3)$ are polynomials with respect to $a$ and $z$, $h_i\ (1 = 1, 2, 3)$ are polynomials with respect to $z$, and their expressions are omitted.

Let

Theorem 3.2. Assume that (3.5) and (3.6) hold, then $E_z$ is a weak focus of order at most $3$.

Proof. Note that $B_3 < 0$ in $\mathcal{M}$, thus

We can easily verify that $P_1 < 0$, $P_3 < 0$, $P_4 < 0$ and $25z^2+20z-94 < 0$ in $\mathcal{M}$. Then by Lemma 2 in Chen and Zhang ^[24], we obtain the following decomposition

where

Next, we prove $V(L_1, L_2, L_3)\cap \mathcal{M} = \emptyset$ by three steps.

Step 1. Proving that $V_1\cap \mathcal{M} = \emptyset$. We can easily verify that $P_2 < 0$ when $\frac{1}{4} \leq z < \frac{1}{2}$, and when $0 < z < \frac{1}{4}$, in $\mathcal{M}$, $P_2 = 0$ has a unique positive root

Note that $5-a_1z > 0$ when $0 < z < \frac{1}{4}$, which leads to

Thus, $V_1\cap \mathcal{M} = \emptyset$.

Step 2. Proving that $V_2\cap \mathcal{M} = \emptyset$. Note that

By the command "realroot$(h_1, 10^{-30})$", $h_1(z)$ exists only one positive real root isolation interval $[\underline{z}, \overline{z}]\subset(0, \frac{1}{2})$, where

For the real root interval $[\underline{z}, \overline{z}]$ of $h_1(z)$, using the real root isolation algorithm of multivariate polynomial systems ^[25], we can obtain a unique positive real root isolation interval $[\underline{a}, \overline{a}]$ in $\mathcal{M}$, where

Note that

This means $f_1(s, a, z) > 0$ in the real root isolation interval $[\underline{z}, \overline{z}]\times[\underline{a}, \overline{a}]$ of $\{h_1, g_1\}$. Thus $V_2\subseteq V(f_1, g_1, h_1) = \emptyset$, which implies $V_2\cap \mathcal{M} = \emptyset$.

Step 3. Proving that $V_3\cap \mathcal{M} = \emptyset$. By Maple software, there is

which implies $V(h_2, h_3) = \emptyset$. Therefore, $V_3\cap \mathcal{M} = \emptyset$.

To sum up, we show that $V(f_1, f_2, f_3)\cap \mathcal{M} = \emptyset$. Thus $E_z$ is a weak focus of order at most $3$. The proof is completed.

Define

From Theorem 3.2, we have the following theorem.

Theorem 3.3. Assume that (3.5) and (3.6) hold, then $E_z$ is a weak focus of order up to $3$. Moreover,

$(1)$ if $(a, s, z)\in\mathcal{M}_1$, then $E_z$ is a weak focus of order $1$;

$(2)$ if $(a, s, z)\in\mathcal{M}_2$, then $E_z$ is a weak focus of order $2$;

$(3)$ if $(a, s, z)\in\mathcal{M}_3$, then $E_z$ is a weak focus of order $3$.

To verify $\mathcal{M}_3\neq\emptyset$, we let $z = \frac{1}{20}$. Again from $f(z) = 0$ and ${\rm Tr }(J_{E_{z}}) = 0$, we have

By $b > 0, \ h > 0$ and ${\rm Det}(J_{E_{z}}) > 0$, we define

We can obtain the following theorem.

Theorem 3.4. When $z = \frac{1}{20}$, assume that (3.7) and (3.8) hold, and there exist $a_0\in[\underline{a}_0, \overline{a}_0]$ and $s_0\in[\underline{s}, \overline{s}]$ given respectively in (3.9) and (3.10), then $E_z$ is a weak focus of order $3$ and system (1.5) can undergo a degenerate Hopf bifurcation of codimension $3$ near $E_z$.

Proof. When $z = \frac{1}{20}$, by the command "realroot$(g_1, 10^{-30})$", $g_1(a, \frac{1}{2})$ exists a positive real root isolation interval $[\underline{a}_0, \overline{a}_0]$ in $\mathcal{M}_*$, where

For the real root interval $[\underline{a}_0, \overline{a}_0]$ of $g_1(a, \frac{1}{2})$, using the real root isolation algorithm of multivariate polynomial systems ^[25], we can obtain a unique positive real root isolation interval $[\underline{s}, \overline{s}]$ in $\mathcal{M}_*$, where

Hence, under the conditions (3.7) and (3.8), when $z = \frac{1}{20}$, there exist $a_0\in[\underline{a}_0, \overline{a}_0]$ and $s_0\in[\underline{s}, \overline{s}]$ such that $f_1(a_0, s_0, \frac{1}{20}) = f_2(a_0, s_0, \frac{1}{20}) = 0$. Therefore, according to Theorem 3.2, $E_z$ is a weak focus of order exactly $3$.

Finally,

where the expression of $J(a, s)$ is too long and is omitted here.

For the polynomial $f(x, y)$, let $f^+(x, y)$ and $f^-(x, y)$ be the summations of the positive and negative terms in $f(x, y)$, respectively (^[25]).

According to Theorem 2.2 in ^[25], we have

Note that $a^2s+40as+a+400s-360\neq0$, thus $\left|\frac{\partial ({\rm Tr}(J_{E_{z}}), L_1, L_2)}{\partial (a, s, z)}\right|\neq0$, for $(a, s, z) = (a_0, s_0, \frac{1}{20})$. Therefore, system (1.5) can undergo a degenerate Hopf bifurcation of codimension $3$ near $E_z$. The proof is completed.

4.   Numerical simulations

In this section, we give some numerical simulations by Matcont to support the bifurcation phenomena of system (1.5).

Choose $a = 0.142$, $b = 14.2$, and $s = 0.1$, and present the bifurcation diagram in the $(h, x)$-plane by Matcont. Figure 2 shows that when $0 < h < 0.0088222644$, the positive equilibrium $E_{12}$ is stable and system (1.5) has no limit cycle; when $0.0088222644 < h < 0.0088313334$, system (1.5) admits an unstable limit cycle around $E_{12}$; when $0.0088313334 < h < 0.0088955378$, system (1.5) has two limit cycles, where the inner one is stable and the outer one is unstable; and when $h > 0.0088955378$, the positive equilibrium $E_{12}$ is unstable and system (1.5) has no limit cycle. It is shown that when the constant-yield prey harvesting $h$ is relatively small, the constant-yield prey harvesting does not change the dynamic behaviors of the system, that is, both populations will coexist in fixed sizes $E_{12}$; as $h$ increases, system (1.5) first generates an unstable limit cycle and then two limit cycles (the inner one is stable); but when $h$ is relatively large, the limit cycles will disappear first and then the positive equilibria, and finally there will be no positive equilibrium point.

In Figure 3, we present the corresponding phase portraits of system (1.5) with $a = 0.142$, $b = 14.2$, and $s = 0.1$, as $h$ varies. When $h = 0.008$, from Figure 3(a), the positive equilibrium $E_{12}$ is stable, that is, when the initial values (i.e. initial densities of both populations) lie inside the region to the right of the stable manifold of the saddle $E_{11}$, both populations will coexist in fixed sizes $E_{12}$; when $h = 0.00883$, from Figure 3(b), system (1.5) admits a unstable limit cycle, that is, when the initial values lie inside the unstable limit cycle, both populations will coexist in fixed sizes $E_{12}$; when $h = 0.00888$, from Figure 3(c), system (1.5) admits two limit cycles (the inner one is stable), that is, when the initial values lie inside the outer limit cycle, both populations will coexist in the stable limit cycle; otherwise, both populations will not coexist in the above three figures. When $h = 0.0098$, from Figure 3(d), all the nonnegative equilibria are unstable, that is, both populations will not coexist. Hence the relatively small harvesting can make both populations coexist in fixed sizes or in a periodic orbit; that is, the small harvesting can affect the populations of both species to some extent; however, the large harvesting is detrimental to the dynamic behaviors of the system.

5.   Conclusions

In this paper, we discussed the stability and bifurcation of a Leslie–Gower predator–prey model with Smith growth and constant-yield harvesting. We proved that system (1.5) admits at most two boundary equilibria, both of which are unstable. For system (1.5) without Smith growth, that is system (1.3), Zhu and Lan ^[21] investigated the stability of equilibria. Also, system (1.3) undergoes the saddle-node bifurcation and supercritical and subcritical Hopf bifurcations. Compared with the results in ^[21], we presented some different dynamics of the system, and found that both the Smith growth and constant-yield harvesting play important roles in the dynamics of system (1.5), more precisely, we proved that the degenerate positive equilibrium of system (1.5) is a cusp of codimension 2 and system (1.5) undergoes cusp-type Bogdanov–Takens bifurcation of codimension 2. Further, using the resultant elimination method and the real root isolation algorithm of multivariate polynomial systems ^[25], we proved that system (1.5) has a weak focus of order at most 3 and undergoes a degenerate Hopf bifurcation of codimension 3. By numerical simulation, we illustrated that system (1.5) has two limit cycles (the inner one is stable). Also, we showed that the small harvesting is conducive to the survival of prey and predators, but the large harvesting is detrimental to the coexistence of the populations. Therefore, our results revealed that Smith growth and constant-yield harvesting can lead to richer dynamic behaviors.

Motivated by Lampart and Zapom$\check{e}$l ^[26], whether there are some new numerical methods to derive some new dynamical behaviors. Also, we can regard the presence of more complex phenomena such as intermittency (see ^[27]) of the system in the future.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgements

This work was supported by the Natural Science Foundation of Fujian Province (2021J01613, 2021J011032, 2023J05250) and the Scientific Research Foundation of Minjiang University (MJY22027).

Conflict of interest

The authors declare there is no conflict of interest.