Modeling effects of matrix heterogeneity on population persistence at the patch-level

1.   Introduction

Seasonal influenza is an acute respiratory infection caused by influenza viruses. Worldwide, these annual epidemics are estimated to result in about 3 to 5 million cases of severe illness, and about 290,000 to 650,000 respiratory deaths ^[1]. This infection can have an endemic, epidemic or pandemic behavior.

There were, three major flu pandemics during the 20th century, the so-called Spanish flu (H1N1) in 1918 was the most devastating pandemic. It has been estimated that the Spanish flu claimed around 40–50 million deaths (as much as 3% of the total population), and it also infected 20–40% of the whole population. In 1957–1958, the Asian flu or bird flu pandemic (H2N2) caused more than two million deaths ^[2]. Unlike the Spanish flu, this time the infection-causing virus was detected earlier due to the advancement of science and technology. A vaccine was made available but with limited supply. After a decade (in 1968), a flu pandemic (H3N2) that originated again in Hong Kong hit mankind. That flu pandemic also claimed one million lives. In 2009, the H1N1 swine flu is one of the more publicized pandemics that attracted the attention of all scientists and health professionals in the world and made them very much concerned. However, the pandemic did not result in great casualties like before. As of July 2010, only about 18,000 related deaths had been reported ^[2]. Besides the 4 influenza pandemics since 1918, annual seasonal influenza epidemics have spread among nations on smaller scales. There are many methods of preventing the spread of infectious disease, one of them is vaccination. Vaccination is the administration of agent-specific, but relatively harmless, antigenic components that in vaccinated individuals can induce protective immunity against the corresponding infectious agent ^[3].

Influenza causes serious public-health problems around the world, therefore, we need to understand transmission mechanisms and control strategies. Mathematical models also provided insight into the severity of past influenza epidemics. Some models were used to investigate the three most devastating historical pandemics of influenza in the 20th century ^[4,5,6]. There are a lot of pathogens with several circulating strains.

An important factor when analyzing the dynamics of a disease is the way in which it is transmitted from an infected individual to a healthy one. The incidence rate of a disease is defined as the number of susceptible individuals that become infected per unit of time. It measures the number of new cases of a disease in a period of time. There are different types of incidence functions that have been used in literature in order to model the force of infection of a disease. For example, Rahman and Zou ^[2] used the bilinear incidence rate $\beta SI$. However, there are more realistic incidence rates than the bilinear incidence rate, For instance, Capasso and his co-workers observed in the seventies ^[7] that the incidence rate may increase more slowly as $I$ increases, so they proposed a saturated incidence rate $\frac{\beta I S}{1+\zeta I}$.

Baba and Hincal ^[8] studied an epidemic model consisting of three strains of influenza ($I_1$, $I_2$, and $I_3$) where we have vaccine for strain 1 ($V_1$) only, and force of infection $\frac{\beta SI} {1 + \zeta S}$ for strain 2. Baba et al. ^[9] studied an epidemic model consisting of two strains of influenza ($I_1$ and $I_2$) where force of infection $\frac{\beta SI_2} {1 + \zeta I_2^2}$ for strain 2. As models with more general incidence functions are considered, the dynamics of the system become more complicated. Models with incidence functions of the form $g(I)h(S)$ have been studied, such as ^[10]. In the most general case, the transmission of the disease may be given by a non-factorable function of $S$ and $I$.

In this paper, our purpose is to study model considered in ^[2] modifying the force of infection in the compartments $I_1$ and $I_2$, by extending the incidence function to a more general form $F(S, I)$, which is based on the incidence rate studied in ^[11]. Our main aim is to mathematically analyze the effect of the vaccine for strain 1, the general incidence rate of strain 1 ($F_1(S, I_1)$), and the general incidence rate of strain 2 ($F_2(S, I_2)$) on the dynamics of the model (2.2).

This paper is organized as follows. In section 2.1, we formulate the model. In section 3.1, we investigate the disease dynamics described by the model. In section 3.2, we calculate the basic reproduction number. In section 3.3, we establish the existence of equilibrium points. In section 3.4, we study the stability of the model. In section 3.5, provides some numeric simulations to illustrate our main theoretical results. The paper ends with some remarks.

2.   Materials and method

2.1.   The model

Rahman and Zou ^[2] proposed a two-strain model with a single vaccination, namely.

where $\lambda = r+\mu$, $\alpha_1 = \gamma_1+v_1+\mu$, $\alpha_2 = \gamma_2+v_2+\mu$. The compartments are $S(t)$, $V_1(t)$, $I_1(t)$, $I_2(t)$ and $R(t)$ which denote the population of susceptible, vaccine of strain 1, infective with respect to strain 1, infective with respect to strain 2 and removed individuals at time t, respectively. We assume that all the parameters are positive constants that can be interpreted as follows:

● $\Lambda$ is the birth rate.

● $\mu$ is the death rate.

● $r$ is the rate of vaccination with strain 1.

● $k$ is the transmission coefficient of vaccinated individuals to strain 2.

● $\beta_1$ is the transmission coefficient of susceptible individuals to strain 1.

● $\beta_2$ is the transmission coefficient of susceptible individuals to strain 2.

● $\frac{1}{\gamma_1}$ is the average infection period of strain 1.

● $\frac{1}{\gamma_2}$ is the average infection period of strain 2.

● $v_1$ is the infection-induced death rate of strain 1.

● $v_2$ is the infection-induced death rate of strain 2.

The modification of the model (2.1) is given by the following system:

Whose state space is $\mathbb{R}^5_+ = \{ (S, V_1, I_1, I_2, R) : S \geq 0, V_1 \geq 0, I_1\geq 0, I_2 \geq 0, R\geq 0 \}$ and subject to the initial conditions $S(0) = S_0 \geq 0$, $V_1(0) = V_{10} \geq 0$, $I_1(0) = I_{10} \geq 0$, $I_2(0) = I_{20} \geq 0$ and $R(0) = R_0 \geq 0$.

We make the following hypotheses on $F_i$, $i = 1, 2.$:

H1) $F_i(S, I_i) = I_if_i(S, I_i)$ with $F_i$, $f_i\in \mathbf{C}^2(\mathbb{R}^2_+ \to \mathbb{R}_+)$ and $F(0, I_i) = F(S, 0) = 0$ for all $S, I_i\geq 0.$

H2) $\frac{\partial f_i}{\partial S}(S, I_i) > 0$ and $\frac{\partial f_i}{\partial I_i}(S, I_i)\leq 0$ for all $S, I_i\geq 0.$

H3) $\lim\limits_{I_i\to 0^+} \frac{F_i(S, I_i)}{I_i}$ exists and is positive for all $S > 0.$

The first of this hypotheses is a basic requirement for any biologically feasible incidence rate, since the disease cannot spread when the number of susceptible or infected individuals is zero.

As for (H2), the condition $\frac{\partial f_i}{\partial S} (S, I_i) > 0$ ensures the monotonicity of $f_i (S, I_i)$ on S, while $\frac{\partial f_i}{\partial I_i} (S, I) \leq 0$ suggests that $\frac{F_i(S, I_i)}{I_i}$ is non-increasing with respect to $I_i$. In the case when $f_i$ monotonically increases with respect to both variables and is concave with respect to $I_i$, the hypothesis (H2) naturally holds. Concave incidence functions have been used to represent the saturation effect in the transmission rate when the number of infected is very high and exposure to the disease is virtually certain.

(H3) is needed only to ensure that the basic reproduction number is well defined. Some examples of incidence functions studied in the literature that satisfy (H1)–(H3) are as follows:

(C1) F(S, I) = $\beta SI$ ^[2].

(C2) F(S, I) = $\frac{\beta SI} {1 + \zeta S}$, where $\zeta\geq 0$ describes the psychological effect of general public towards the infective ^[8].

(C3) F(S, I) = $\frac{\beta SI} {1 + \zeta I^2}$, where $\zeta\geq 0$ measures the psychological or inhibitory effect of the population ^[9].

A more thorough list can be found in ^[11]. It should be noted that model (2.2) extends as well as generalizes many special cases.

3.   Results

3.1.   Disease dynamics described by the model

Lemma 1. Under the initial value $(S_0, V_{10}, I_{10}, I_{20}, R_0)\in \mathbb{R}^5_+$ the system (2.2) has a unique positive and bounded solution in $\mathbb{R}^5_+$ for $t > 0$. All solutions ultimately enter and remain in the following bounded and positively invariant region

Proof. The right hand side of system (2.2) is continuous and satisfies the Lipschitz condition in $\mathbb{R}^5_+$. Then the system (2.2) has a unique solution $(S(t, V_1(t), I_1(t), I_2(t), R(t))$ in $[0, t_m)$ for some $t_m > 0$. Adding all equations in (2.2), the total population $N = S+V_1+I_1+I_2+R$ satisfies:

The comparison theorem implies that $\lim\limits_{t \to \infty} \sup N(t)\leq \frac{\Lambda}{\mu}$. Hence $N(t)$ is bounded and so are all components $S(t)$, $V_1(t)$, $I_1(t)$, $I_2(t)$ and $R(t)$. This in turn shows that the solution exists globally, i.e. for all $t\geq 0$. Consequently, the solutions $S(t)$, $V_1(t)$, $I_1(t)$, $I_2(t)$, $R(t)$ of (2.2) are ultimately bounded in the positively invariant region $\Omega$.

Let $(S(t), \ V_1(t), \ I_1(t), \ I_2(t), \ R(t))$ be a solution of system (2.2) with positive initial conditions. Assume by contradiction that there exists $t > 0$ such that $S(t)\leq 0$, $V_1(t)\leq 0$, $I_1(t)\leq 0$, $I_2(t)\leq 0$ or $R(t)\leq 0$. By continuity of solutions, this implies that there is a minimal $t_0 > 0$ such that $S(t_0)$, $V_1(t_0)$, $I_1(t_0)$, $I_2(t_0)$ or $R(t_0)$ is zero.

If $S(t_0) = 0$, then $\dot{S} = \Lambda > 0$ at $t_0$, so $S$ is increasing in a neighbourhood $(t_0-\epsilon, \ t_0+\epsilon)$ of $t_0$. Thus $S(t_0-\frac{\epsilon}{2}) < S(t_0) = 0$, and since $S(0) > 0$ and $S(t_0-\frac{\epsilon}{2}) < 0$, there exists a $t_1\in\left(0, \ t_0-\frac{\epsilon}{2} \right)$ with $S(t_1) = 0$. But $t_1 < t_0$, which contradicts the minimality of $t_0$, then $S(t_0) > 0$.

If $V_1(t_0) = 0$, then $\dot{V_1} = rS(t_0) > 0$ at $t_0$. So $V_1$ is increasing in a neighbourhood $(t_0-\epsilon, \ t_0+\epsilon)$ of $t_0$. Thus $V_1(t_0-\frac{\epsilon}{2}) < V_1(t_0) = 0$, and since $V_1(0) > 0$ and $V_1(t_0-\frac{\epsilon}{2}) < 0$, there exists a $t_1\in\left(0, \ t_0-\epsilon/2\right)$ with $V_1(t_1) = 0$. But $t_1 < t_0$, which contradicts the minimality of $t_0$, then $V_1(t_0) > 0$.

If $I_1(t_0) = 0$, then $\dot{I_1} = 0$ at $t_0$. On the other hand, any solution with $I_1(0) = 0$ satisfies $I(t) = 0$ for all $t > 0$. Since $I_1(0) > 0$ and $I_1(t_0) = 0$, this contradicts the uniqueness of solutions. Similar contradictions are obtained if we assume that $I_2(t_0) = 0$ or $R(t_0) = 0$. Thus we conclude that the solutions of (2.2) are positive for all $t > 0$.

Since the equation for $\dot{R}$ is actually decoupled from the rest in Eq (2.2), we only need to consider dynamics of the following four-dimensional sub-system:

3.2.   Basic reproduction number

The basic reproduction number is a dimensionless quantity denoted by $\mathcal{R}_0$. It is defined as the expected number of secondary infection cases caused by a single typical infective case during its entire period of infectivity in a wholly susceptible population. Then, referring to the method of ^[12].

Then

where $E_0 = (S^0, V_{1}^0, 0, 0) = \left(\frac{\Lambda}{\lambda}, \frac{r \Lambda}{\mu \lambda}, 0, 0 \right)$. The matrix F is non-negative and is responsible for new infections (transmission matrix), while the V is invertible and is referred to as the transition matrix for the model (3.1). It follows that,

where $\sigma_i = \frac{\partial F_i(S^0, 0)}{\partial I_i}$, for $i = 1, 2$. Thus, the basic reproduction number can be calculated as

where $\rho(A)$ denotes the spectral radius of a matrix A. Let

Then

Therefore $\mathcal{R}_1, \mathcal{R}_2\leq \mathcal{R}_0.$

3.3.   Existence of equilibrium solutions

The four possible equilibrium points for the system (3.1) are: Disease-free equilibrium, single-strain $(I_1)$-infection, single-strain $(I_2)$-infection and endemic equilibrium. The system (3.1) has disease-free equilibrium $E_0 = \left(\frac{\Lambda}{\lambda}, \frac{r \Lambda}{\mu \lambda}, 0, 0 \right)$ for all parameter values. We will now prove the existence of the other equilibrium points. First we will show some lemmas.

Lemma 2. For i = 1, 2.

Also:

Proof. By H1)

Then

By H2) $\frac{\partial f_i(S, I_i)}{\partial I_i}\leq 0$, then:

Lemma 3. For model (3.1), the closed set $\Omega_1 = \left\lbrace (S, V_1, I_1, I_2)\in \Omega | S \leq S^0 \ \text{and} \ V_1\leq V_1^0 \right\rbrace$ is a positively invariant set.

Proof. As $\Omega$ is a positively invariant set for model (3.1), it will be enough to show that if $S = S^0$, then $\dot{S}\leq 0$ and if $S\leq S^0$ and $V_1 = V_1^0$, then $\dot{V_1}\leq 0.$

If $S = S^0$, then

If $S\leq S^0$ and $V_1 = V_1^0$, Then

Lemma 4. For i = 1, 2.

Proof. By (H1)

Then

Remark 1. By (H2) given $a$ and $b$, if $S \leq a$ and $I_i \geq b$, then $f_i(S, I_i)\leq f_i(a, b)$, $i = 1, 2$.

Theorem 1. (1) The model (3.1) admits a unique single-strain $(I_1)$-infection equilibrium $E_1 = (\bar{S}, \bar{V_1}, \bar{I_1}, 0)$ if and only if $\mathcal{R}_1 > 1$.

(2) The model (3.1) admits a unique single-strain $(I_2)$-infection equilibrium $E_2 = (\tilde{S}, \tilde{V_1}, 0, \tilde{I_2})$ if and only if $\mathcal{R}_2 > 1$.

Proof. $(1)$ If $I_2 = 0$ and $\mathcal{R}_1 > 1$, we consider the system

By (3.3) and (3.4)

Substituting in (3.2).

Note that $\bar{S}\geq 0$ if and only if $\bar{I_1}\leq \frac{\Lambda}{\alpha_1}$. $\bar{I_1}$ being determined by the positive roots of the equation.

Then

And

Therefore $G(\bar{I_1}) > 0$ by $I_1$ sufficiently small. Also

Then Eq (3.5) has a positive root. Also if $E_1$ exists then

Note that $\bar{S} < S^0$. Then by Lemma 2 and remark 1

Then $\mathcal{R}_1 > 1.$

Next, we shall show that $\bar{I_1}$ is unique. From (3.4), it follows that

Using (H2) and Lemma 2, we have that $\frac{-\alpha_1}{\lambda} \frac{\partial F_1(\bar{S}, \bar{I_1})}{\partial S} < 0$ and $\bar{I_1}\frac{\partial f_1(\bar{S}, \bar{I_1})}{\partial I_1}\leq 0$. Furthermore, it can be found that

Which implies that $G(\bar{I_1})$ strictly decreases at any of the zero points of (3.5). Let us suppose that (3.5) has more than one positive root. Without loss of generality, we choose the one, denoted by $\bar{I_1}^*$, that is the nearest to $\bar{I_1}$. Because of the continuity of $G(\bar{I_1})$, we must have $G'(\bar{I_1}^*)\geq 0$, which results in a contradiction with the strictly decreasing property of $G(\bar{I_1})$ at all the zero points.

$(2)$ If $I_1 = 0$ and $\mathcal{R}_2 > 1$, we consider the system

By (3.7) and (3.8)

Substituting in (3.6).

Note that $\tilde{S}\geq 0$ if and only if $\tilde{I_2}\leq \frac{\Lambda}{\alpha_2}$. $\tilde{I_2}$ being determined by the positive roots of the equation.

Then

And

Therefore $H(\tilde{I_2}) > 0$ by $\tilde{I_2}$ sufficiently small. Also

Then Eq (3.9) has a positive root. Also if $E_2$ exists then

Note that by H1 we have $-F_2(\tilde{S}, \tilde{I_2}) < 0$, then $\Lambda-\lambda \tilde{S} > 0$, therefore $\tilde{S} < S^0$ and $\tilde{V_1} < V_1^0$. Then by Lemma 2 and remark 1

Then $\mathcal{R}_2 > 1.$

Next, we shall show that $\tilde{I_2}$ is unique. From (3.8), it follows that

Furthermore, it can be found that

If $-\alpha_2 r \mu-\alpha_2 \mu^2+k\Lambda r\leq 0$, then $H'(\tilde{I_2}) < 0$ which implies that $H(\tilde{I_2})$ strictly decreases at any of the zero points of (3.9). Let us suppose that (3.9) has more than one positive root. Without loss of generality, we choose the one, denoted by $\tilde{I_2}^*$, that is the nearest to $\tilde{I_2}$. Because of the continuity of $H(\tilde{I_2})$, we must have $H'(\tilde{I_2}^*)\geq 0$, which results in a contradiction with the strictly decreasing property of $H(\tilde{I_2})$ at all the zero points.

If $-\alpha_2 r \mu-\alpha_2 \mu^2+k\Lambda r > 0$. Next, we show that $\tilde{I_2}\notin \left[0, \frac{-r\alpha_2-\alpha_2\mu +\sqrt{r\alpha_2(r\alpha_2+\alpha_2 \mu+k\Lambda)}}{\alpha_2k}\;\;\right)$. Note that

Then

If $\tilde{I_2} \in \left[0, \frac{-r\alpha_2-\alpha_2\mu +\sqrt{r\alpha_2(r\alpha_2+\alpha_2 \mu+k\Lambda)}}{\alpha_2k}\right)$, then $\tilde{S}'(\tilde{I_2}) > 0$, therefore $\tilde{S}\geq S^0$, which results in a contradiction, since $\tilde{S} < S^0$.

Thus $\tilde{I_2} \in \left[\frac{-r\alpha_2-\alpha_2\mu +\sqrt{r\alpha_2(r\alpha_2+\alpha_2 \mu+k\Lambda)}}{\alpha_2k}, \frac{\Lambda}{\alpha_2}\right]$, which implies that $H(\tilde{I_2})$ strictly decreases at any of the zero points of (3.9). Let us suppose that (3.9) has more than one positive root in $\left[\frac{-r\alpha_2-\alpha_2\mu +\sqrt{r\alpha_2(r\alpha_2+\alpha_2 \mu+k\Lambda)}}{\alpha_2k}, \frac{\Lambda}{\alpha_2}\right]$. Without loss of generality, we choose the one, denoted by $\tilde{I_2}^*$, that is the nearest to $\tilde{I_2}$. Note that $H'(\tilde{I_2}^*) < 0$ and $H'(\tilde{I_2}) < 0$. Because of the continuity of $H(\tilde{I_2})$, we must have $H'(\tilde{I_2}^*)\geq 0$, which results in a contradiction.

The model (3.1) can have endemic infection equilibrium $E_3 = (S^*, V_1^*, I_1^*, I_2^*)$. To find $E_3$, we consider the system

By (3.11), (3.12) and (3.13)

Substituting in (3.10).

Note that $S^* \geq 0$ if and only if $I_1^* \leq \frac{\Lambda-\alpha_2 I_2^*}{\alpha_1}$ and $I_2^* \leq \frac{\Lambda-\alpha_1 I_1^*}{\alpha_2}$. $\bar{I_2}$ being determined by the positive roots of the equation.

$I_1^*$ being determined by the positive roots of the equation.

3.4.   Stability of equilibrium

In this section we will study the local and global stability of the equilibrium points.

Theorem 2. The disease-free equilibrium $E_0 = \left(\frac{\Lambda}{\lambda}, \frac{r \Lambda}{\mu \lambda}, 0, 0 \right)$ is unstable if $\mathcal{R}_0 > 1$ while it is locally asymptotically stable if $\mathcal{R}_0 < 1$.

Proof. The Jacobian matrix of the model we get is the following one

Then Eq (3.14) at the disease-free equilibrium $E_0$ is

Thus the eigenvalues of the above Eq (3.15) are

From (3.16), if $\mathcal{R}_0 < 1$, then $\lambda_3, \lambda_4 < 0$ and we obtain that the disease-free equilibrium $E_0$ of Model (3.1) is locally asymptotically stable. If $\mathcal{R}_0 > 1$, then the disease-free equilibrium loses its stability.

Theorem 3. Let $\bar{\mathcal{R}_2} = \frac{1}{\alpha_2}\frac{\partial F_2(\bar{S}, 0)}{\partial I_2}+\frac{k \bar{V_1}}{\alpha_2}$. The equilibrium $E_1$ is unstable if $\bar{\mathcal{R}_2} > 1$ while it is locally asymptotically stable if $\bar{\mathcal{R}_2} < 1$.

Proof. Then Eq (3.14) at the equilibrium $E_1$ is

where

The last equality regarding $A_{33}$, is because Eq (3.4) implies that $f_1(\bar{S}, \bar{I_1})-\alpha_1 = 0$. The corresponding characteristic polynomial is

Then an eigenvalue is $A_{44}$ and the remaining ones satisfy

where

Note that

Then $a_1$, $a_0 > 0$ and

Applying the Routh–Hurwitz criterion, we see that all roots of $x^3+a_2 x^2+a_1 x+ a_0$ have negative real parts. If $\bar{\mathcal{R}_2} > 1$, then $A_{44} > 0$ therefore $E_1$ is unstable and if $\bar{\mathcal{R}_2} < 1$, then $A_{44} < 0$ therefore $E_1$ is stable.

Remark 2. $\bar{S}\leq S^0$ and $\bar{V_1}\leq V_1^0$, then $\bar{\mathcal{R}_2}\leq \mathcal{R}_2$, therefore if $\mathcal{R}_2 < 1$ then $\bar{\mathcal{R}_2} < 1$.

Theorem 4. Let $\tilde{\mathcal{R}_1} = \frac{1}{\alpha_1}\frac{\partial F_1(\tilde{S}, 0)}{\partial I_1}$. If $\frac{\partial F_2(\tilde{S}, \tilde{I_2})}{\partial I_2}\leq 0$ the equilibrium $E_2$ is unstable if $\tilde{\mathcal{R}_1} > 1$ while it is locally asymptotically stable if $\tilde{\mathcal{R}_1} < 1$.

Proof. Then Eq (3.14) at the equilibrium $E_1$ is

Where

The last equality regarding $B_{44}$, is because Eq (3.8) implies that $k \tilde{V_1}-\alpha_2 = -f_2(\tilde{S}, \tilde{I_2})$. The corresponding characteristic polynomial is

Then (3.18) has an eigenvalue equal to $B_{33}$ and the remaining ones satisfy

where

Note that

And

Then $b_1$, $b_0 > 0$. Also

Applying the Routh–Hurwitz criterion, we see that all roots of $x^3+b_2 x^2+b_1 x+ b_0$ have negative real parts. If $\tilde{\mathcal{R}_1} > 1$, then $B_{33} > 0$ therefore $E_2$ is unstable and if $\tilde{\mathcal{R}_1} < 1$, then $B_{33} < 0$ therefore $E_2$ is stable.

Remark 3. $\tilde{S}\leq S^0$, then $\tilde{\mathcal{R}_1}\leq \mathcal{R}_1$, therefore if $\mathcal{R}_1 < 1$ then $\tilde{\mathcal{R}_1} < 1$.

Remark 4. If $\frac{\partial F_2(\tilde{S}, \tilde{I_2})}{\partial I_2} > 0$, then $b_i > 0$ $i = 0, 1, 2$.

Remark 5. The Theorem 4 is valid for $\frac{\partial F_2(\tilde{S}, \tilde{I_2})}{\partial I_2} > 0$ if $b_2b_1-b_0 > 0$.

Theorem 5. If $\bar{\mathcal{R}_2} > 1$ and $\tilde{\mathcal{R}_1} > 1$ then system (3.1) is uniformly persistent.

Proof. The result follows from an application of Theorem 4.6 in ^[13], with $X_1 = \text{int} (\mathbb{R}^4_+)$ and $X_2 = \text{bd}(\mathbb{R}^4 _+)$ this choice is in accordance with the conditions stated in the theorem. Now, note that by of Lemma 1 there exists a compact set $\Omega$ in which all solution of system (3.1) initiated in $\mathbb{R}^4_+$ ultimately enter and remain forever after. The condition $\left(C_ {4.2} \right)$ is easily verified for this set $\Omega_1$. On other hand, we denote the omega limit set of the solution $x(t, x_0)$ of system (3.1) starting in $x_0 \in \mathbb{R}^4_+$ by $w(x_0)$. Note that $w(x_0)$ is bounded (Lemma 1), we need to determine the following set:

From the system equations (3.1) it follows that all solutions starting in $\text{bd}(\mathbb{R}^4_+)$ but not on the $I_1$ axis or on the $I_2$ axis leave $\text{bd}(\mathbb{R}^4 _+)$. This implies that

Furthermore, we see that $\Omega_2 = \{E_0, E_1, E_2\}$, then $\bigcup_{i = 1}^3 \{E_i\}$ is a covering of $\Omega_2$, which is isolated (since $E_i$ $(i = 1, 2, 3)$ is a saddle point) and acyclic. Finally we need to prove that $E_i$ (i = 1, 2, 3) is a weak repeller for $X_ 1$ to end the prove.

By definition $E_i$ is a weak repeller for $X_1$ if for every solution $(S(t), V_1(t), I_1(t), I_2(t))$ starting in $(S_0, V_{10}, I_{10}, I_{20}) \in X_1$

We will first show that $E_0$ is a weak repeller for $X_1$. Since $\bar{\mathcal{R}_2} > 1$ and $\tilde{\mathcal{R}_1} > 1$, then $\mathcal{R}_2 = \frac{1}{\alpha_2}\left(f_2(S^0, 0)+kV^0\right) > 1$ and $\mathcal{R}_1 = \frac{1}{\alpha_1}\left(f_1(S^0, 0)\right) > 1$, therefore $f_2(S^0, 0)+kV^0-\alpha_2 > 0$ and $f_1(S^0, 0)-\alpha_1 > 0$. Because of the continuity of $f_2(S, I_2)+kV_1-\alpha_2$ and $f_1(S, I_1)-\alpha_1$, there exists a sufficiently small constant $\eta_2 > 0$, such that $f_1(S^0-\eta_2, \eta_2)-\alpha_1 > 0$ and $f_2(S^0-\eta_2, \eta_2)+k (V_1^0-\eta_2)-\alpha_1 > 0$.

Now, we suppose that $E_0$ is not a weak repeller for $X_1$, i.e., there exists a solution $(S(t), V_1(t), I_1(t), I_2(t))$ starting in $(S_0, V_{10}, I_{10}, I_{20}) \in X_1$ such that

Then exists $T_1 > 0$ such that for every $\eta_1 > 0$

Let $\eta_1 = \eta_2$, then for $t\geq T_1$.

and

By comparison principle, we have

Note that $f_1(S^0-\eta_2, \eta_2)-\alpha_1 > 0$, $f_2(S^0-\eta_2, \eta_2)+k (V_1^0-\eta_2)-\alpha_1 > 0$, $I_1(T_1) > 0$ and $I_2(T_1) > 0$, which implies that $\lim\limits_{t \to \infty} I_1 = \lim\limits_{t \to \infty} I_2 = \infty$, this gives a contradiction. Then $E_0$ is a weak repeller for $X_1$.

Similarly it is shown that $E_1$ and $E_2$ are weak repeller for $X_1$. Then we conclude that system (3.1) is uniformly persistent.

Further, it is proved in ^[14] uniform persistence implies the existence of an interior equilibrium point. Therefore, we have established the following.

Theorem 6. The model (3.1) admits a endemic equilibrium $E_3 = (S^*, V_1^*, I_1^*, I_2^*)$ if $\bar{\mathcal{R}_2} > 1$ and $\tilde{\mathcal{R}_1} > 1$.

Theorem 7. If $c_1c_2-c_3 > 0$ and $c_1c_2c_3-c_3^2-c_1^2c_4 > 0$, where

Then $E_3$ is locally asymptotically stable.

Proof. Then Eq (3.14) at the equilibrium $E_3$ is

Where

The corresponding characteristic polynomial is

Note that $c_1 > 0$,

then $c_2 > 0$, If $C_{14}\geq 0$ then $c_3 > 0$ and $c_4 > 0$, while if $C_{14} < 0$ we have that

and let

then

Then $c_3 > 0$ and $c_4 > 0$. If $c_1c_2-c_3 > 0$ and $c_1c_2c_3-c_3^2-c_1^2c_4 > 0$ by Routh–Hurwitz criterion, we see that all roots of $x^4+c_1x^3+c_2x^2+c_3x+c_4$ have negative real parts, then $E_3$ is locally asymptotically stable.

3.4.1.   Global stability of equilibria

In this section, we study the global properties of the equilibria. We use Lyapunov function to show the global stabilities. Such Lyapunov functions all take advantage of the properties of the function.

which is positive in $\mathbb{R}_+$ except at $x = 1$, where it vanishes.

Theorem 8. The DFE $E_0$ is globally asymptotically stable if,

Proof. Consider the Lyapunov function

Since $I_1, I_2 > 0$, then $V (S, V_1, I_1, I_2) \geq 0$ and $V (S, V_1, I_1, I_2)$ attains zero at $I_1 = I_2 = 0$.

Now, we need to show $\dot{V}\leq 0$.

For $S\leq S^0$ and $V_1\leq V_1^0$

Furthermore, $\frac{\text{d}V}{\text{d}t} = 0$ if and only if $I_1 = I_2 = 0$, so the largest invariant set contained in $\left\lbrace(S, V_1, I_1, I_2)\in \Omega_1 | \frac{dV}{dt} = 0 \right\rbrace$ is the hyperplane $I_1 = I_2 = 0$, By LaSalle's invariant principle, this implies that all solution in $\Omega_1$ approach the hyperplane $I_1 = I_2 = 0$ as $t\to \infty$. Also, All solution of (3.1) contained in such plane satisfy $\dot{S} = \Lambda-\lambda S$, $\dot{V_1} = rS-\mu V_1$, which implies that $S\to \frac{\Lambda}{\lambda}$ and $V_1 \to \frac{r \Lambda}{\mu \lambda}$ as $t \to \infty$, that is, all of these solution approach $E_0$. Therefore we conclude that $E_0$ is globally asymptotically stable in $\Omega_1$.

Now we will show that every solution $(S(t), V_1(t), I_1(t), I_2(t))\in \mathbb{R}^4_+$, where $t\to \infty$ $(S(t), V_1(t), I_1(t), I_2(t))\in \Omega_1$, let $(S(t), V_1(t), I_1(t), I_2(t))\in \mathbb{R}^4_+$. Then

By the comparison principle $\lim\limits_{t \to \infty} \sup S(t)\leq \frac{\Lambda}{\lambda} = S^0$. Then $S(t)\leq S^0$ for $t$ sufficiently large.

Also if $S(t)\leq S^0$.

By the comparison principle $\lim\limits_{t \to \infty} \sup V(t)\leq \frac{r S^0}{\mu} = V_1^0$. Therefore $E_0$ is globally asymptotically stable.

From now on, we assume that

$H4)$ For $i = 1, 2.$ $f_i(S, I_i) = Sg_i(S, I_i)$.

Lemma 5. Let $a > 0$ be a constant, for $i = 1, 2$ if $\frac{\partial F_i(S, I_i)}{\partial I_i}\geq 0$ for all $I_i$, then

Proof. Note that

If $a\geq I_i$, then

If $a\leq I_i$, then

Therefore

Theorem 9. Suppose that $\frac{\partial F_1(S, I_1)}{\partial I_1}\geq 0$ for all $I_1$, then $E_1$ is globally asymptotically stable if,

Proof. Consider the Lyapunov function

Since $I_2 > 0$, then $V (S, V_1, I_1, I_2) \geq 0$ and $V (S, V_1, I_1, I_2)$ attains zero at $I_2 = 0$. Now, we need to show $\dot{V}\leq 0$.

For $S\leq S^0$ and $V_1\leq V_1^0$

Furthermore, $\frac{\text{d}V}{\text{d}t} = 0$ if and only if $I_2 = 0$. Suppose that $(S(t), V_1(t), I_1(t), I_2(t))$ is a solution of (3.1) contained entirely in the set $M = \{(S(t), V_1(t), I_1(t), I_2(t)) \in \Omega_1|\dot{V} = 0 \}$. Then, $\dot{I_2} = 0$ and, from the above inequalities, we have $I_2 = 0$. Thus, the largest positively invariant set contained in $M$ is the plane $I_2 = 0$. By LaSalle’s invariance principle, this implies that all solutions in approach the plane $I_2 = 0$ as $t\to \infty$. On the other hand, solutions of (3.1) contained in such plane satisfy

Now we will show that $S(t)\to \bar{S}$, $V_1(t)\to \bar{V_1}$ and $I_1(t)\to \bar{I_1}$ Consider the Lyapunov function

Note that $1-\frac{F_1(\bar{S}, \bar{I_1})}{F_1(\chi, \bar{I_1})} = \frac{\bar{I_1}(f_1(\chi, \bar{I_1})-f_1(\bar{S}, \bar{I_1}))}{F_1(\chi, \bar{I_1})}$, by H2) $f_1(S, \bar{I_1})-f_1(\bar{S}, \bar{I_1})\geq 0$ if $S\geq \bar{S}$ and $f_1(S, \bar{I_1})-f_1(\bar{S}, \bar{I_1})\leq 0$ if $S\leq \bar{S}$, then $\int^S_{\bar{S}} \left(1-\frac{F_1(\bar{S}, \bar{I_1})}{F_1(\chi, \bar{I_1})} \right)\text{d}\chi\geq 0$ for all $S$. Therefore, $V (S, V_1, I_1) \geq 0$ and $V (S, V_1, I_1)$ attains zero at $S(t) = \bar{S},$ and $I_1(t) = \bar{I_1}$.

Now, we need to show $\dot{V}\leq 0$.

Note that

and

The last inequality is due to the Lemma 5 and the relation of the geometric and arithmetic means, then $\dot{V}\leq 0$. Furthermore, $\frac{\text{d}V}{\text{d}t} = 0$ if and only if $S = \bar{S}$ and $I_1 = \bar{I_1}$, which implies that $S\to \bar{S}$, $I_1 \to \bar{I_1}$ and $I_2\to 0$ as $t \to \infty$. By LaSalle's invariant principle, this implies that all solutions in $\Omega_1$ approach the plane $S = \bar{S}$, $I_1 = \bar{I_1}$ and $I_2 = 0$ as $t\to \infty$. Also, All solutions of (3.1) contained in such plane satisfy $\dot{V_1} = r\bar{S}-\mu V_1$, which implies that $V_1 \to \frac{r\bar{S}}{\mu} = \bar{V_1}$ as $t \to \infty$, that is, all of these solution approach $E_1$. Therefore we conclude that $E_1$ is globally asymptotically stable in $\Omega_1$.

Corollary 1. If $2-\frac{F_1(\bar{S}, \bar{I_1})}{F_1(S, \bar{I_1})}+\frac{F_1(S, I_1)}{F_1(S, \bar{I_1})} -\frac{I_1}{\bar{I_1}} - \frac{\bar{I_1}f_1(S, I_1)}{F_1(\bar{S}, \bar{I_1})}\leq 0$ and $\mathcal{R}_2 < 1$ then $E_1$ is globally asymptotically stable.

Theorem 10. Suppose that $\frac{\partial F_2(S, I_2)}{\partial I_2}\geq 0$ for all $I_2$, then $E_2$ is globally asymptotically stable if,

Proof. Consider the Lyapunov function

Since $I_1 > 0$, then $V (S, V_1, I_1, I_2) \geq 0$ and $V (S, V_1, I_1, I_2)$ attains zero at $I_1 = 0$. Now, we need to show $\dot{V}\leq 0$.

For $S\leq S^0$

Furthermore, $\frac{\text{d}V}{\text{d}t} = 0$ if and only if $I_1 = 0$. Suppose that $(S(t), V_1(t), I_1(t), I_2(t))$ is a solution of (3.1) contained entirely in the set $M = \{(S(t), V_1(t), I_1(t), I_2(t)) \in \Omega_1|\dot{V} = 0 \}$. Then, $\dot{I_1} = 0$ and, from the above inequalities, we have $I_1 = 0$. Thus, the largest positively invariant set contained in $M$ is the plane $I_1 = 0$. By LaSalle’s invariance principle, this implies that all solutions in approach the plane $I_1 = 0$ as $t\to \infty$. On the other hand, solutions of (3.1) contained in such plane satisfy.

Now we will show that $S(t)\to \tilde{S}$, $V_1(t)\to \tilde{V_1}$ and $I_1(t)\to \tilde{I_1}$ Consider the Lyapunov function

Now, we need to show $\dot{V}\leq 0$.

Note that

The last inequality is due to the Lemma 5 and the relation of the geometric and arithmetic means, then $\dot{V}\leq0$. Furthermore, $\dot{V} = 0$ if and only if $S = \tilde{S}$, $I_2 = \tilde{I_2}$ and $V_1 = \tilde{V_1}$. Therefore $E_2$ is globally asymptotically stable.

Remark 6. Note that if $\frac{\partial g_2(S, I_2)}{\partial S}\geq 0$, then

Corollary 2. If $2-\frac{F_2(\tilde{S}, \tilde{I_2})}{F_2(S, \tilde{I_2})}+\frac{F_2(S, I_2)}{F_2(S, \tilde{I_2})} -\frac{I_2}{\tilde{I_2}} - \frac{\tilde{I_2}f_2(S, I_2)}{F_2(\tilde{S}, \tilde{I_2})}\leq 0$, $\mathcal{R}_1 < 1$ and $2-\frac{F_2(\tilde{S}, \tilde{I_2})}{F_2(S, \tilde{I_2})}+\frac{S F_2(\tilde{S}, \tilde{I_2})}{\tilde{S} F_2(S, \tilde{I_2})}-\frac{V_1}{\tilde{V_1}}-\frac{S\tilde{V_1}}{\tilde{S}V_1}\leq 0$ then $E_2$ is globally asymptotically stable.

Theorem 11. $E_3$ is globally asymptotically stable if

Proof. Assume $E_3$ exists. Consider the Lyapunov function

Where $g(x) = x -1 -ln(x)$. Then $V (S, V_1, I_1, I_2) \geq 0$ and $V (S, V_1, I_1, I_2)$ attains zero at $E_3$.

Now, we need to show $\dot{V}\leq 0$.

By the relation of geometric and arithmetic means, we conclude $\dot{V}\leq 0$, with equality holding only at the equilibrium $E_3$. Therefore $E_3$ is globally asymptotically stable.

3.5.   Numerical simulations

In this section, we present some numerical simulations of the solutions for system (3.1) to verify the results obtained in section 3.3 and give examples to illustrate theorems in section 3.4. In system (3.1), we set:

In this case

Parameters and units are arbitrary and have been used for illustration purposes only. Anyway, when considering a realistic scenario such values could be derived from statistical data.

$•$ Example 6.1. In system (3.1), we set $\beta_1 = 0.00003$, $r = 0.1$, $\beta_2 = 0.0002$, $\zeta_1 = 0.7$ and $\zeta_2 = 0.9$. Then $S^0\approx 1667$, $V^0\approx 8333$, $\mathcal{R}_1\approx0.2632$ $\mathcal{R}_2\approx0.7947$. By Theorem 8, we see that the disease-free equilibrium $E_0$ is globally asymptotically stable. Numerical simulation illustrates our result (see Figure 1).

$•$ Example 6.2. In system (3.1), we set $\beta_1 = 0.0002$, $r = 0.1$, $\beta_2 = 0.0002$, $\zeta_1 = 0$ and $\zeta_2 = 0.9$. Then $\bar{S}\approx 950$, $\bar{V_1}\approx 4750$, $\bar{I_1}\approx 453$, $\mathcal{R}_1\approx1.7544$, $\mathcal{R}_2\approx0.7947$. By Theorem 9, we see that the $E_1$ is globally asymptotically stable. Numerical simulation illustrates our result (see Figure 2).

$•$ Example 6.3. In system (3.1), we set $\beta_1 = 0.00003$, $r = 0.1$, $\beta_2 = 0.0002$, $\zeta_1 = 0.7$ and $\zeta_2 = 0.001$. Then $\tilde{S}\approx 1317$, $\tilde{V_1}\approx 4814$, $\tilde{I_2}\approx 368$, $\mathcal{R}_1\approx 0.2632$, $\mathcal{R}_2\approx 1.3889$ and $2-\frac{F_2(\tilde{S}, \tilde{I_2})}{F_2(S, \tilde{I_2})}+\frac{S F_2(\tilde{S}, \tilde{I_2})}{\tilde{S} F_2(S, \tilde{I_2})}-\frac{V_1}{\tilde{V_1}}-\frac{S\tilde{V_1}}{\tilde{S}V_1}\leq 0$ (see Figure 3). By corollary 2, we see that the $E_2$ is globally asymptotically stable. Numerical simulation illustrates our result (see Figure 4).

$•$ Example 6.4. In system (3.1), we set $\beta_1 = 0.0002$, $r = 0.01$, $\beta_2 = 0.0002$, $\zeta_1 = 0.0001$ and $\zeta_2 = 0.0001$. Then $\mathcal{R}_1\approx 7.0175$, $\mathcal{R}_2\approx 4.1270$, $\tilde{S}\approx 1134$, $\bar{S}\approx 5310$, $\bar{V_1} \approx 2655$, $\bar{R_2}\approx 3.555$ and $\tilde{R_1}\approx 1.194$. Then by Theorem 6, $E_3 = (S^*, V_1^*, I_1^*, I_2^*)$ exists ($S^*\approx 1133$, $V_1^*\approx 320$, $I_1^* \approx 44$, $I_2^* \approx 774$), Also $c_1\approx 0.2501$ $c_2\approx 0.0171$ $c_3\approx 3.4759\times 10^{-04}$ $c_4\approx 3.4759\times 3.9242 10^{-06}$, $c_1 c_2-c_3^2\approx 0.0043$ and $c_1c_2c_3-c_3^2-c_1^2c_4\approx 1.1218e\times 10^{-06}$ by Theorem 7, $E_3$ is locally asymptotically stable. Also $E_3$ satisfies $F_1(S^*, I_1^*)\left(2 - \frac{S^*}{S}- \frac{S g_1(S, I_1)}{S^* g_1(S^*, I_1^*)} \right)+F_2(S^*, I_2^*) \left(2-\frac{S^*}{S}-\frac{S g_2(S, I_2)}{S^* g_2(S^*, I_2^*)}\right)+\mu S^* \left(2- \frac{S^*}{S}-\frac{S}{S^*} \right) +I_1 \left(S^* g_1(S, I_1)-\alpha_1 \right)+I_2 \left({S}^* g_2(S, I_2)+ k {V_1}^*-\alpha_2 \right) < 0$. By Theorem 11, we see that the $E_3$ is globally asymptotically stable. Numerical simulation illustrates our result (see Figure 5).

4.   Conclusions and discussions

In this paper, we studied a system of ordinary differential equations to model the disease dynamics of two strains of influenza with only one vaccination for strain 1 being implemented, and general incidence rate for strain 1 and strain 2. We obtained four equilibrium points:

● $E_0$ disease-free equilibrium, $I_1$ and $I_2$ are both zero.

● $E_1$ single-strain-infection equilibria, $I_2$ are zero.

● $E_2$ single-strain-infection-equilibria, $I_1$ are zero.

● $E_3$ double-strain-infection equilibrium, $I_1$ and $I_2$ are both positive.

We have investigated the topics of existence and non-existence of equilibrium points and their stabilities. We also used the next-generation matrix method to obtain two threshold quantities $\mathcal{R}_1$ and $\mathcal{R}_2$, called the basic reproduction ratios for strain 1 and 2 respectively. It was shown that the global stability of each of the equilibrium points depends on the magnitude of these threshold quantities. More precisely, we have proved the following:

● If $\mathcal{R}_0 < 1$ the disease free equilibrium $E_0$ is globally asymptotically stable and if $\mathcal{R}_0 > 1$, then $E_0$ is unstable.

● If $\mathcal{R}_1 > 1$ the model (3.1) admits a unique single-strain-infection-equilibria $E_1$. Also if $\mathcal{R}_2 < 1$ then $E_1$ is globally asymptotically stable and if $\bar{\mathcal{R}_2} > 1$, then $E_1$ is unstable.

● If $\mathcal{R}_2 > 1$ the model (3.1) admits a unique single-strain-infection equilibria $E_2$. Also if $\mathcal{R}_1 < 1$ and $2-\frac{F_2(\tilde{S}, \tilde{I_2})}{F_2(S, \tilde{I_2})}+\frac{S F_2(\tilde{S}, \tilde{I_2})}{\tilde{S} F_2(S, \tilde{I_2})}-\frac{V_1}{\tilde{V_1}}-\frac{S\tilde{V_1}}{\tilde{S}V_1} < 0$, then $E_2$ is globally asymptotically stable and if $\tilde{\mathcal{R}_1} > 1$, then $E_2$ is unstable.

● If $\bar{\mathcal{R}_2} > 1$ and $\tilde{\mathcal{R}_1} > 1$ the model (3.1) admits a double strain infection equilibrium $E_3$. Also if $F_1(S^*, I_1^*)\left(2 - \frac{S^*}{S}- \frac{S g_1(S, I_1)}{S^* g_1(S^*, I_1^*)} \right)+F_2(S^*, I_2^*) \left(2-\frac{S^*}{S}-\frac{S g_2(S, I_2)}{S^* g_2(S^*, I_2^*)}\right)$ $+r S^*\left(3-\frac{S^*}{S}-\frac{V_1}{V_1^*}-\frac{S V_1^*}{S^* V_1}\right)+\mu S^* \left(2- \frac{S^*}{S}-\frac{S}{S^*} \right) +I_1 \left(S^* g_1(S, I_1)-\alpha_1 \right)+I_2 \left({S}^* \right.$ $\left.g_2(S, I_2)+k {V_1}^*-\alpha_2 \right) < 0$. Then $E_3$ is globally asymptotically stable.

In order to discuss the meaning of our mathematical results, let us rewrite the two key indirect parameters $\mathcal{R}_1$ and $\mathcal{R}_2$ in terms of the rate of vaccination ($r$), the incidence rate of strain 1 ($F_1(S, I_1)$) and the incidence rate of strain 2 ($F_2(S, I_2)$) as shown below:

Also, the derivative of $\mathcal{R}_2$ with respect to $r$ is,

Note that $\mathcal{R}_1(r)$ is decreasing and $\mathcal{R}_2(r)$ depends on $\frac{\partial f_2\left(\frac{\Lambda}{\mu}, 0\right)}{\partial S}$. Now we will analyse some cases of incidence rate.

(C1) $F_i(S, I) = \beta_i SI_i$, then $\frac{\partial f_2\left(\frac{\Lambda}{\mu}, 0\right)}{\partial S} = \beta_i$.

(C2) $F_i(S, I) = \frac{\beta_i SI_i} {1 + \zeta_i S}$, then $\frac{\partial f_2\left(\frac{\Lambda}{\mu}, 0\right)}{\partial S} = \frac{\beta_i} {1 + \zeta_i \left(\frac{\Lambda}{r+\mu} \right)}$.

(C3) $F_i(S, I) = \frac{\beta_i SI_i} {1 + \zeta_i I_i^2}$, then $\frac{\partial f_2\left(\frac{\Lambda}{\mu}, 0\right)}{\partial S} = \beta_i$.

Note that for (C1) and (C3), $\mathcal{R}_2(r)$ is increasing if $\beta_i < k$, $\mathcal{R}_2(r)$ is decreasing if $\beta_i > k$ and $\mathcal{R}_2(r)$ is constant if $\beta_i = k$. For (C2), $\mathcal{R}_2(r)$ is increasing if $\beta_i\leq k$ ($\zeta\neq 0$). If $\beta_i > k$ $\mathcal{R}_2(r)$ is increasing if $\frac{\zeta_i k \Lambda}{\beta_i-k}-\mu < r$ and decreasing if $\frac{\zeta_i k \Lambda}{\beta_i-k}-\mu > r$.

Furthermore, if the force of infection of strain 1 is (C2), then $\mathcal{R}_1 = \frac{\beta_1}{\alpha_1(1 + \zeta_1 S^0)}$, note that $\mathcal{R}_1$ is decreasing in $\zeta_1$. If the force of infection of strain 2 is (C2), then $\mathcal{R}_2 = \frac{\beta_2 \Lambda}{\alpha_2(\lambda+\zeta_2 \Lambda)}+\frac{k r \Lambda}{\alpha_2 \mu \lambda}$, note that $\mathcal{R}_2$ is decreasing in $\zeta_2$.

With the above information and the results in section 3.4, we conclude that the vaccination is always beneficial for controlling strain 1, its impact on strain 2 depends on the force of infection of strain 2. For example, if the force of infection of strain 2 is (C2), the impact of vaccination depends on values of $\beta_2$, $k$ and $\zeta_2$. If $\zeta_2 = 0$ and $\beta_2 > k$ it plays a positive role and if $\zeta_2 = 0$ and $\beta_2 < k$, it has a negative impact in controlling strain 2. This is reasonable because larger $k$ (than $\beta_2$) means that vaccinated individuals are more likely to be infected by strain 2 than those who are not vaccinated, and thus, is helpful to strain 2. Smaller $k$ (than $\beta_2$) implies the opposite. If $\zeta_2 \neq 0$ and $\beta_2 \leq k$, it plays a negative role and if $\zeta_2 \neq 0$ and $\beta_2 > k$, not necessarily has a positivity impact in controlling strain 2. This is reasonable because larger $k$ (than $\beta_2$) means that vaccinated individuals are more likely to be infected by strain 2 than those who are not vaccinated, but if $\zeta_2$ is large it means that the population is taking precautions to avoid the infection of strain 2. Also, we conclude that $\zeta_1$ (of the force of infection (C2)) is always beneficial for controlling strain 1 and $\zeta_2$ (of the force of infection (C2)) is always beneficial for controlling strain 2, it means that it is very important that people are taking precautions not to get infected.

Acknowledgments

This work was supported by Sistema Nacional de Investigadores (15284) and Conacyt-Becas.

Conflict of interest

No conflict of interest.