A Model for the Transmission of Covid-19 in a Mass Gathering

Mario Cavani

Instituto Superior de Formación Docente Salomé Ureña, Recinto “Félix Evaristo Mejía”, Santo Domingo, República Dominicana.

DOI: 10.4236/jamp.2020.811194                

Abstract

Recently the pandemic disease Covid-19 has spread rapidly around the world, the disease has some similar characteristics with a previous endemic disease, the Middle East Respiratory Syndrome Coronavirus (MERS-CoV), which has received recently the attention of many researchers. Here, a two patch SEIR model is introduced in order to explain the importance of the impact of the mass gathering of susceptible individuals when the infected individuals and the exposed individuals are taking separately its role in the dynamics of those diseases that involve the two groups. One or more endemic equilibrium points could appear when the reproduction number exceeds the unity, and the model undergoes a forward bifurcation. A re-formulation of the model as an optimal control problem will permit to evaluate the impact of adequate control strategies for the disease.

Keywords

Covid-19, MERS-CoV, Epidemics, Differential Equations, Optimal Control

1. Introduction

This article proposes a model for virus transmission such as those that have been reported in the current Covid-19 pandemic or in the MERS-CoV epidemic. Two populations are considered and a SEIR model is proposed in each population. It is assumed that at the time of the meeting of the two groups there are no people in quarantine nor is it well known that the infection is in process. The disease evolves in such a way that the viruses are efficiently transmitted by mean of close person-to-person contact, it is presumed that this occurs through droplets of respiratory secretions, or by contact with contaminated surfaces. A previous model given in [1] analyzes the lives of millions of pilgrims going to Mecca for Hajj and Umrah (due to the associated mass meetings) are threatened by MERS-CoV.

Recently, in December 2019, a new virus, SARS-CoV-2 appears for the first time, and an outbreak epidemic disease has evolved. The virus has spread rapidly around the world and a very complicated pandemic disease has come out named Covid-19. A similitude in relation to the transmission dynamics of MERS-CoV during Hajj and Umrah is established. The situation has to do with a similar process of mass gathering in Italy. Specifically, in this opportunity of the Covid-19 pandemic on his early stages, according to [2], a mass gathering that spread the disease in the city of Bergamo. This city is located in a province of the same name, with a total population of just over 1.1 million. On March 10, 2020, the day of Atalanta’s return game against Valencia, the number of confirmed coronavirus cases in the province increased to 1472. In all of Lombardy, the region containing Bergamo and also Milan, there have already been 468 deaths. The feeling of an accelerated crisis was present before Atalanta flew to Spain. Approximately 40,000 fans made the 35-mile journey from Bergamo to Milan for the Valencia game. The children were pulled out of school, with mischievous notes from parents to teachers explaining that their offspring needed to participate in a “cultural-historical” moment in international headlines. Just playing in such a game, for an Atalanta media club, was cause for celebration. The “Corriere della Sera” newspaper has published on February 19, 2020 that 40,000 Bergamaschi (people from Bergamo) went to San Siro for Atalanta-Valencia. On buses, cars and train, all of them were heading to a mass gathering where the transmission of the virus would be extremely high. A Spanish journalist who traveled to the game would later become the second infected person in the Valencia region, according to a report by the Associated Press. More than a third of Valencia players would finally test positive for COVID-19 and, as of Friday, Spain has become the European nation with the highest number of cases, on May 2020 beating Italy.

In this article, a two patch SEIR model is considered for the study of the transmission dynamics of MERS-CoV during Hajj and Umrah, as for as Covid-19 during the Bergamo and San Siro Atalanta-Valencia game. The model undergoes a forward bifurcation and, on the other hand, will permit to evaluate optimal control strategies.

2. Formulation of the Model

The model is built on the premise that at the time of the meeting of the two groups there are still no people in quarantine nor is it well known that infection with SARS Cov 2 is in process. For the analysis of the model, some ideas are given in [1]. The dynamics of the transmission of viruses in the population is assumed that there are two groups that interact with each other. Hereafter named Group 1 (for example the local residents of Mecca, or Milan for the Covid-19 case) and Group 2 (for example pilgrims visiting Mecca from other countries or other parts of Saudi Arabia, or the people from Bergamo, Valencia (Spain) visiting Milan for the Atalanta-Valencia soccer game). In this model, some populations are allowed to travel in either direction between the two groups.

The total population of Group 1 at time t, denoted by $N_1\left(t\right)$, is sub-divided into four mutually exclusive compartments of susceptible $S_1\left(t\right)$, exposed (i.e., latently infected, and showing no clinical symptoms of MERS-CoV or Covid-19) $E_1\left(t\right)$, non-quarantined symptomatic (i.e., infected with clinical symptoms of the diseases) $I_1\left(t\right)$ and recovered $R_1\left(t\right)$ individuals. Similarly, the total population of Group 2 at time t, denoted by $N_2\left(t\right)$, is sub-divided into susceptible $S_2\left(t\right)$, exposed $E_2\left(t\right)$, symptomatic $I_2\left(t\right)$ and recovered $R_2\left(t\right)$ individuals. Thus, the total population at time t, $N\left(t\right)$, is given by

$N\left(t\right)=N_1\left(t\right)+N_2\left(t\right),$ (1)

where

$N_1\left(t\right)=S_1\left(t\right)+E_1\left(t\right)+I_1\left(t\right)+R_1\left(t\right)$ (2)

and

$N_2\left(t\right)=S_2\left(t\right)+E_2\left(t\right)+I_2\left(t\right)+R_2\left(t\right).$ (3)

It should be emphasized that, in this paper, individuals in the exposed $E_i$, symptomatic $I_i$ classes $i=1,2$ are infected with the disease, and can transmit it to susceptible individuals. Furthermore, it should be mentioned that the model does not explicitly account for screening.

Table 1 describes the populations and parameters by mean of which the model will be described. The population of susceptible individuals in Group 1, $S_1$, is increased by the recruitment of individuals at a rate ${\pi }_1$, by migration from Group 2 to Group 1 at a rate $m_{21}$. The susceptible individuals in Group 1 acquire the viruses infection, following effective contact with infected individuals in either group (i.e., those in the $E_i$ and $I_i$ classes, where $i=1,2$, at the rate ${\lambda }_1+{\lambda }_2$, where

${\lambda }_i=\frac{{\beta }_i\left({\tau }_e{E}_i+I_i\right)}{S_i+E_i+I_i+R_i}$ (4)

It should be mentioned that the model considers infection transfer during mass gatherings and therefore assumes homogeneous mixing. In 4 for $i=1,2$, ${\beta }_i$ is the effective contact rate for a susceptible in Group i. Furthermore, ${\tau }_e$, with $0<{\tau }_e\le 1$ is the modification parameter accounting for the assumed reduction of infectiousness of exposed individuals in Group i in relation to individuals in the symptomatic classes $I_i$.

Denoting by the derivative with respect to time t, the equations are the following:

It is supposed that population $S_1$ is reduced by migration to Group 2 at a rate $m_{12}$ and by natural death at a rate ${\mu }_1$ ; it is assumed that individuals in all epidemiological compartments suffer natural death at the rate ${\mu }_1$. Thus,

${S^{\prime }}_1={\pi }_1+m_{21}{S}_2-\left({\lambda }_1+{\lambda }_2\right)S_1-m_{12}{S}_1-{\mu }_1{S}_1\mathrm{.}$

The population of exposed individuals in Group 1, $E_1$ is generated by the

Table 1. Description of variables and parameters.

susceptible individuals in Group 1 at the rate ${\lambda }_1+{\lambda }_2$ and by the individuals that migrate from Group 2 at a rate of $a_{21}$. Exposed individual may migrate to Group 2 at a rate of $a_{12}$, develop clinical symptoms of the disease at a rate ${\alpha }_1$ and suffer natural death. Thus,

${E^{\prime }}_1=\left({\lambda }_1+{\lambda }_2\right)S_1+a_{21}{E}_2-a_{12}{E}_1-{\alpha }_1{E}_1-{\mu }_1{E}_1\mathrm{.}$

The population of symptomatic individuals in Group 1, $I_1$ is generated by the exposed population at the rate ${\alpha }_1$ and decreased, natural recovery at a rate ${\kappa }_1$. Hence,

${I^{\prime }}_1={\alpha }_1{E}_1-{\kappa }_1{I}_1-{\mu }_1{I}_1\mathrm{.}$

The population of recovered individuals in Group 1, $R_1$ is generated by the recovery of asymptomatic individuals, infected ${\kappa }_1$, it is further increased by the migration of recovered individuals from Group 2 to Group 1 at a rate $n_{21}$. This population is decreased by migration to Group 2 at a rate $n_{12}$ and natural death. Thus,

${R^{\prime }}_1={\kappa }_1{I}_1+n_{21}{R}_2-n_{12}{R}_1-{\mu }_1{R}_1\mathrm{.}$

The equations for the disease dynamics in Group 2 are derived in a similar fashion. Thus, taking into account the above derivations and assumptions, the equations for the two-group model for the transmission dynamics of the diseases during a mass gathering, is given by

$\begin{array}{l}{S^{\prime }}_1={\pi }_1+m_{21}{S}_2-\left({\lambda }_1+{\lambda }_2\right)S_1-m_{12}{S}_1-{\mu }_1{S}_1,\\ {E^{\prime }}_1=\left({\lambda }_1+{\lambda }_2\right)S_1+a_{21}{E}_2-a_{12}{E}_1-{\alpha }_1{E}_1-{\mu }_1{E}_1,\\ {I^{\prime }}_1={\alpha }_1{E}_1-{\kappa }_1{I}_1-{\mu }_1{I}_1,\\ {R^{\prime }}_1={\kappa }_1{I}_1+n_{21}{R}_2-n_{12}{R}_1-{\mu }_1{R}_1,\end{array}$

$\begin{array}{l}{S^{\prime }}_2={\pi }_2+m_{12}{S}_1-\left({\lambda }_1+{\lambda }_2\right)S_2-m_{21}{S}_2-{\mu }_2{S}_2,\\ {E^{\prime }}_2=\left({\lambda }_1+{\lambda }_2\right)S_2+a_{12}{E}_1-a_{21}{E}_2-{\alpha }_2{E}_2-{\mu }_2{E}_2,\\ {I^{\prime }}_2={\alpha }_2{E}_2-{\kappa }_2{I}_2-{\mu }_2{I}_2,\\ {R^{\prime }}_2={\kappa }_2{I}_2+n_{12}{R}_1-n_{21}{R}_2-{\mu }_2{R}_2.\end{array}$ (5)

where ${\lambda }_i\mathrm{,}\left(i=\mathrm{1,2}\right)$ are given by (4). The model proposed here is in the same line of the models that appear in the literature such as [3] [4] [5], and characterized by:

1) Includes two groups instead of a single group/patch as is considered in [3] [4].

2) Allows for disease transmission by latently-infected individuals $E_i$ which is is not considered in [5].

3) Including heterogeneity in the transmission based on groups. That is, susceptible individuals can acquire infection from infected individuals in Group 1 at rate ${\lambda }_1$ or Group 2 at rate ${\lambda }_2$, defined in (4), so the infection rate is given as ( ${\lambda }_1+{\lambda }_2$ ), not considered in [3] [4] [5].

4) The population of the infected asymptomatic individuals is considered here as included in the exposed group.

The flow diagram of the model is represented in Figure 1.

Basic Properties

The parameters of Equation (5) are supposed to be non-negative and it is well known that the model is well behaved, i.e., solutions with positive initial conditions stay positive for all $t>0$ and can be stated following the methods of [6].

Theorem 1. All of the solutions of the model giving by Equation (5) with positive initial conditions remain positive for all time $t>0$.

Let $\mathcal{D}$ be the following set,

$\begin{array}{l}\mathcal{D}=\{\left(S_1,E_1,I_1,R_1,S_2,E_2,I_2,R_2\right)\in {ℝ}^8:\begin{array}{c}\\ \end{array}\\ S_1+E_1+I_1+R_1+S_2+E_2+I_2+R_2<\frac{{\pi }_1+{\pi }_2}{\min \left\{{\mu }_1,{\mu }_2\right\}}\}.\end{array}$

For the region defined with the set $\mathcal{D}$ the following lemma holds.

Lemma 2. The region described by the set $\mathcal{D}$ is positively invariant for the system given by Equation (5).

Figure 1. Flow diagram of the model (5).

Proof. By mean of the notations given in (1), (2), (3) and afterwards of adding Equation (5), the following equality can be obtained

$\frac{\text{d}N\left(t\right)}{\text{d}t}={\pi }_1+{\pi }_2-{\mu }_1{N}_1\left(t\right)-{\mu }_2{N}_2(\; t\; )$

and therefore,

$\frac{\text{d}N\left(t\right)}{\text{d}t}\le {\pi }_1+{\pi }_2-\min \left\{{\mu }_1,{\mu }_2\right\}N\left(t\right).$

This implies that for $N\left(t\right)\ge \frac{{\pi }_1+{\pi }_2}{\min \left\{{\mu }_1,{\mu }_2\right\}}$, $\frac{\text{d}N\left(t\right)}{\text{d}t}\le 0$. Hence by standard comparison arguments [7] [8] the following inequality holds

$N\left(t\right)\le N\left(0\right){\text{e}}^{-\min \left\{{\mu }_1,{\mu }_2\right\}t}+\frac{{\pi }_1+{\pi }_2}{\min \left\{{\mu }_1,{\mu }_2\right\}}\left(1-{\text{e}}^{-\min \left\{{\mu }_1,{\mu }_2\right\}t}\right).$

Therefore the region $\mathcal{D}$ is positively invariant which means that solutions that enter in $\mathcal{D}$ remain there for all subsequent time t. QED

It is well known that systems that have the property as in the previous lemma are said to be mathematically and epidemiologically well-posed [9].

3. Existence and Asymptotic Stability of Equilibria

In this section the equilibria of the system (5) will be studied.

3.1. Local Asymptotic Stability of Disease-Free Equilibrium

For the system given by Equation (5) the infected compartments are

$E_1\mathrm{,}{E}_2\mathrm{,}{I}_1\mathrm{,}{I}_2\mathrm{.}$

The equilibrium solution for

$E_1^{*}=E_2^{*}=I_1^{*}=I_2^{*}=R_1^{*}=R_2^{*}=0.$

is given by

${\mathcal{E}}_0=\left(S_1^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,0,0,0,}{S}_2^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,0,0,0}\right)\mathrm{,}$ (6)

where

$S_1^{\mathrm{<mo>\; *\; </mo>}}=\frac{\left(m_{21}+{\mu }_2\right){\pi }_1+m_{21}{\pi }_2}{m_{12}{\mu }_2+m_{21}{\mu }_1+{\mu }_1{\mu }_2}\mathrm{,}$ (7)

and

$S_2^{\mathrm{<mo>\; *\; </mo>}}=\frac{\left(m_{12}+{\mu }_1\right){\pi }_2+m_{12}{\pi }_1}{m_{12}{\mu }_2+m_{21}{\mu }_1+{\mu }_1{\mu }_2}$ (8)

is the disease-free equilibrium for the system. In order to establish the local-asymptotic stability of the DFE equilibrium (6), the next operator theory, whose methods and notations are given in [10], will be used here. In the case of the model investigated here, the non-negative, $\mathcal{F}$, and the M-matrix, $\mathcal{V}$, are given respectively by

$\mathcal{F}=\left(\begin{array}{cccc}{\beta }_1{\tau }_e& {\beta }_1& {\beta }_2{\tau }_e{S}^{\ast }& {\beta }_2{S}^{\ast }\\ 0& 0& 0& 0\\ \frac{{\beta }_1{\tau }_e}{S^{\ast }}& \frac{{\beta }_1}{S^{\ast }}& {\beta }_2{\tau }_e& {\beta }_2\\ 0& 0& 0& 0\end{array}\right)$

and

$\mathcal{V}=\left(\begin{array}{cccc}{a}_1& 0& -a_{21}& 0\\ -{\alpha }_1& s_1& 0& 0\\ -a_{12}& 0& a_2& 0\\ 0& 0& -{\alpha }_2& s_2\end{array}\right)\mathrm{,}$

where

$S^{*}=\frac{S_1^{*}}{S_2^{*}},s_i={\kappa }_i+{\mu }_i,\left(i=1,2\right),$

$a_1=a_{12}+{\alpha }_1+{\mu }_1,a_2=a_{21}+{\alpha }_2+{\mu }_2$

and putting

$a=a_1{a}_2-a_{12}{a}_{21}>0$

the inverse matrix of $\mathcal{V}$ is given by

${\mathcal{V}}^{-1}=\left(\begin{array}{cccc}-\frac{a_2}{a}& 0& -\frac{a_{21}}{a}& 0\\ -a_2\frac{{\alpha }_1}{a{s}_1}& \frac{1}{s_1}& -a_{21}\frac{{\alpha }_1}{a{s}_1}& 0\\ -\frac{a_{12}}{a}& 0& -\frac{a_1}{a}& 0\\ -a_{12}\frac{{\alpha }_2}{a{s}_2}& 0& -a_1\frac{{\alpha }_2}{a{s}_2}& \frac{1}{s_2}\end{array}\right)\mathrm{.}$

Hence the matrix $\mathcal{F}{\mathcal{V}}^{-1}$ is given by

$\mathcal{F}{\mathcal{V}}^{-1}=\left(\begin{array}{cccc}-u& \frac{1}{s_1}{\beta }_1& -v& \frac{1}{s_2}{S}^{\ast }{\beta }_2\\ 0& 0& 0& 0\\ -w& \frac{1}{s_1{S}^{\ast }}{\beta }_1& -r& \frac{1}{s_2}{\beta }_2\\ 0& 0& 0& 0\end{array}\right)$

where, $u\mathrm{,}v\mathrm{,}w\mathrm{,}r$ are defined by the following relations

$au=a_2{\beta }_1\left(c_1+{\tau }_e\right)+a_{12}{S}^{\ast }{\beta }_2\left(c_2+{\tau }_e\right)\mathrm{,}$

$av=a_{21}{\beta }_1\left(c_1+{\tau }_e\right)+a_1{S}^{\ast }{\beta }_2\left(c_2+{\tau }_e\right)\mathrm{,}$

$aw=a_{12}{\beta }_2\left(c_2+{\tau }_e\right)+\frac{a_2}{S^{\ast }}{\beta }_1\left(c_1+{\tau }_e\right)\mathrm{,}$

$ar=a_1{\beta }_2\left(c_2+{\tau }_e\right)+\frac{a_{21}}{S^{\ast }}{\beta }_1\left(c_1+{\tau }_e\right)\mathrm{.}$

As is established in [10], the reproduction number ${\mathcal{R}}_v$ of the system is given by the spectral radius of the next generation matrix $\mathcal{F}{\mathcal{V}}^{-1}$. The eigenvalues of this matrix are given by

$\lambda =0,\lambda =-\frac{1}{2}r-\frac{1}{2}u-\frac{1}{2}\sqrt{4vw+{\left(r-u\right)}^2},\lambda =-\frac{1}{2}r-\frac{1}{2}u+\frac{1}{2}\sqrt{4vw+{\left(r-u\right)}^2}.$

Clearly, the spectral radius of the matrix $\mathcal{F}{\mathcal{V}}^{-1}$ is therefore given by

${\mathcal{R}}_v=\frac{1}{2}r+\frac{1}{2}u+\frac{1}{2}\sqrt{4vw+{\left(r-u\right)}^2}.$

It can be seen that the following relations hold:

Lemma 3.

1) The following equality holds, $vw=ru$.

2) The reproduction number, ${\mathcal{R}}_v$, is given by ${\mathcal{R}}_v=r+u$.

3) As a consequence of the previous relation, ${\mathcal{R}}_v\le 1$ (respectively ${\mathcal{R}}_v>1$ ) is equivalent to $r+u\le 1$ (respectively $r+u>1$ ).

Then according with the theory of the next generation operator, see [10], the following theorem holds

Theorem 4. The DFE equilibrium (6) of the model given by the system (5) is locally-asymptotically stable for ${\mathcal{R}}_v<1$ , and unstable for ${\mathcal{R}}_v>1$.

From a epidemiological point of view the previous theorem says that a small influx of infected individuals does not produce a large Coronavirus outbreak in a mass gatherings for values of ${\mathcal{R}}_v<1$ and therefore, the disease is subject to control whenever the sizes of the initial sub-populations of the model belong to the attraction basin of the equilibrium DFE (6).

3.2. Existence of the Endemic Equilibrium Point

By mean of

${\mathcal{E}}_1=\left(S_1^{\mathrm{<mo>\; *\; </mo>\; <mo>\; *\; </mo>}}\mathrm{,}{E}_1^{\mathrm{<mo>\; *\; </mo>\; <mo>\; *\; </mo>}}\mathrm{,}{I}_1^{\mathrm{<mo>\; *\; </mo>\; <mo>\; *\; </mo>}}\mathrm{,}{R}_1^{\mathrm{<mo>\; *\; </mo>\; <mo>\; *\; </mo>}}\mathrm{,}{S}_2^{\mathrm{<mo>\; *\; </mo>\; <mo>\; *\; </mo>}}\mathrm{,}{E}_2^{\mathrm{<mo>\; *\; </mo>\; <mo>\; *\; </mo>}}\mathrm{,}{I}_2^{\mathrm{<mo>\; *\; </mo>\; <mo>\; *\; </mo>}}\mathrm{,}{R}_2^{\mathrm{<mo>\; *\; </mo>\; <mo>\; *\; </mo>}}\right)$

represents an endemic equilibrium point of the system (5), where all of the infected components are non-zero. In such a case the infection forces are given by

${\lambda }_i^{**}=\frac{{\beta }_i\left({\tau }_e{E}_i^{**}+I_i^{**}\right)}{S_i^{**}+E_i^{**}+I_i^{**}+R_i^{**}},i=1,2.$ (9)

In order to obtain the coordinates of ${\mathcal{E}}_1$ consider the equations for ${S^{\prime }}_1$ and ${S^{\prime }}_2$ of system (5), for the coordinates of ${\mathcal{E}}_1$ the following equations hold,

$\left({\lambda }^{**}+m_1\right)S_1-m_{21}{S}_2={\pi }_1,$

$-m_{12}{S}_1+\left({\lambda }^{\mathrm{<mo>\; *\; </mo>\; <mo>\; *\; </mo>}}+m_2\right)S_2={\pi }_2\mathrm{,}$

where

$m_1=m_{12}+{\mu }_1,m_2=m_{21}+{\mu }_2,$

and

${\lambda }^{**}={\lambda }_1^{**}+{\lambda }_2^{**}.$ (10)

Thus, the values of $S_1^{\mathrm{<mo>\; *\; </mo>\; <mo>\; *\; </mo>}}\mathrm{,}{S}_2^{\mathrm{<mo>\; *\; </mo>\; <mo>\; *\; </mo>}}$ are given by

$S_1^{\mathrm{<mo>\; *\; </mo>\; <mo>\; *\; </mo>}}=\frac{\left({\lambda }^{\mathrm{<mo>\; *\; </mo>\; <mo>\; *\; </mo>}}+m_2\right){\pi }_1+m_{21}{\pi }_2}{{\left({\lambda }^{\mathrm{<mo>\; *\; </mo>\; <mo>\; *\; </mo>}}\right)}^2+\left(m_1+m_2\right)\left({\lambda }^{\mathrm{<mo>\; *\; </mo>\; <mo>\; *\; </mo>}}\right)+m}\mathrm{,}$

$S_2^{\mathrm{<mo>\; *\; </mo>\; <mo>\; *\; </mo>}}=\frac{\left({\lambda }^{\mathrm{<mo>\; *\; </mo>\; <mo>\; *\; </mo>}}+m_1\right){\pi }_2+m_{12}{\pi }_1}{{\left({\lambda }^{\mathrm{<mo>\; *\; </mo>\; <mo>\; *\; </mo>}}\right)}^2+\left(m_1+m_2\right)\left({\lambda }^{\mathrm{<mo>\; *\; </mo>\; <mo>\; *\; </mo>}}\right)+m}\mathrm{,}$

where $m=m_1{m}_2-m_{12}{m}_{21}>0$. In the same way considering the equations for ${E^{\prime }}_1$ and ${E^{\prime }}_2$ of system (5) gives for $E_1^{\mathrm{<mo>\; *\; </mo>\; <mo>\; *\; </mo>}}\mathrm{,}{E}_2^{\mathrm{<mo>\; *\; </mo>\; <mo>\; *\; </mo>}}$ the following values

$E_1^{**}=\frac{a_2{\lambda }^{**}}{a}{S}_1^{**}+\frac{a_{21}{\lambda }^{**}}{a}{S}_2^{**},$

$E_2^{**}=\frac{a_{12}{\lambda }^{**}}{a}{S}_1^{**}+\frac{a_1{\lambda }^{**}}{a}{S}_2^{**}.$

Finally from the equations for ${I^{\prime }}_1$, ${I^{\prime }}_2$, ${R^{\prime }}_1$ and ${R^{\prime }}_2$ of system (5) gives the following values

$I_1^{**}=\frac{{\alpha }_1{E}_1^{**}}{s_1}=\frac{{\alpha }_1}{s_1}\left(\frac{a_2{\lambda }^{**}}{a}{S}_1^{**}+\frac{a_{21}{\lambda }^{**}}{a}{S}_2^{**}\right),$

$I_2^{**}=\frac{{\alpha }_2{E}_2^{**}}{s_2}=\frac{{\alpha }_2}{s_2}\left(\frac{a_{12}{\lambda }^{**}}{a}{S}_1^{**}+\frac{a_1{\lambda }^{**}}{a}{S}_2^{**}\right),$

$R_1^{**}=\frac{{\lambda }^{**}}{a}{d}_{12}{S}_1^{**}+\frac{{\lambda }^{**}}{a}{d}_{21}{S}_2^{**},$

$R_2^{**}=\frac{{\lambda }^{**}}{a}{f}_{12}{S}_1^{**}+\frac{{\lambda }^{**}}{a}{f}_{21}{S}_2^{**},$

where, $c_i=\frac{{\alpha }_i}{s_i},i=1,2$,

$d_{12}=\frac{{\kappa }_1{n}_2{c}_1{a}_2}{n}+\frac{{\kappa }_2{n}_{21}{c}_2{a}_{12}}{n},d_{21}=\frac{{\kappa }_1{n}_2{c}_1{a}_{21}}{n}+\frac{{\kappa }_2{n}_{21}{c}_2{a}_1}{n},$

$f_{12}=\frac{{\kappa }_1{n}_{12}{c}_1{a}_2}{n}+\frac{{\kappa }_2{n}_1{c}_2{a}_{12}}{n},f_{21}=\frac{{\kappa }_1{n}_{12}{c}_1{a}_{21}}{n}+\frac{{\kappa }_2{n}_1{c}_2{a}_1}{n},$

and

$n_1=n_{12}+{\mu }_1,n_2=n_{21}+{\mu }_2,n=n_1{n}_2-n_{12}{n}_{21}.$

Now from the relation (10) substituting all of the previous values in terms of ${\lambda }^{\mathrm{<mo>\; *\; </mo>\; <mo>\; *\; </mo>}}$ gives

${\lambda }_1^{**}=\frac{{\beta }_1\left(c_1+{\tau }_e\right){\lambda }^{**}\left(w_{21}{\lambda }^{**}+r_1\right)}{a{v}_{21}+M_{1d}{\lambda }^{**}+M_{2d}{\left({\lambda }^{**}\right)}^2}$

and

${\lambda }_2^{**}=\frac{{\beta }_2\left(c_2+{\tau }_e\right){\lambda }^{**}\left(w_{12}{\lambda }^{**}+r_2\right)}{a{v}_{12}+M_{1f}{\lambda }^{**}+M_{2f}{\left({\lambda }^{**}\right)}^2},$

where

$v_{12}=m_1{\pi }_2+m_{12}{\pi }_1,v_{21}=m_2{\pi }_1+m_{21}{\pi }_2,$

$w_{12}=a_1{\pi }_2+a_{12}{\pi }_1,w_{21}=a_2{\pi }_1+a_{21}{\pi }_2,r_1=a_1{v}_{12}+a_{12}{v}_{21},r_2=a_2{v}_{21}+a_{21}{v}_{12},$

$M_{1d}=r_1+a{\pi }_1+c_1{r}_1+d_{12}{v}_{12}+d_{21}{v}_{21},M_{2d}=w_{21}+c_1{w}_{21}+d_{21}{\pi }_1+d_{12}{\pi }_2,$

$M_{1f}=r_2+a{\pi }_2+c_2{r}_2+f_{12}{v}_{12}+f_{21}{v}_{21},M_{2f}=w_{12}+c_2{w}_{12}+f_{21}{\pi }_1+f_{12}{\pi }_2.$

Substituting all of these values in terms of ${\lambda }^{\mathrm{<mo>\; *\; </mo>\; <mo>\; *\; </mo>}}$ in Equation (4) gives either ${\lambda }^{**}=0$ (which implies that ${\lambda }_1^{**}=0$ and ${\lambda }_2^{**}=0$ that corresponds to the DFE case), or

$A_4{\left({\lambda }^{\mathrm{<mo>\; *\; </mo>\; <mo>\; *\; </mo>}}\right)}^4+A_3{\left({\lambda }^{\mathrm{<mo>\; *\; </mo>\; <mo>\; *\; </mo>}}\right)}^3+A_2{\left({\lambda }^{\mathrm{<mo>\; *\; </mo>\; <mo>\; *\; </mo>}}\right)}^2+A_1\left({\lambda }^{\mathrm{<mo>\; *\; </mo>\; <mo>\; *\; </mo>}}\right)+A_0=0.$ (11)

Here,

$A_0=a^2{v}_{12}{v}_{21}-a{\beta }_1\left(c_1+{\tau }_e\right)r_1{v}_{12}-a{\beta }_2\left(c_2+{\tau }_e\right)r_2{v}_{21}\mathrm{,}$

$\begin{array}{c}{A}_1=\left(a{v}_{21}-{\beta }_1\left(c_1+{\tau }_e\right)r_1\right)M_{1f}+\left(a{v}_{12}-{\beta }_2\left(c_2+{\tau }_e\right)r_2\right)M_{1d}\\ -a\left({\beta }_1\left(c_1+{\tau }_e\right)w_{21}{v}_{12}+{\beta }_2\left(c_2+{\tau }_e\right)w_{12}{v}_{21}\right)\mathrm{,}\end{array}$

$\begin{array}{c}{A}_2=\left(a{v}_{21}-{\beta }_1\left(c_1+{\tau }_e\right)r_1\right)M_{2f}+\left(a{v}_{12}-{\beta }_2\left(c_2+{\tau }_e\right)r_2\right)M_{2d}\\ -\left({\beta }_1\left(c_1+{\tau }_e\right)w_{21}{M}_{1f}+{\beta }_2\left(c_2+{\tau }_e\right)w_{12}{M}_{1d}-M_{1d}{M}_{1f}\right)\mathrm{,}\end{array}$

$A_3=\left(M_{1d}-{\beta }_1\left(c_1+{\tau }_e\right)w_{21}\right)M_{2f}+\left(M_{1f}-{\beta }_2\left(c_2+{\tau }_e\right)w_{12}\right)M_{2d},$

$A_4=M_{2d}{M}_{2f}$

Note that $A_4$ is always positive. The other coefficients may change of sign. The following lemma implies the negativity of $A_0$.

Lemma 5. Suppose that $u+r>1$. If the following condition holds

$a<{\beta }_1{a}_2\left(c_1+{\tau }_e\right)+{\beta }_2{a}_1\left(c_2+{\tau }_e\right),$ (12)

then $A_0<0$.

Proof. This proof need long calculations, the idea is to note that it is possible to write as,

$A_0=\left(u+r-1\right)P\left(v_{12},v_{21}\right),$

where the hypothesis (12) implies that the polynomial

$P\left(v_{12}\mathrm{,}{v}_{21}\right)=b{v}_{12}^2+c{v}_{21}^2+d{v}_{12}{v}_{21}$

has negative coefficients $b\mathrm{,}c\mathrm{,}d$. And therefore $A_0<0$. QED

The previous lemmas implies that if ${\mathcal{R}}_v>1$ and (12) hold, then $A_0<0$. The possibilities for the sign of the coefficients of the coefficients and positive roots of the polynomial (11) are given by Table 2 with $A_4>0$ and $A_0<0$ where the conclusions are in accord with the Descartes criteria for the existence of positive real roots of polynomials.

With all of these observations subsists the following result.

Theorem 6. The model given by the system (5)[i.]

1) has a unique EEP if ${\mathcal{R}}_v>1$, the condition given in (12) and any of the Cases 1, 2, 3, 4 of Table 2 hold;

2) has one or multiple EEP if the conditions given in (12) and any of the Cases 5, 6, 7, 8 of Table 2 hold;

3) may have more than one EEP if the conditions given in (12) and any of the Cases 2 to 8 of Table 3 hold.

From the epidemiological point of view this result implies that the viruses disease will establish itself in the population, when ${\mathcal{R}}_v>1$ whenever the initial states of the sub-populations of the system (5) belong to the basin of attraction of the unique endemic equilibrium ${\mathcal{E}}_1$. In the case when ${\mathcal{R}}_v>1$ and the endemic equilibrium of the system (5) is stable.

3.3. Stability Analysis

When ${\mathcal{R}}_v<1$, the equilibrium ${\mathcal{E}}_0$ is asymptotically stable and for values of the populations in the attraction basin of this point, the control and elimination of the disease is guaranteed. For ${\mathcal{R}}_v>1$ the equilibrium ${\mathcal{E}}_0$ lost its stability and

Table 2. Table of change of sign and positive roots for $A_0<0$.

Table 3. Table of change of sign and positive roots for $A_0>0$.

a forward bifurcation occurs. That means that an equilibrium of positive coordinates ${\mathcal{E}}_1$ wins stability what drives the populations to an endemic equilibrium point situation, where all of the coordinates are positive. The forward bifurcation term indicates too that the model cannot undergoes a backward bifurcation where ${\mathcal{E}}_0$ and ${\mathcal{E}}_1$ can coexist in a stable way. The following theorem is a consequence of the Centre Manifold theory as appear in [11], and the Theorem 4.1 of [12].

Theorem 7. The model described for mean of the system (5) exhibits a forward bifurcation at ${\mathcal{R}}_v=1$ in the sense previously described.

Proof. It is convenient to make the following change of variables:

$\begin{array}{l}{x}_1=S_1,x_2=E_1,x_3=I_1,x_4=R_1\\ x_5=S_2,x_6=E_2,x_7=I_2,x_8=R_2.\end{array}$

Let $X=col\left(x_1,x_2,\cdots ,x_7,x_8\right)$. Then the model (5) can be re-written in vector shape doing $F\left(X\right)=col\left(f_1\left(X\right)\mathrm{,}{f}_2\left(X\right)\mathrm{,}\cdots \mathrm{,}{f}_7\left(X\right)\mathrm{,}{f}_8\left(X\right)\right)$, as follows:

$f_1\left(X\right)={\pi }_1+m_{21}{x}_5-\frac{{\beta }_1\left({\tau }_e{x}_2+x_3\right)x_1}{x_1+x_2+x_3+x_4}-\frac{{\beta }_2\left({\tau }_e{x}_6+x_7\right)x_1}{x_5+x_6+x_7+x_8}-m_{12}{x}_1-{\mu }_1{x}_1,$

$f_2\left(X\right)=\frac{{\beta }_1\left({\tau }_e{x}_2+x_3\right)x_1}{x_1+x_2+x_3+x_4}+\frac{{\beta }_2\left({\tau }_e{x}_6+x_7\right)x_1}{x_5+x_6+x_7+x_8}+a_{21}{x}_6-a_{12}{x}_2-{\alpha }_1{x}_2-{\mu }_1{x}_2,$

$f_3\left(X\right)={\alpha }_1{x}_2-{\kappa }_1{x}_3-{\mu }_1{x}_3,$

$f_4\left(X\right)={\kappa }_1{x}_3+n_{21}{x}_8-n_{12}{x}_4-{\mu }_1{x}_4,$

$f_5\left(X\right)={\pi }_1+m_{12}{x}_1-\frac{{\beta }_1\left({\tau }_e{x}_2+x_3\right)x_5}{x_1+x_2+x_3+x_4}-\frac{{\beta }_2\left({\tau }_e{x}_6+x_7\right)x_5}{x_5+x_6+x_7+x_8}-m_{21}{x}_5-{\mu }_2{x}_5,$

$f_6\left(X\right)=\frac{{\beta }_1\left({\tau }_e{x}_2+x_3\right)x_5}{x_1+x_2+x_3+x_4}+\frac{{\beta }_2\left({\tau }_e{x}_6+x_7\right)x_5}{x_5+x_6+x_7+x_8}+a_{12}{x}_2-a_{21}{x}_6-{\alpha }_1{x}_6-{\mu }_2{x}_6,$

$f_7\left(X\right)={\alpha }_1{x}_6-{\kappa }_1{x}_7-{\mu }_1{x}_7,$

$f_8\left(X\right)={\kappa }_2{x}_7+n_{12}{x}_4-n_{21}{x}_8-{\mu }_2{x}_8.$

Therefore system (5) has been converted in the following

$\frac{\text{d}{x}_1}{\text{d}t}={\pi }_1+m_{21}{x}_5-\frac{{\beta }_1\left({\tau }_e{x}_2+x_3\right)x_1}{x_1+x_2+x_3+x_4}-\frac{{\beta }_2\left({\tau }_e{x}_6+x_7\right)x_1}{x_5+x_6+x_7+x_8}-m_1{x}_1,$

$\frac{\text{d}{x}_2}{\text{d}t}=\frac{{\beta }_1\left({\tau }_e{x}_2+x_3\right)x_1}{x_1+x_2+x_3+x_4}+\frac{{\beta }_2\left({\tau }_e{x}_6+x_7\right)x_1}{x_5+x_6+x_7+x_8}+a_{21}{x}_6-a_1{x}_2\mathrm{,}$

$\frac{\text{d}{x}_3}{\text{d}t}={\alpha }_1{x}_2-s_1{x}_3\mathrm{,}$

$\frac{\text{d}{x}_4}{\text{d}t}={\kappa }_1{x}_3+n_{21}{x}_8-n_1{x}_4\mathrm{,}$ (13)

$\frac{\text{d}{x}_5}{\text{d}t}={\pi }_2+m_{12}{x}_1-\frac{{\beta }_1\left({\tau }_e{x}_2+x_3\right)x_5}{x_1+x_2+x_3+x_4}-\frac{{\beta }_2\left({\tau }_e{x}_6+x_7\right)x_5}{x_5+x_6+x_7+x_8}-m_2{x}_5\mathrm{,}$

$\frac{\text{d}{x}_6}{\text{d}t}=\frac{{\beta }_1\left({\tau }_e{x}_2+x_3\right)x_5}{x_1+x_2+x_3+x_4}+\frac{{\beta }_2\left({\tau }_e{x}_6+x_7\right)x_5}{x_5+x_6+x_7+x_8}+a_{12}{x}_2-a_2{x}_6\mathrm{,}$

$\frac{\text{d}{x}_7}{\text{d}t}={\alpha }_1{x}_6-s_2{x}_7\mathrm{,}$

$\frac{\text{d}{x}_8}{\text{d}t}={\kappa }_2{x}_7+n_{12}{x}_4-n_2{x}_8\mathrm{.}$

where $a_i,s_i,m_i,n_i,i=1,2$, are the same given in the previous section. With these notations the model given by the system (5) may be thought as

$X^{\prime }=F\left(X\right)\mathrm{.}$ (14)

The Jacobian matrix, $JF\left({\mathcal{E}}_0\right)$ of the function $F\left(X\right)$ evaluated at the DFE critical point, ${\mathcal{E}}_0$ is given by

$JF\left({\mathcal{E}}_0\right)=\left(\begin{array}{cc}{A}_{7\times 7}& B_{7\times 7}\\ C_{7\times 7}& D_{7\times 7}\end{array}\right)$ (15)

where,

$A=\left(\begin{array}{cccc}-m_1& -{\beta }_1{\tau }_e& -{\beta }_1& 0\\ 0& {\beta }_1{\tau }_e-a_1& {\beta }_1& 0\\ 0& {\alpha }_1& -s_1& 0\\ 0& 0& {\kappa }_1& -n_1\end{array}\right)$

$B=\left(\begin{array}{cccc}{m}_{21}& -{\beta }_2{\tau }_e{S}^{\mathrm{<mo>\; *\; </mo>}}& -{\beta }_2{S}^{\mathrm{<mo>\; *\; </mo>}}& 0\\ 0& {\beta }_2{\tau }_e{S}^{\mathrm{<mo>\; *\; </mo>}}+a_{21}& {\beta }_2{S}^{\mathrm{<mo>\; *\; </mo>}}& 0\\ 0& 0& 0& 0\\ 0& 0& 0& n_{21}\end{array}\right)$

$C=\left(\begin{array}{cccc}{m}_{12}& -{\beta }_1{\tau }_e{\left(S^{\mathrm{<mo>\; *\; </mo>}}\right)}^{-1}& -{\beta }_1{\left(S^{\mathrm{<mo>\; *\; </mo>}}\right)}^{-1}& 0\\ 0& {\beta }_1{\tau }_e{\left(S^{\mathrm{<mo>\; *\; </mo>}}\right)}^{-1}+a_{12}& {\beta }_1{\left(S^{\mathrm{<mo>\; *\; </mo>}}\right)}^{-1}& 0\\ 0& 0& 0& 0\\ 0& 0& 0& n_{12}\end{array}\right)$

and,

$D=\left(\begin{array}{cccc}-m_2& -{\beta }_2{\tau }_e& -{\beta }_2& 0\\ 0& {\beta }_2{\tau }_e-a_2& {\beta }_2& 0\\ 0& {\alpha }_2& -s_2& 0\\ 0& 0& {\kappa }_2& -n_2\end{array}\right)\mathrm{.}$

Let ${\mathcal{R}}_v=1$ and suppose that ${\beta }_1$ is chosen as a bifurcation parameter. The transformed system (14), with ${\beta }_1={\beta }_1^{\mathrm{<mo>\; *\; </mo>}}$, has a hyperbolic equilibrium point (i.e., the corresponding associated linear system has a simple eigenvalue with zero real part, and all other eigenvalues have negative real part), so that the center manifold theory [12] can be used to analyze the dynamics of system (14) near ${\beta }_1={\beta }_1^{*}$ It can be shown that the right eigenvector of ${JF\left({\mathcal{E}}_0\right)|}_{{\beta }_1={\beta }_1^{\mathrm{<mo>\; *\; </mo>}}}$, denoted by $w$, is given by $w=col\left(w_1\mathrm{,}{w}_2\mathrm{,}\cdots \mathrm{,}{w}_7\mathrm{,}{w}_8\right)$, with, (it should be noted that the expressions for $w_1\mathrm{,}\cdots \mathrm{,}{w}_8$ are very lengthy and they are not reported here)

$\begin{array}{l}{w}_1<0,w_2>0,w_3>0,w_4>0,w_5<0,w_6>0,w_7>0,w_8>0.\hfill \end{array}$

Similarly, ${JF\left({\mathcal{E}}_0\right)|}_{{\beta }_1={\beta }_1^{\mathrm{<mo>\; *\; </mo>}}}$ has a left eigenvector, $v$, given by $v=\left(v_1\mathrm{,}{v}_2\mathrm{,}\cdots \mathrm{,}{v}_7\mathrm{,}{v}_8\right)$, with

$\begin{array}{l}{v}_1=0,v_2>0,v_3>0,v_4=0,v_5=0,v_6>0,v_7>0,v_8=0.\hfill \end{array}$

The computation of the required bifurcation coefficients a and b gives the following results according with Theorem 4.1 of [12].

$\begin{array}{c}a={\displaystyle \underset{k,i,j=1}{\overset{14}{\sum }}}{v}_k{w}_i{w}_j\frac{{\partial }^2{f}_k\left(0,0\right)}{\partial x_i\partial x_j}\\ =-\frac{4{\beta }_1\left({\tau }_e{w}_2+w_3\right)v_2}{{\left(S_1^{*}\right)}^2}\left(-w_3{S}_1^{*}\right)\\ -\frac{2{\beta }_1\left({\tau }_e{w}_2+w_3\right)v_6}{{\left(S_1^{*}\right)}^2}\left[S_2^{*}\left(w_2+w_3+w_4\right)-w_5\left(S_1^{*}\right)\right]\\ -\frac{2{\beta }_2\left({\tau }_e{w}_6+w_7\right)v_2}{{\left(S_2^{*}\right)}^2}\left[S_1^{*}\left(w_6+w_7+w_8\right)-w_1\left(S_2^{*}\right)\right]\end{array}$

$\begin{array}{l}-\frac{2{\beta }_2\left({\tau }_e{w}_6+w_7+\right)v_2}{{\left(S_2^{*}\right)}^2}\left[-w_3\left(S_2^{*}\right)\right]\\ -\frac{2{\beta }_2\left({\tau }_e{w}_6+w_7\right)v_6}{{\left(S_2^{*}\right)}^2}\left[S_2^{*}\left(w_6+w_7+w_8\right)\right]\\ <0,\end{array}$

and

$b={\displaystyle \underset{k,i=1}{\overset{14}{\sum }}}{v}_k{w}_i\frac{{\partial }^2{f}_k\left(0,0\right)}{\partial x_i\partial {\beta }_1^{*}}=\frac{\left(v_2{S}_1^{*}\right)\left({\tau }_e{w}_2+w_3\right)}{S_1^{*}}>0.$

The sign of the parameters $a<0,b>0$. When the bifurcation parameter changes from negative to positive, the equilibrium point ${\mathcal{E}}_0$ changes its stability from stable to unstable. Correspondingly a negative unstable equilibrium becomes positive and locally asymptotically stable. QED

4. Optimal Control

In this section an optimal control based in the exposed strategies considered in the system (5) will be carried out. The system will be reformulated as an optimal control problem in the following way,

$\begin{array}{l}\frac{\text{d}{S}_1}{\text{d}t}={\pi }_1+m_{21}{S}_2-\left({\lambda }_1+{\lambda }_2\right)S_1-u_1\left(t\right)S_1-m_{12}{S}_1-{\mu }_1{S}_1,\\ \frac{\text{d}{E}_1}{\text{d}t}=\left({\lambda }_1+{\lambda }_2\right)S_1+u_1\left(t\right)S_1+a_{21}{E}_2-a_{12}{E}_1-{\alpha }_1{E}_1-{\mu }_1{E}_1,\\ \frac{\text{d}{I}_1}{\text{d}t}={\alpha }_1{E}_1-{\kappa }_1{I}_1-{\mu }_1{I}_1,\\ \frac{\text{d}{R}_1}{\text{d}t}={\kappa }_1{I}_1+n_{21}{R}_2-n_{12}{R}_1-{\mu }_1{R}_1,\end{array}$

$\begin{array}{l}\frac{\text{d}{S}_2}{\text{d}t}={\pi }_2+m_{12}{S}_1-\left({\lambda }_1+{\lambda }_2\right)S_2-u_2\left(t\right)S_2-m_{21}{S}_2-{\mu }_2{S}_2,\\ \frac{\text{d}{E}_2}{\text{d}t}=\left({\lambda }_1+{\lambda }_2\right)S_2+u_2\left(t\right)S_2+a_{12}{E}_1-a_{21}{E}_2-{\alpha }_2{E}_2-{\mu }_2{E}_2,\\ \frac{\text{d}{I}_2}{\text{d}t}={\alpha }_2{E}_2-{\kappa }_2{I}_2-{\mu }_2{I}_2,\\ \frac{\text{d}{R}_2}{\text{d}t}={\kappa }_2{I}_2+n_{12}{R}_1-n_{21}{R}_2-{\mu }_2{R}_2,\end{array}$ (16)

where ${\lambda }_i,i=1,2$ are given by (4) and the functions $u_i\left(t\right),i=1,2$ are bounded and Lebesgue integrable. The optimal control problem looks for minimizing the number of exposed individuals. The problem reduce to minimize the functional

$J\left(u_1\left(t\right),u_2\left(t\right)\right)={\displaystyle {\int }_0^T}\left[d_1{E}_1+X_1{u}_1^2\left(t\right)+{\bar{d}}_1{E}_2+X_2{u}_2^2\left(t\right)\right]\text{d}t$

subject to the differential Equation (16), the coefficients $d_i,d_i,X_i,i=1,2$ are balancing cost factors. Specifically the problem consists of to find optimal control functions $\left(u_1^{†}\left(t\right)\mathrm{,}{u}_2^{†}\right)$ satisfying the Lagrange formulation (see [13] )

$J\left(u_1^{†}\left(t\right)\mathrm{,}{u}_2^{†}\left(t\right)\right)=\underset{U\mathrm{<mo>\; :\; </mo>}{x}^u\left(t_0\right)=x_0}{\inf }J\left(u_1\left(t\right)\mathrm{,}{u}_2\left(t\right)\right)\mathrm{,}$ (17)

$U=\left\{\left(u_1\mathrm{,}{u}_2\right)\mathrm{<mo>\; :\; </mo>}{u}_i\in \mathcal{L}\left(\left[\mathrm{0,}T\right]\right),0\le u_i\le {\bar{u}}_i\mathrm{,}i=\mathrm{1,2}\right\}\mathrm{,}$

where $\mathcal{L}\left(\left[\mathrm{0,}T\right]\right)$ is the space of the Lebesgue integrable functions on $\left[\mathrm{0,}T\right]$, the ${\bar{u}}_i$ are positive constants that may depend on the community (Groups 1 or 2) local resources, $x^u$ denotes the dependence of the state variable $x\in {ℝ}^8$ of the system (16) on the control variables and $x_0\in {ℝ}^8$ represents the initial conditions of the state variable x. The choice of the cost functional is very common in the subject, see [14] [15] [16].

4.1. Existence of an Optimal Control

The following theorem guarantees that the system (16) is controllable this open the possibility to take a strategy to prevent the spread of the disease.

Theorem 8. There exists a pair $\left(u_1^{†}\left(t\right)\mathrm{,}{u}_2^{†}\left(t\right)\right)$ which is an optimal control corresponding to an unique solution

$\left(S_1^{†}\mathrm{,}{E}_1^{†}\mathrm{,}{I}_1^{†}\mathrm{,}{R}_1^{†}\mathrm{,}{S}_2^{†}\mathrm{,}{E}_2^{†}\mathrm{,}{I}_2^{†}\mathrm{,}{R}_2^{†}\right)$

of the system (16), that minimizes the functional $J\left(u_1\mathrm{,}{u}_2\right)$ over the control set U, subject to the constraint state of system (16).

Proof. The optimal control there exists because of the integrand of the functional J, is convex on the closed, convex control set U, a priori boundedness of the state system in relation to the state variables (see [17] ). QED

4.2. The Optimality System

By mean of the Pontryagin’s Maximum Principle [18], necessary and sufficient conditions are derived to be satisfied by the control functions and the corresponding states. Introducing the notation $f_i$ to denote the right hand side of each one of the equations in the system (16) for the i-th state variable of system (16), the Hamiltonian of the Lagrange form (17) is given by:

$\begin{array}{l}\mathcal{H}=d_1{E}_1+d_2{I}_1+X_1{u}_1^2\left(t\right)+{\bar{d}}_1{E}_2+{\bar{d}}_2{I}_2\\ +X_2{u}_2^2\left(t\right)+{\displaystyle \underset{i=1}{\overset{8}{\sum }}}{\Lambda }_i{f}_i.\end{array}$

Therefore,

$\begin{array}{c}\mathcal{H}=d_1{E}_1+d_2{I}_1+X_1{u}_1^2\left(t\right)+{\bar{d}}_1{E}_2+{\bar{d}}_2{I}_2+X_2{u}_2^2\left(t\right)\\ +{\Lambda }_1\left[{\pi }_1+m_{21}{S}_2-\left({\lambda }_1+{\lambda }_2\right)S_1-u_1\left(t\right)S_1-m_{12}{S}_1-{\mu }_1{S}_1\right]\\ +{\Lambda }_2\left[\left({\lambda }_1+{\lambda }_2\right)S_1+u_1\left(t\right)S_1+a_{21}{E}_2-a_{12}{E}_1-{\alpha }_1{E}_1-{\mu }_1{E}_1\right]\\ +{\Lambda }_3\left[{\alpha }_1{E}_1-{\kappa }_1{I}_1-{\mu }_1{I}_1\right]+{\Lambda }_4\left[{\kappa }_1{I}_1+n_{21}{R}_2-n_{12}{R}_1-{\mu }_1{R}_1\right]\end{array}$

$\begin{array}{l}+{\Lambda }_5\left[{\pi }_2+m_{12}{S}_1-\left({\lambda }_1+{\lambda }_2\right)S_2-u_2\left(t\right)S_2-m_{21}{S}_2-{\mu }_2{S}_2\right]\\ +{\Lambda }_6\left[\left({\lambda }_1+{\lambda }_2\right)S_2+u_2\left(t\right)S_2+a_{12}{E}_1-a_{21}{E}_2-{\alpha }_2{E}_2-{\mu }_2{E}_2\right]\\ +{\Lambda }_7\left[{\alpha }_2{E}_2-{\kappa }_2{I}_2-{\mu }_2{I}_2\right]+{\Lambda }_8\left[{\kappa }_2{I}_2+n_{12}{R}_1-n_{21}{R}_2-{\mu }_2{R}_2\right]\mathrm{,}\end{array}$

where the ${\Lambda }_i$, are the adjoint variables that under certain conditions are guaranteed by the Pontryagin’s Maximum Principle [18] as is established in the following theorem.

Theorem 9. For each optimal control $\left(u_1^{†}\left(t\right)\mathrm{,}{u}_2^{†}\left(t\right)\right)$ and its corresponding state solution of the system (16),

$\left(S_1^{†}\mathrm{,}{E}_1^{†}\mathrm{,}{I}_1^{†}\mathrm{,}{R}_1^{†}\mathrm{,}{S}_2^{†}\mathrm{,}{E}_2^{†}\mathrm{,}{I}_2^{†}\mathrm{,}{R}_2^{†}\right)\mathrm{,}$

there exist adjoint variables ${\Lambda }_i,i=1,\cdots ,8$, that satisfy the following adjoint system of differential equations:

$\begin{array}{c}{{\Lambda }^{\prime }}_1={\Lambda }_1\left({\lambda }_1^{†}+{\lambda }_2^{†}-\frac{{\lambda }_1^{†}}{{\underset{\_}{\lambda }}_1^{†}}{S}_1^{†}+u_1^{†}+m_{12}+{\mu }_1\right)\\ -{\Lambda }_2\left({\lambda }_1^{†}+{\lambda }_2^{†}-\frac{{\lambda }_1^{†}}{{\underset{\_}{\lambda }}_1^{†}}{S}_1^{†}+u_1^{†}\right)-{\Lambda }_5\left(\frac{{\lambda }_1^{†}}{{\underset{\_}{\lambda }}_1^{†}}{S}_2^{†}+m_{12}\right)-{\Lambda }_6\frac{{\lambda }_1^{†}}{{\underset{\_}{\lambda }}_1^{†}}{S}_2^{†},\end{array}$

$\begin{array}{c}{{\Lambda }^{\prime }}_2={\Lambda }_1\left(\frac{{\beta }_1{\tau }_e-{\lambda }_1^{†}}{{\underset{\_}{\lambda }}_1^{†}}{S}_1^{†}\right)-{\Lambda }_2\left(\frac{{\beta }_1{\tau }_e-{\lambda }_1^{†}}{{\underset{\_}{\lambda }}_1^{†}}{S}_1^{†}+a_{12}+{\alpha }_1+{\mu }_1\right)\\ -{\alpha }_1{\Lambda }_3+{\Lambda }_5\left(\frac{{\beta }_1{\tau }_e-{\lambda }_1^{†}}{{\underset{\_}{\lambda }}_1^{†}}{S}_2^{†}\right)-{\Lambda }_6\left(\frac{{\beta }_1{\tau }_e-{\lambda }_1^{†}}{{\underset{\_}{\lambda }}_1^{†}}{S}_2^{†}+a_{12}\right)-d_1,\end{array}$

$\begin{array}{c}{{\Lambda }^{\prime }}_3={\Lambda }_1\left(\frac{{\beta }_1-{\lambda }_1^{†}}{{\underset{\_}{\lambda }}_1^{†}}{S}_1^{†}\right)-{\Lambda }_2\left(\frac{{\beta }_1-{\lambda }_1^{†}}{{\underset{\_}{\lambda }}_1^{†}}{S}_1^{†}\right)+{\Lambda }_3\left({\kappa }_1+{\mu }_1\right)\\ -\kappa {\Lambda }_4+{\Lambda }_5\left(\frac{{\beta }_1-{\lambda }_1^{†}}{{\underset{\_}{\lambda }}_1^{†}}{S}_2^{†}\right)-{\Lambda }_6\left(\frac{{\beta }_1-{\lambda }_1^{†}}{{\underset{\_}{\lambda }}_1^{†}}{S}_2^{†}\right)-d_2,\end{array}$

${{\Lambda }^{\prime }}_4=-{\Lambda }_1\frac{{\lambda }_1^{†}}{{\underset{\_}{\lambda }}_1^{†}}{S}_1^{†}+{\Lambda }_2\frac{{\lambda }_1^{†}}{{\underset{\_}{\lambda }}_1^{†}}{S}_1^{†}+{\Lambda }_4\left(n_{12}+{\mu }_1\right)-{\Lambda }_5\frac{{\lambda }_1^{†}}{{\underset{\_}{\lambda }}_1^{†}}{S}_2^{†}+{\Lambda }_6\frac{{\lambda }_1^{†}}{{\underset{\_}{\lambda }}_1^{†}}{S}_2^{†}-n_{12}{\Lambda }_8,$

$\begin{array}{c}{{\Lambda }^{\prime }}_5=-{\Lambda }_1\left(\frac{{\lambda }_2^{†}}{{\underset{\_}{\lambda }}_2^{†}}{S}_2^{†}+m_{21}\right)+{\Lambda }_2\frac{{\lambda }_2^{†}}{{\underset{\_}{\lambda }}_2^{†}}{S}_1^{†}\\ +{\Lambda }_5\left({\lambda }_1^{†}+{\lambda }_2^{†}-\frac{{\lambda }_2^{†}}{{\underset{\_}{\lambda }}_2^{†}}{S}_2^{†}+m_{21}+u_2^{†}+{\mu }_2\right)-{\Lambda }_6\left({\lambda }_1^{†}+{\lambda }_2^{†}-\frac{{\lambda }_2^{†}}{{\underset{\_}{\lambda }}_2^{†}}{S}_2^{†}\right),\end{array}$

$\begin{array}{c}{{\Lambda }^{\prime }}_6={\Lambda }_1\left(\frac{{\beta }_2{\tau }_e-{\lambda }_2^{†}}{{\underset{\_}{\lambda }}_2^{†}}{S}_1^{†}\right)-{\Lambda }_2\left(\frac{{\beta }_2{\tau }_e-{\lambda }_2^{†}}{{\underset{\_}{\lambda }}_2^{†}}{S}_1^{†}+a_{21}\right)-{\Lambda }_5\frac{{\lambda }_2^{†}}{{\underset{\_}{\lambda }}_2^{†}}{S}_2^{†}\\ +{\Lambda }_6\left(\frac{{\lambda }_2^{†}}{{\underset{\_}{\lambda }}_2^{†}}{S}_2^{†}+a_{21}+{\alpha }_2+{\mu }_2\right)-{\alpha }_2{\Lambda }_7-{\bar{d}}_1,\end{array}$

$\begin{array}{c}{{\Lambda }^{\prime }}_7={\Lambda }_1\left(\frac{{\beta }_2-{\lambda }_{2}dag}{S_1^{†}}\right)+{\Lambda }_5\left(\frac{{\beta }_2-{\lambda }_2^{†}}{{\underset{\_}{\lambda }}_2^{†}}{S}_2^{†}\right)-{\Lambda }_6\left(\frac{{\beta }_2-{\lambda }_2^{†}}{{\underset{\_}{\lambda }}_2^{†}}{S}_2^{†}\right)\\ +{\Lambda }_7\left({\kappa }_2+{\mu }_2\right)-{\kappa }_2{\Lambda }_8-{\bar{d}}_2\end{array}$

${{\Lambda }^{\prime }}_8=-{\Lambda }_1\frac{{\lambda }_2^{†}}{{\underset{\_}{\lambda }}_2^{†}}{S}_1^{†}+{\Lambda }_2\frac{{\lambda }_2^{†}}{{\underset{\_}{\lambda }}_2^{†}}{S}_1^{†}-n_{21}{\Lambda }_4-{\Lambda }_5\frac{{\lambda }_2^{†}}{{\underset{\_}{\lambda }}_2^{†}}{S}_2^{†}-{\Lambda }_6{\beta }_2\frac{{\lambda }_2^{†}}{{\underset{\_}{\lambda }}_2^{†}}{S}_2^{†}+{\Lambda }_8\left(n_{21}+{\mu }_2\right),$

with transversality conditions ${\Lambda }_i\left(T\right)=0,i=1,\cdots ,8$. In addition, the optimal controls $u_1^{†}\mathrm{,}{u}_2^{†}$, are given by

$u_1^{†}=\min \left\{{\bar{u}}_1,\max \left\{0,\frac{S_1{\Lambda }_1}{2X_1}\right\}\right\},u_2^{†}=\min \left\{{\bar{u}}_2,\max \left\{0,\frac{S_2{\Lambda }_5}{2X_2}\right\}\right\}.$

Proof. According with the Pontryagin’s Maximum Principle, by mean of the equations:

${{\Lambda }^{\prime }}_1=-\frac{\partial \mathcal{H}}{\partial S_1},\cdots ,{{\Lambda }^{\prime }}_8=-\frac{\partial \mathcal{H}}{\partial R_2},$

the ${\Lambda }_i$ are obtained, with zero final time (transversality) conditions. By other hand, the optimal control set is obtained solving in the interior of the control set U the following equations:

$\frac{\partial \mathcal{H}}{\partial u_1}=0,\frac{\partial \mathcal{H}}{\partial u_2}=0.$

QED

5. Conclusions

A two-group deterministic model is proposed and analyzed for the transmission dynamics of Covid-19 or MERS-CoV (coronaviruses), taking into account fundamentally that the epidemic is in a beginning stage and the groups converge in a mass gathering. The main results are the following:

1) The model has one disease-free equilibrium which is locally asymptotically stable when the reproduction number ${\mathcal{R}}_v<1$. Epidemiologically it implies that a small influx of infected individuals does not produce a large Coronavirus outbreak in mass gatherings, therefore, the disease is subject to control whenever the sizes of the initial sub-populations of the model belong to the attraction basin of the equilibrium DFE (6). When the reproduction number ${\mathcal{R}}_v=1$ the equilibrium DFE (6) lost stability and appears an endemic equilibrium point of positive coordinates that gains stability for ${\mathcal{R}}_v>0$. According to the proof of the Theorem 7, the model exhibits exclusively the phenomenon of forward bifurcation and no backward bifurcation are possible.

2) According to the Pontryagin’s Maximum Principle, the system (5) has a pair of optimal controls for the susceptible populations in both groups. This result suggests that with other actions such as quarantine and vaccination the disease could prevent the disease from spreading.

3) The results obtained in this investigation are purely theoretical. In a research in progress related to the model proposed here, a comparison with experimental results will be established, as well as consideration of delays in the system.

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Malik, T.M., Alsaleh, A.A., Gumel, A.B. and Safi, M.A. (2015) Optimal Strategies for Controlling the MERS Coronavirus during a Mass Gathering. Global Journal of Pure and Applied Mathematics, 11, 4831-4865.
[2]	Bandini, N. (2020) How Atalanta’s Feel-Good Champions League Story Became a “Biological Bomb” for Coronavirus in Italy, Spain. https://www.espn.com/soccer/italian-serie-a/story/4081211/how-atalan-tas-feel-good- champions-league-story-became-a-biological-bomb-for-co-ronavirus-in-italyspain
[3]	Safi, M.A. and Gumel, A.B. (2011) Mathematical Analysis of a Disease Transmission Model with Quarantine, Isolation and an Imperfect Vaccine. Computers & Mathematics with Applications, 61, 3044-3070. https://doi.org/10.1016/j.camwa.2011.03.095
[4]	Safi, M.A. and Gumel, A.B. (2013) Dynamics of a Model with Quarantine-Adjusted Incidence and Quarantine of Susceptible Individuals. Journal of Mathematical Analysis and Applications, 399, 565-575. https://doi.org/10.1016/j.jmaa.2012.10.015
[5]	Safi, M.A. and Gumel, A.B. (2015) Dynamics Analysis of a Quarantine Model in Two Patches. Mathematical Methods in the Applied Sciences, 38, No. 2.
[6]	Thieme, H.R. (2003) Mathematics in Population Biology. Princeton University Press, Princeton. https://doi.org/10.1515/9780691187655
[7]	Lakshmikantham, V., Leela, S. and Martynyuk, A.A. (1989) Stability Analysis of Nonlinear Systems. Marcel Dekker Inc., New York and Basel. https://doi.org/10.1142/1192
[8]	Smith, H.L. and Waltman, P. (1995) The Theory of the Chemostat. Cambridge University Press, Cambridge. https://doi.org/10.1017/CBO9780511530043
[9]	Hethcote, H.W. (2000) The Mathematics Of Infectious Diseases. SIAM Review, 42, 599-653. https://doi.org/10.1137/S0036144500371907
[10]	Van den Driessche, P. and Watmough, J. (2002) Reproduction Number and Sub-Threshold Endemic Equilibria for Compartmental Models of Disease Transmission. Mathematical Biosciences, 180, 29-48. https://doi.org/10.1016/S0025-5564(02)00108-6
[11]	Carr, J. (1981) Applications of Centre Manifold Theory. Springer-Verlag, New York. https://doi.org/10.1007/978-1-4612-5929-9
[12]	Castillo-Chavez, C. and Song, B. (2004) Dynamical Models of Tuberculosis and Their Applications. Mathematical Biosciences and Engineering, 1, 361-404. https://doi.org/10.3934/mbe.2004.1.361
[13]	Ouzi, N. (2010) Deterministic and Stochastic Control, Application to Finance. Lecture Notes, Departement de Mathematiques Appliques, Ecole Polytechnique, Paris.
[14]	Blayneh, K.W., Gumel, A.B., Lenhart, S. and Clayton, T. (2010) Backward Bifurcation and Optimal Control in Transmission Dynamics of West Nile Virus. Bulletin of Mathematical Biology, 72, 1006-1028. https://doi.org/10.1007/s11538-009-9480-0
[15]	Fister, K.R., Lenhart, S. and McNally, J.S. (1998) Optimizing Chemoterapy in an HIV Model. Electronic Journal of Differential Equations, 32, 1-12.
[16]	Jung, E., Lenhart, S. and Feng, Z. (2002) Optimal Control of Treatments in a Two-Strain Tuberculosis Model. Discrete & Continuous Dynamical Systems-B, 2, 473-482. https://doi.org/10.3934/dcdsb.2002.2.473
[17]	Fleming, W.H. and Rishel, R.W. (1975) Deterministic and Stochastic Optimal Control. Springer-Verlag, New York. https://doi.org/10.1007/978-1-4612-6380-7
[18]	Pontryagin, L.S., Boltyanskii, V.G., Gamkrelidze, R.V. and Mishchenko, E.F. (1986) The Mathematical Theory of Optimal Process. Gordon and Breach Science Publishers, Philadelphia.