COVID-19 and dengue co-infection in Brazil: optimal control and cost-effectiveness analysis

Abstract

A mathematical model for the co-interaction of COVID-19 and dengue transmission dynamics is formulated and analyzed. The sub-models are shown to be locally asymptotically stable when the respective reproduction numbers are below unity. Using available data sets, the model is fitted to the cumulative confirmed daily COVID-19 cases and deaths for Brazil (a country with high co-endemicity of both diseases) from February 1, 2021 to September 20, 2021. The fitting was done using the fmincon function in the Optimization Toolbox of MATLAB. Parameters denoting the COVID-19 contact rate, death rate and loss of infection acquired immunity to COVID-19 were estimated using the two data sets. The model is then extended to include optimal control strategies. The appropriate conditions for the existence of optimal control and the optimality system for the co-infection model are established using the Pontryagin’s Principle. Different control strategies and their cost-effectiveness analyses were considered and simulated for the model, which include: controls against incident dengue and COVID-19 infections, control against co-infection with a second disease and treatment controls for both dengue and COVID-19. Highlights of the simulation results show that: (1) dengue prevention strategy could avert as much as 870,000 new COVID-19 infections; (2) dengue only control strategy or COVID-19 only control strategy significantly reduces new co-infection cases; (3) the strategy implementing control against incident dengue infection is the most cost-effective in controlling dengue and COVID-19 co-infections.

1 Introduction

Dengue and Coronavirus disease 2019 (COVID-19) may share clinical and laboratory features [1]. Since 2018, increase in the number of dengue cases in at risk regions to arbovirus outbreaks such as the Reunion Island has been highlighted [1]. The 2019 Coronavirus disease (COVID-19), now a global pandemic, is a respiratory disease caused by the severe acute respiratory syndrome coronavirus 2 (SARSCoV- 2 [2]. In tropical and sub-tropical areas of the world where arboviruses (viral or bacterial infections) and COVID-19 may coexist due to the geographical overlap of the two diseases, clinical diagnosis is difficult, and patients should be tested for both viruses. Verduyn et al.  [1] reported the first confirmed case of co-infection of dengue fever and COVID-19 in a French overseas department located in the Indian Ocean. A comprehensive review of the data on plausible co-infection in a single individuals of dengue and COVID-19 has been reported [3]. Co-epidemics can create a high burden on communities and the health system in the affected areas [4]. Although COVID-19 and dengue are caused by different viruses, the symptomatic appearance of both infections is quite identical and may be hard to distinguish [5]. These similarity of dengue and COVID-19 symptoms could lead to misdiagnosis of one disease for the other and therefore minimizing the extend of co-infection of the two diseases [6].

Due to clinical characteristics and underlying co-morbidities similar to COVID-19 [2, 4, 7], Saddique et al. [5] reported that co-infection of COVID-19 and dengue is an emerging public health concern in dengue endemic countries and investigated what role dengue co-infection plays on the severity and outcome of COVID-19 patients [2]. Co-infection of COVID-19 with vector-borne disease such as malaria is a public health threat [8]. In fact, co-infection of COVID-19 and dengue has already been reported from Asian countries and the Americas, see [9,10,11] and the references therein. These studies highlight high mortality rate in dengue and COVID-19 co-infected patients that may lead to adverse consequences [9]. Dengue viruses circulate throughout the year in Maldives, a dengue holoendemic country [12]. Because clinical and epidemiological criteria may not be sufficient to differentiate COVID-19 and dengue infection, several co-infected patients may be misdiagnosed [13], which could potentially lead to minimizing the extend of the co-infection. Paradoxically, there has been a decrease in dengue cases in Guangzhou (China), mainly attributable to the impact of COVID-19 lockdown in early 2020 [14]. However, the emergence of COVID-19 and dengue co-infection warrants further investigations, at least at the population level to understand the potential of COVID-19 and dengue outbreaks, which could be exacerbated during the post-monsoon months with elevated dengue infections. Also, the dynamics and outcome of a disease may be altered when co-infection with another disease is present [12]. However, non-severe dengue may be more symptomatic than COVID-19 in a co-epidemic dengue endemic settings [7].

Modeling is often the timely option for informing quick decision-making, and to this effect, mathematical modeling for public health purposes has become more refined and used to provide framework for understanding the dynamics of infectious diseases, especially when direct experiments are not possible [8, 15,16,17,18,19,20,21,22,23,24,25]. More recently, Rehman et al. [26] studied a fractional order model for COVID-19, comparing the behavior of the model using different derivatives (Caputo, Caputo–Fabrizio and Atangana–Baleanu) and showed that Caputo presented better results in the form of stability as compared to the other two operators. We formulate and analyze a robust mathematical model for the co-infection of COVID-19 and dengue transmission dynamics, with optimal control and cost-effectiveness analyses. Using available data sets, the proposed model is fitted to the cumulative confirmed daily COVID-19 cases and deaths for Brazil (a country with high co-endemicity of both diseases), and some important parameters are also estimated.

The organization of the rest of the paper is as follows. The proposed co-dynamic model is formulated in Sect. 2 and theoretically analyzed in Sect. 3. By applying Pontryagin’s maximum principle, optimal control of the model to mitigate the spread of both diseases and cost-effectiveness of the interventions is presented in Sect. 5. Numerical simulations performed to support theoretical results and cost-effectiveness analysis are presented in Sect. 6. The conclusion is provided in Sect. 7.

2 The model

Consider a homogeneously mixed population, i.e., individuals in the population have equal probability of contact with each other. Using a deterministic compartmental modeling approach to describe the disease transmission dynamics, at any time t, the total population $N_{\!\!\textsc { h}}$ is subdivided into several epidemiological states depending on individuals health status: susceptible humans $S_{\!\!\textsc { h}}$, infectious individuals with dengue $I_{\!\!\textsc { hd}}$, individuals who have recovered from dengue $R_{\!\!\textsc { hd}}$, infectious individuals with COVID-19 $I_{\!\!\textsc { hc}}$, individuals who have recovered from COVID-19 $R_{\!\!\textsc { hd}}$, infectious individuals with co-infected with dengue and COVID-19 $I_{\!\!\textsc { dc}}$.

The mosquito vector population is given by $N_{\!\!\textsc { v}}$ comprises the susceptible vectors $S_{\!\!\textsc { v}}$, and the infectious vectors with dengue $I_{\!\!\textsc { vd}}$. All the model parameters and their description are provided in 1, while the flows between all the model variables (compartments) are shown in Fig. 1.

The model has the following assumptions:

• i. individuals infected with COVID-19 infection are susceptible to infection with dengue and vice versa.

• ii. co-infected infected individuals can transmit either COVID-19 or dengue but not the mixed infections at the same time,

• iii. co-infected infected individuals can recover either from COVID-19 or dengue but not from the mixed infections at the same time,

• iv. Rate of transmissibility for singly infected and co-infected individuals are assumed same.

Individuals are recruited into the population through birth or immigration at the rate $\Omega _{\!\!\textsc { h}}$. Susceptible humans, $S_{\!\!\textsc { h}}$ acquire COVID-19, following effective contacts with either singly or co-infected individuals with COVID-19 at the rate:

$$\begin{aligned} \lambda _{\!\!\textsc { c}} = \frac{\Lambda _{\!\!\textsc { hc}}(I_{\!\!\textsc { hc}} + I_{\!\!\textsc { dc}})}{N_{\!\!\textsc { h}}}. \end{aligned}$$

(1)

Similarly, the population $S_{\!\!\textsc { h}}$ is reduced due to infection with dengue at the rate:

$$\begin{aligned} \lambda _{\!\!\textsc { d}} = \frac{\Lambda _{\!\!\textsc { vd}}I_{\!\!\textsc { vd}}}{N_{\!\!\textsc { h}}}. \end{aligned}$$

(2)

The parameters $\Lambda _{\!\!\textsc { hc}}$ and $\Lambda _{\!\!\textsc { vd}}$ denote the effective contact rate for the acquisition of COVID-19 and dengue, respectively. The variables in the expressions are defined in Table 1.

Following from the assumptions above, the COVID-19-dengue co-infection model is given by the following system of equations (the flow diagram of the model is presented in Fig. 1, and related parameters of the model are given in Table 1.

Fig. 1

Compartment diagram of the human component of the model

Table 1 Description of the variables and parameters

From Fig. 1, we establish the following system of nonlinear ordinary differential equation describing the dynamics of dengue and COVID-19 co-infection.

$$\begin{aligned} \frac{{\mathrm{d}}S_{\!\!\textsc { h}}}{{\mathrm{d}}t}= & {} \omega _{\!\!\textsc { h}} - \left( \frac{\Lambda _{\!\!\textsc { vd}}I_{\!\!\textsc { vd}}}{N_{\!\!\textsc { h}}} + \frac{\Lambda _{\!\!\textsc { hc}}(I_{\!\!\textsc { hc}} + I_{\!\!\textsc { dc}})}{N_{\!\!\textsc { h}}} \right) S_{\!\!\textsc { h}}- \varrho _{\!\!\textsc { h}} S_{\!\!\textsc { h}} + \eta _{\!\!\textsc { hd}}R_{\!\!\textsc { hd}} + \eta _{\!\!\textsc { hc}}R_{\!\!\textsc { hc}}, \nonumber \\ \frac{{\mathrm{d}}I_{\!\!\textsc { hd}}}{{\mathrm{d}}t}= & {} \frac{\Lambda _{\!\!\textsc { vd}}I_{\!\!\textsc { vd}}}{N_{\!\!\textsc { h}}} (S_{\!\!\textsc { h}} + R_{\!\!\textsc { hc}}) -(\alpha _{\!\!\textsc { hd}} +\varrho _{\!\!\textsc { h}} + \varphi _{\!\!\textsc { hd}}) I_{\!\!\textsc { hd}} - \vartheta _{\!\!\textsc { 1}}\frac{\Lambda _{\!\!\textsc { hc}}(I_{\!\!\textsc { hc}} + I_{\!\!\textsc { dc}})}{N_{\!\!\textsc { h}}}I_{\!\!\textsc { hd}} + \alpha _{\!\!\textsc { hc}}I_{\!\!\textsc { dc}},\nonumber \\ \frac{{\mathrm{d}}R_{\!\!\textsc { hd}}}{{\mathrm{d}}t}= & {} \alpha _{\!\!\textsc { hd}} I_{\!\!\textsc { hd}} - \varrho _{\!\!\textsc { h}} R_{\!\!\textsc { hd}} - \eta _{\!\!\textsc { hd}} R_{\!\!\textsc { hd}} - \frac{\Lambda _{\!\!\textsc { hc}}(I_{\!\!\textsc { hc}} + I_{\!\!\textsc { dc}})}{N_{\!\!\textsc { h}}} R_{\!\!\textsc { hd}},\nonumber \\ \frac{{\mathrm{d}}I_{\!\!\textsc { hc}}}{{\mathrm{d}}t}= & {} \frac{\Lambda _{\!\!\textsc { hc}}(I_{\!\!\textsc { hc}} + I_{\!\!\textsc { dc}})}{N_{\!\!\textsc { h}}} (S_{\!\!\textsc { h}} + R_{\!\!\textsc { hd}}) - (\alpha _{\!\!\textsc { hc}} + \varrho _{\!\!\textsc { h}} + \varphi _{\!\!\textsc { hc}}) I_{\!\!\textsc { hc}} - \vartheta _{\!\!\textsc { 2}}\frac{\Lambda _{\!\!\textsc { vd}}I_{\!\!\textsc { vd}}}{N_{\!\!\textsc { h}}}I_{\!\!\textsc { hc}} + \alpha _{\!\!\textsc { hd}}I_{\!\!\textsc { dc}}, \nonumber \\ \frac{{\mathrm{d}}R_{\!\!\textsc { hc}}}{{\mathrm{d}}t}= & {} \alpha _{\!\!\textsc { hc}} I_{\!\!\textsc { hc}} - \varrho _{\!\!\textsc { h}} R_{\!\!\textsc { hc}} - \eta _{\!\!\textsc { hc}} R_{\!\!\textsc { hc}} - \frac{\Lambda _{\!\!\textsc { vd}}I_{\!\!\textsc { vd}}}{N_{\!\!\textsc { h}}} R_{\!\!\textsc { hc}},\nonumber \\ \frac{{\mathrm{d}}I_{\!\!\textsc { dc}}}{{\mathrm{d}}t}= & {} \vartheta _{\!\!\textsc { 1}}\frac{\Lambda _{\!\!\textsc { hc}}(I_{\!\!\textsc { hc}} + I_{\!\!\textsc { dc}})}{N_{\!\!\textsc { h}}}I_{\!\!\textsc { hd}} + \vartheta _{\!\!\textsc { 2}}\frac{\Lambda _{\!\!\textsc { vd}}I_{\!\!\textsc { vd}}}{N_{\!\!\textsc { h}}}I_{\!\!\textsc { hc}} - (\varrho _{\!\!\textsc { h}} + \varphi _{\!\!\textsc { hd}} + \varphi _{\!\!\textsc { hc}} + \alpha _{\!\!\textsc { hd}} + \alpha _{\!\!\textsc { hc}})I_{\!\!\textsc { dc}}, \nonumber \\ \frac{{\mathrm{d}}S_{\!\!\textsc { vd}}}{{\mathrm{d}}t}= & {} \omega _{\!\!\textsc { d}} - \frac{\Lambda _{\!\!\textsc { hd}}(I_{\!\!\textsc { hd}} + I_{\!\!\textsc { dc}})}{N_{\!\!\textsc { h}}}S_{\!\!\textsc { vd}} - \varrho _{\!\!\textsc { v}}S_{\!\!\textsc { vd}}, \nonumber \\ \frac{{\mathrm{d}}I_{\!\!\textsc { vd}}}{{\mathrm{d}}t}= & {} \frac{\Lambda _{\!\!\textsc { hd}}(I_{\!\!\textsc { hd}} + I_{\!\!\textsc { dc}})}{N_{\!\!\textsc { h}}}S_{\!\!\textsc { vd}} - \varrho _{\!\!\textsc { v}} I_{\!\!\textsc { vd}}, \end{aligned}$$

(3)

with initial conditions

$$\begin{aligned} \begin{array}{ll} S_{\!\!\textsc { h}}(0) \ge 0,\; I_{\!\!\textsc { hd}}(0) \ge 0,\; R_{\!\!\textsc { hd}}(0) \ge 0,\; I_{\!\!\textsc { hc}}(0) \ge 0,\; R_{\!\!\textsc { hc}}(0) \ge 0, I_{\!\!\textsc { dc}}(0) \ge 0, S_{\!\!\textsc { vd}}(0) \ge 0,\; I_{\!\!\textsc { vd}}(0) \ge 0. \end{array}\nonumber \\ \end{aligned}$$

(4)

3 Model analysis

The main focus of our study is on investigating the impact of optimal control on dengue-COCIVD-19 co-dynamics. For this reason, the basic analysis of the dengue-only and COVID-19-only sub-models will focus on deriving the infection threshold parameter that governs the stability of the model equilibria.

3.1 Invariant regions

Since the above model monitors human and mosquito populations, it is assumed that all the state variables and parameters of the model are non-negative for all time $t\ge 0$. The COVID-19 and dengue transmission model (3) will therefore be analyzed in a feasible region $\Omega$.

Lemma 3.1

Solutions of model system (3) are contained in the region $\Omega = \Omega _{\!\!\textsc { h}}\times \Omega _{\!\!\textsc { v}}.$

Proof

Let

$$\begin{aligned} N_{\!\!\textsc { h}}=S_{\!\!\textsc { h}}+ I_{\!\!\textsc { hd}}+ R_{\!\!\textsc { hd}}+ I_{\!\!\textsc { hc}}+ R_{\!\!\textsc { hc}}+ I_{\!\!\textsc { dc}}, \end{aligned}$$

and

$$\begin{aligned} N_{\!\!\textsc { v}}=S_{\!\!\textsc { vd}}+ I_{\!\!\textsc { vd}}. \end{aligned}$$

Assume that $(S_{\!\!\textsc { h}}(t), I_{\!\!\textsc { hd}}(t), R_{\!\!\textsc { hd}}(t), I_{\!\!\textsc { hc}}(t), R_{\!\!\textsc { hc}}(t), I_{\!\!\textsc { dc}}(t)) \in {\mathbb {R}}_+^6$ is a solution of the system with non-negative initial conditions. Then, by summing all the equations of the human-only component of the system (3) we have

$$\begin{aligned} \displaystyle \dot{N_{\!\!\textsc { h}}}=\omega _{\!\!\textsc { h}}-\varrho _{\!\!\textsc { h}} N_{\!\!\textsc { h}}-\varphi _{\!\!\textsc { hd}}I_{\!\!\textsc { hd}} -\varphi _{\!\!\textsc { hc}}I_{\!\!\textsc { hc}}-(\varphi _{\!\!\textsc { hc}}+\varphi _{\!\!\textsc { hd}})I_{\!\!\textsc { dc}} \le \omega _{\!\!\textsc { h}}-\varrho _{\!\!\textsc { h}} N_{\!\!\textsc { h}}, \ \forall t\ge 0. \end{aligned}$$

Thus, on applying Birkhoff and Rota’s Theorem on differential inequality [30], as $t\rightarrow \infty$ we obtain $0\le N_{\!\!\textsc { h}}\le \dfrac{\omega _{\!\!\textsc { h}}}{\varrho _{\!\!\textsc { h}}}$. Therefore, all feasible solutions of the human-only component of the system (3) enters the region

$$\begin{aligned} \Omega _{\!\!\textsc { h}} = \left\{ ( S_{\!\!\textsc { h}}, I_{\!\!\textsc { hd}}, R_{\!\!\textsc { hd}}, I_{\!\!\textsc { hc}}, R_{\!\!\textsc { hc}}, I_{\!\!\textsc { dc}} ) \in {\mathbb {R}}_+^6 : N(t) \leqslant \dfrac{\omega _{\!\!\textsc { h}}}{\varrho _{\!\!\textsc { h}}} \right\} . \end{aligned}$$

Similarly, it can be shown that

$$\begin{aligned} \Omega _{\!\!\textsc { v}} = \left\{ ( S_{\!\!\textsc { vd}}, I_{\!\!\textsc { vd}}) \in {\mathbb {R}}_+^2 : N(t) \leqslant \dfrac{\omega _{\!\!\textsc { d}}}{\varrho _{\!\!\textsc { v}}} \right\} . \end{aligned}$$

Thus, for $t \ge 0$, all possible solutions of (3) will enter the region $\Omega = \Omega _{\!\!\textsc { h}}\times \Omega _{\!\!\textsc { v}},$ which is positively invariant under the flow induced by the model system (3). Also, using the theory of permanence, it can be shown that all solutions on the boundary of $\Omega$ eventually enter the interior of $\Omega$ [31], and the usual existence, uniqueness and continuation results hold. Hence, the model system (3) is well-posed mathematically and epidemiologically, and it is sufficient to consider the dynamics of the flow generated by the model (3) in $\Omega$. Note that the proof of the boundedness of solutions uses the Gronwall’s inequality, see [32]. $\square$

3.2 Analysis of the model without controls

Before analyzing the dynamics of the full model (3), we first analyze the two sub-models namely: COVID-19-only and dengue-only models.

3.3 COVID-19-only model

The COVID-19-only model is obtained by setting $I_{\!\!\textsc { hd}}= R_{\!\!\textsc { hd}} = I_{\!\!\textsc { dc}}=S_{\!\!\textsc { vd}}= I_{\!\!\textsc { vd}}=0$ in (3). Thus, we have,

$$\begin{aligned} \begin{aligned} \frac{{\mathrm{d}}S_{\!\!\textsc { h}}}{{\mathrm{d}}t}&= \omega _{\!\!\textsc { h}} - \frac{\Lambda _{\!\!\textsc { hc}}I_{\!\!\textsc { hc}}}{N_{\!\!\textsc { h}}} S_{\!\!\textsc { h}}- \varrho _{\!\!\textsc { h}} S_{\!\!\textsc { h}} + \eta _{\!\!\textsc { hc}}R_{\!\!\textsc { hc}}, \\ \frac{{\mathrm{d}}I_{\!\!\textsc { hc}}}{{\mathrm{d}}t}&= \frac{\Lambda _{\!\!\textsc { hc}}I_{\!\!\textsc { hc}}}{N_{\!\!\textsc { h}}} S_{\!\!\textsc { h}} - (\alpha _{\!\!\textsc { hc}} + \varrho _{\!\!\textsc { h}} + \varphi _{\!\!\textsc { hc}}) I_{\!\!\textsc { hc}}, \\ \frac{{\mathrm{d}}R_{\!\!\textsc { hc}}}{{\mathrm{d}}t}&= \alpha _{\!\!\textsc { hc}} I_{\!\!\textsc { hc}} - \varrho _{\!\!\textsc { h}} R_{\!\!\textsc { hc}} - \eta _{\!\!\textsc { hc}} R_{\!\!\textsc { hc}}, \\ \end{aligned} \end{aligned}$$

(5)

where, now, the total human population is given by, $N_{\!\!\textsc { h}}=S_{\!\!\textsc { h}}+ I_{\!\!\textsc { hc}}+ R_{\!\!\textsc { hc}}.$ By adding up all the equations of the system (5), we have

$$\begin{aligned} \displaystyle \dot{N_{\!\!\textsc { h}}}=\omega _{\!\!\textsc { h}}-\varrho _{\!\!\textsc { h}} N_{\!\!\textsc { h}}-\varphi _{\!\!\textsc { hc}}I_{\!\!\textsc { hc}} \le \omega _{\!\!\textsc { h}}-\varrho _{\!\!\textsc { h}} N_{\!\!\textsc { h}}. \end{aligned}$$

Consider the region

$$\begin{aligned} \Omega _{\!\!\textsc { hc}} = \left\{ ( S_{\!\!\textsc { h}}, I_{\!\!\textsc { hc}}, R_{\!\!\textsc { hc}} ) : N(t) \leqslant \dfrac{\omega _{\!\!\textsc { h}}}{\varrho _{\!\!\textsc { h}}}\right\} . \end{aligned}$$

Note that the region $\Omega _{\!\!\textsc { hc}}$ is positively invariant [33], and it is sufficient to consider the dynamics of the dengue only sub-model (5) in $\Omega _{\!\!\textsc { hc}}$.

The COVID-19-only model (5) has a DFE given by,

$$\begin{aligned} \epsilon _{\!\!\textsc { hc}}^0=( S_{\!\!\textsc { h}}, I_{\!\!\textsc { hc}}, R_{\!\!\textsc { hc}} )=\left( \dfrac{\omega _{\!\!\textsc { h}}}{\varrho _{\!\!\textsc { h}}},0,0\right) . \end{aligned}$$

One measure of the potential for disease spread in a population is the threshold parameter know as the reproduction number, $R_{0_C}$, which governs the local stability of the DFE of the COVID-19-only model. Using the approach in [34], the associated next generation matrices are given by

$$\begin{aligned}F= \left[ {\begin{array}{cc} \Lambda _{\!\!\textsc { hc}} &{} 0 \\ 0 &{} 0 \\ \end{array}}\right] , \end{aligned}$$

and

$$\begin{aligned}V= \left[ {\begin{array}{cc} \alpha _{\!\!\textsc { hc}} + \varrho _{\!\!\textsc { h}} + \varphi _{\!\!\textsc { hc}} &{} 0 \\ -\alpha _{\!\!\textsc { hc}} &{} \varrho _{\!\!\textsc { h}} + \eta _{\!\!\textsc { hc}} \\ \end{array}}\right] . \end{aligned}$$

The associated basic reproduction number $R_{0_C}$ is given by

$$\begin{aligned} R_{0_C}=\rho (FV^{-1})=\frac{\Lambda _{\!\!\textsc { hc}}}{\alpha _{\!\!\textsc { hc}} + \varrho _{\!\!\textsc { h}} + \varphi _{\!\!\textsc { hc}}}, \end{aligned}$$

(6)

and the following result holds.

Lemma 3.2

The DFE of the COVID-19-only model (5) is locally asymptotically stable if $R_{0_C}<1$ and unstable if $R_{0_C}>1$.

Proof

The stability of the DFE $\epsilon _{\!\!\textsc { hc}}^0=(\dfrac{\omega _{\!\!\textsc { h}}}{\varrho _{\!\!\textsc { h}}},0,0)$ of the COVID-19-only model (5) is obtained from the eigenvalues of the characteristic polynomial, which states that the equilibrium is stable if the eigenvalues of the characteristic polynomial are all negative. For $\epsilon _{\!\!\textsc { hc}}^0$, the Jacobian matrix of the system is obtained as

$$\begin{aligned}J(\epsilon _{\!\!\textsc { hc}}^0)= \left[ {\begin{array}{ccc} - \varrho _{\!\!\textsc { h}} &{} -\Lambda _{\!\!\textsc { hc}}&{} \eta _{\!\!\textsc { hc}}\\ 0&{}\Lambda _{\!\!\textsc { hc}}-(\alpha _{\!\!\textsc { hc}} + \varrho _{\!\!\textsc { h}} + \varphi _{\!\!\textsc { hc}})&{}0 \\ 0&{}\alpha _{\!\!\textsc { hc}} &{} -(\varrho _{\!\!\textsc { h}} + \eta _{\!\!\textsc { hc}}) \\ \end{array}}\right] . \end{aligned}$$

The eigenvalues of characteristic polynomial are given by

$$\begin{aligned} a_1= & {} -\varrho _{\!\!\textsc { h}}, \;\;a_2=-(\varrho _{\!\!\textsc { h}} + \eta _{\!\!\textsc { hc}}), \; \text{ and } \; a_3=\Lambda _{\!\!\textsc { hc}}-(\alpha _{\!\!\textsc { hc}} + \varrho _{\!\!\textsc { h}} + \varphi _{\!\!\textsc { hc}}) \\= & {} \frac{\Lambda _{\!\!\textsc { hc}}}{(\alpha _{\!\!\textsc { hc}} + \varrho _{\!\!\textsc { h}} + \varphi _{\!\!\textsc { hc}})} - 1 = R_{0_C} - 1. \end{aligned}$$

Hence, the DFE $\epsilon _{\!\!\textsc { hc}}^0$ of the COVID-19-only sub-model (5) is locally asymptotically stable if $R_{0_C}<1$. For $R_{0_C} > 1$, prevalence of COVID-19 approaches an endemic equilibrium. $\square$

We will skip the proof of the global stability of the endemic equilibrium (EE) of model system (5) which was carried out in details in [8], where the COVID-19 only model is shown not to undergo the phenomenon of backward bifurcation, consequently, the EE of the model system (5) is globally asymptotically stable.

3.4 Dengue-only model

The dengue-only sub-model is obtained by setting $I_{\!\!\textsc { hc}}= R_{\!\!\textsc { hc}} = I_{\!\!\textsc { dc}}=0$ in (3). Thus, we have

$$\begin{aligned} \begin{aligned} \frac{{\mathrm{d}}S_{\!\!\textsc { h}}}{{\mathrm{d}}t}&= \omega _{\!\!\textsc { h}} - \frac{\Lambda _{\!\!\textsc { vd}}I_{\!\!\textsc { vd}}}{N_{\!\!\textsc { h}}} S_{\!\!\textsc { h}}- \varrho _{\!\!\textsc { h}} S_{\!\!\textsc { h}} + \eta _{\!\!\textsc { hd}}R_{\!\!\textsc { hd}}, \\ \frac{{\mathrm{d}}I_{\!\!\textsc { hd}}}{{\mathrm{d}}t}&= \frac{\Lambda _{\!\!\textsc { vd}}I_{\!\!\textsc { vd}}}{N_{\!\!\textsc { h}}} S_{\!\!\textsc { h}} -(\alpha _{\!\!\textsc { hd}} +\varrho _{\!\!\textsc { h}} + \varphi _{\!\!\textsc { hd}}) I_{\!\!\textsc { hd}}, \\ \frac{{\mathrm{d}}R_{\!\!\textsc { hd}}}{{\mathrm{d}}t}&= \alpha _{\!\!\textsc { hd}} I_{\!\!\textsc { hd}} - \varrho _{\!\!\textsc { h}} R_{\!\!\textsc { hd}} - \eta _{\!\!\textsc { hd}} R_{\!\!\textsc { hd}}, \\ \frac{{\mathrm{d}}S_{\!\!\textsc { vd}}}{{\mathrm{d}}t}&= \omega _{\!\!\textsc { d}} - \frac{\Lambda _{\!\!\textsc { hd}}I_{\!\!\textsc { hd}}}{N_{\!\!\textsc { h}}}S_{\!\!\textsc { vd}} - \varrho _{\!\!\textsc { v}}S_{\!\!\textsc { vd}}, \\ \frac{{\mathrm{d}}I_{\!\!\textsc { vd}}}{{\mathrm{d}}t}&= \frac{\Lambda _{\!\!\textsc { hd}}I_{\!\!\textsc { hd}}}{N_{\!\!\textsc { h}}}S_{\!\!\textsc { vd}} - \varrho _{\!\!\textsc { v}} I_{\!\!\textsc { vd}}, \\ \end{aligned} \end{aligned}$$

(7)

where, now, the total human population is $N_{\!\!\textsc { h}}=S_{\!\!\textsc { h}}+ I_{\!\!\textsc { hd}}+ R_{\!\!\textsc { hd}}$ and the vector population is $N_{\!\!\textsc { v}}=S_{\!\!\textsc { vd}}+ I_{\!\!\textsc { vd}}.$ The feasible region for the sub-model system (7) is

$$\begin{aligned} \Omega _{\!\!\textsc { d}} = \left\{ ( S_{\!\!\textsc { h}}, I_{\!\!\textsc { hd}}, R_{\!\!\textsc { hd}}, S_{\!\!\textsc { vd}}, I_{\!\!\textsc { vd}} ) \in {\mathbb {R}}_+^5 : N_{\!\!\textsc { h}} \leqslant \dfrac{\omega _{\!\!\textsc { h}}}{\varrho _{\!\!\textsc { h}}};N_{\!\!\textsc { v}} \leqslant \dfrac{\omega _{\!\!\textsc { d}}}{\varrho _{\!\!\textsc { v}}} \right\} . \end{aligned}$$

It can be shown that the region $\Omega _{\!\!\textsc { d}}$ is positively invariant (so that it is sufficient to consider the dynamics of the model (7) in $\Omega _{\!\!\textsc { d}}$).

The disease-free equilibrium (DFE) of the dengue-only model (7) is given by

$$\begin{aligned} \epsilon _{\!\!\textsc { d}}^0=\left( S_{\!\!\textsc { h}}, I_{\!\!\textsc { hd}}, R_{\!\!\textsc { hd}},S_{\!\!\textsc { vd}},I_{\!\!\textsc { vd}}\right) =\left( \dfrac{\omega _{\!\!\textsc { h}}}{\varrho _{\!\!\textsc { h}}},0,0,\dfrac{\omega _{\!\!\textsc { v}}}{\varrho _{\!\!\textsc { v}}},0\right) , \end{aligned}$$

and its associated next generation matrices are

$$\begin{aligned}F= \left[ {\begin{array}{ccc} 0&{}0&{}\Lambda _{\!\!\textsc { vd}}\\ 0&{}0&{}0\\ \frac{\Lambda _{\!\!\textsc { hd}}\omega _{\!\!\textsc { v}}\varrho _{\!\!\textsc { h}}}{\omega _{\!\!\textsc { h}}\varrho _{\!\!\textsc { v}}}&{}0 &{} 0 \\ \end{array}}\right] , \end{aligned}$$

and

$$\begin{aligned}V= \left[ {\begin{array}{ccc} \alpha _{\!\!\textsc { hd}} + \varrho _{\!\!\textsc { h}} + \varphi _{\!\!\textsc { hd}} &{} 0&{}0 \\ -\alpha _{\!\!\textsc { hd}} &{} \varrho _{\!\!\textsc { h}} + \eta _{\!\!\textsc { hd}}&{}0 \\ 0&{}0&{}\varrho _{\!\!\textsc { v}}\\ \end{array}}\right] , \end{aligned}$$

so that,

$$\begin{aligned} R_{0_D}=\rho (FV^{-1})=\frac{1}{\varrho _{\!\!\textsc { v}}}\sqrt{\frac{\Lambda _{\!\!\textsc { hd}}\Lambda _{\!\!\textsc { vd}}\varrho _{\!\!\textsc { h}}\omega _{\!\!\textsc { v}}}{\omega _{\!\!\textsc { h}}(\alpha _{\!\!\textsc { hd}} + \varrho _{\!\!\textsc { h}} + \varphi _{\!\!\textsc { hd}})}}. \end{aligned}$$

(8)

The threshold parameter $R_{0_D}$ is the geometric mean of the average number of secondary host infections produced by one vector, and the average number of secondary vector infections produced by one host [35, 36]. In fact, the form of the basic reproduction number in a vector-borne disease is generally a geometric mean between infections caused by hosts and infections caused by vectors [36].

Lemma 3.3

The DFE of the Dengue-only model (7) is locally asymptotically stable if $R_{0_D}<1$ and unstable if $R_{0_D}>1$.

Proof

The stability of the DFE $\epsilon _{\!\!\textsc { d}}^0=( S_{\!\!\textsc { h}}, I_{\!\!\textsc { hd}}, R_{\!\!\textsc { hd}},S_{\!\!\textsc { vd}},I_{\!\!\textsc { vd}})=(\dfrac{\omega _{\!\!\textsc { h}}}{\varrho _{\!\!\textsc { h}}},0,0,\dfrac{\omega _{\!\!\textsc { v}}}{\varrho _{\!\!\textsc { v}}},0),$ of the Dengue-only model (7) is obtained from the eigenvalues of the characteristic polynomial, which states that the equilibrium is stable if the eigenvalues of the characteristic polynomial are all negative. For $\epsilon _{\!\!\textsc { d}}^0$, the Jacobian matrix of the system is obtained as

$$\begin{aligned}J(\epsilon _{\!\!\textsc { d}}^0)= \left[ {\begin{array}{ccccc} - \varrho _{\!\!\textsc { h}} &{} 0 &{} 0 &{} 0 &{} -\Lambda _{\!\!\textsc { vd}}\\ 0&{}-(\alpha _{\!\!\textsc { hc}} + \varrho _{\!\!\textsc { h}} + \varphi _{\!\!\textsc { hc}})&{}0 &{} 0 &{} \Lambda _{\!\!\textsc { vd}} \\ 0&{}\alpha _{\!\!\textsc { hc}} &{} -(\varrho _{\!\!\textsc { h}} + \eta _{\!\!\textsc { hc}}) &{} 0 &{} 0\\ 0&{} -\dfrac{ \Lambda _{\!\!\textsc { hd}} \omega _{\!\!\textsc { v}} \varrho _{\!\!\textsc { h}}}{\omega _{\!\!\textsc { h}}\varrho _{\!\!\textsc { h}}} &{} 0 &{} -\varrho _{\!\!\textsc { v}} &{} 0 \\ 0 &{} \dfrac{ \Lambda _{\!\!\textsc { hd}} \omega _{\!\!\textsc { v}} \varrho _{\!\!\textsc { h}}}{\omega _{\!\!\textsc { h}}\varrho _{\!\!\textsc { h}}} &{} 0 &{} 0 &{} -\varrho _{\!\!\textsc { v}} \end{array}}\right] . \end{aligned}$$

The eigenvalues of characteristic polynomial are given by

$$\begin{aligned} a_1=-\varrho _{\!\!\textsc { h}}, \;\;\; \;a_2=-\varrho _{\!\!\textsc { v}}, \;\;\; \;a_3=-(\varrho _{\!\!\textsc { h}} + \eta _{\!\!\textsc { hc}}), \\ a_4= -\dfrac{1}{2} \sqrt{\frac{\omega _{\!\!\textsc { h}}\varrho _{\!\!\textsc { v}}(\tau -\varrho _{\!\!\textsc { v}})^{2} +4\Lambda _{\!\!\textsc { hc}} \Lambda _{\!\!\textsc { hd}}\varrho _{\!\!\textsc { h}} \omega _{\!\!\textsc { v}} }{\omega _{\!\!\textsc { h}} \varrho _{\!\!\textsc { v}}}} -\dfrac{1}{2}( \tau + \varrho _{\!\!\textsc { v}}), \end{aligned}$$

and

$$\begin{aligned} a_5= \dfrac{1}{2} \sqrt{\frac{\omega _{\!\!\textsc { h}}\varrho _{\!\!\textsc { v}}(\tau -\varrho _{\!\!\textsc { v}})^{2} +4\Lambda _{\!\!\textsc { hc}} \Lambda _{\!\!\textsc { hd}}\varrho _{\!\!\textsc { h}} \omega _{\!\!\textsc { v}} }{\omega _{\!\!\textsc { h}} \varrho _{\!\!\textsc { v}}}} -\dfrac{1}{2}( \tau + \varrho _{\!\!\textsc { v}}), \end{aligned}$$

where $\tau =(\alpha _{\!\!\textsc { hd}} + \varrho _{\!\!\textsc { h}} + \varphi _{\!\!\textsc { hd}}).$ It can be shown that all eigenvalues of the characteristic equation have negative real parts if $R_{0_D}<1$. Hence, the DFE $\epsilon _{\!\!\textsc { hd}}^0$ of the Dengue-only model system (7) is locally asymptotically stable if $R_{0_D}<1$.

$\square$

For $R_{0_D} > 1$, prevalence of dengue approaches an endemic equilibrium.

Other analyses of the dengue-only sub-model have been adequately dealt with in [37], where also, the dengue only sub-model is shown to undergo the phenomenon of backward bifurcation under certain conditions.

4 Dengue-COVID-19 full model

The feasible region for system 3 is given by

$$\begin{aligned} \Omega _{CD} = \Omega _{HC}\times \Omega _{\!\!\textsc { hd}}, \end{aligned}$$

with $\Omega _{HC}$ and $\Omega _{\!\!\textsc { hd}}$ are as defined in the previous sections. It can be shown following the approach in [15, 17] that all solutions of the co-infection Dengue-COVID-19 model system 3 with non-negative initial conditions remain non-negative for all time $t\ge 0$. Also, from the theory of permanence [31], all solutions on the boundary of $\Omega _{CD}$ eventually enter the interior of $\Omega _{CD}$. Thus, $\Omega _{CD}$ is positively invariant and attracting under the flow induced by the system 3

4.1 Stability of the disease-free equilibrium

The disease-free equilibrium of the Dengue-COVID-19 3 is given by

$$\begin{aligned} \displaystyle {\mathcal {E}}_{0} = (S_{\!\!\textsc { h}},I_{\!\!\textsc { hd}},R_{\!\!\textsc { hd}}, I_{\!\!\textsc { hc}},R_{\!\!\textsc { hc}},I_{\!\!\textsc { dc}},S_{\!\!\textsc { vd}},I_{\!\!\textsc { vd}})= \Big (\frac{\Lambda _{\!\!\textsc { hd}}}{\varrho _{\!\!\textsc { h}}},0,0,0,0,0,\frac{\Lambda _{\!\!\textsc { vd}}}{\varrho _{\!\!\textsc { v}}},0\Big ). \end{aligned}$$

(9)

Having derived the basic reproduction numbers for the COVID-19 only and Dengue only sub-models using the next generation method in [34], the associated reproduction number for the full model system 3 is given by

$$\begin{aligned} \displaystyle {\mathcal {R}}_{\!\!\textsc { 0cd}} =\max \left\{ {\mathcal {R}}_{\!\!\textsc { 0hc}}, {\mathcal {R}}_{\!\!\textsc { 0hd}} \right\} . \end{aligned}$$

(10)

The following result follows from Theorem 2 in [34].

Theorem 4.1

The DFE of the Dengue-COVID-19 model 3 is locally asymptotically stable if the threshold parameter ${\mathcal {R}}_{\!\!\textsc { 0cd}} < 1$, and unstable if ${\mathcal {R}}_{\!\!\textsc { 0cd}} > 1$.

5 Optimal control model

In this section, we add time variant controls $u_1(t)$, $u_2(t)$, $u_3(t)$, $u_4(t)$ and $u_5(t)$ into (3) to obtain the optimal interventions for the eradication of COVID-19 and dengue. The controls are defined below:

• i. $u_1$: control against incident dengue infection,

• ii. $u_2$: control against incident COVID-19 infection,

• iii. $u_3$: control against co-infection with a second disease,

• iv. $u_4$: dengue treatment control,

• v. $u_5$: COVID-19 treatment control.

The controls $u_1, u_2$ and $u_3$ satisfy $0 \le u_1 < 0.81$ following from the efficacy of the dengue vaccine reported in [29], $0<u_2<0.90$, following the general efficacy of the COVID-19 vaccine [38], $0< u_3 \le 0.75$ taking the average of both the dengue and COVID-19 vaccine efficacies. The dengue and COVID-19 treatment controls $u_4$ and $u_5$ are bounded as follows: $0 < u_4,u_5 \le 0.80$, with the assumption that treatment cannot be 100% effective against either disease. The control system is presented thus:

$$\begin{aligned} \begin{aligned} \frac{{\mathrm{d}}S_{\!\!\textsc { h}}}{{\mathrm{d}}t}&= \omega _{\!\!\textsc { h}} - \left( (1-u_1)\frac{\Lambda _{\!\!\textsc { vd}}I_{\!\!\textsc { vd}}}{N_{\!\!\textsc { h}}} + (1-u_2)\frac{\Lambda _{\!\!\textsc { hc}}(I_{\!\!\textsc { hc}} + I_{\!\!\textsc { dc}})}{N_{\!\!\textsc { h}}} \right) S_{\!\!\textsc { h}}- \varrho _{\!\!\textsc { h}} S_{\!\!\textsc { h}} + \eta _{\!\!\textsc { hd}}R_{\!\!\textsc { hd}} + \eta _{\!\!\textsc { hc}}R_{\!\!\textsc { hc}}, \\ \frac{{\mathrm{d}}I_{\!\!\textsc { hd}}}{{\mathrm{d}}t}&= \frac{\Lambda _{\!\!\textsc { vd}}I_{\!\!\textsc { vd}}}{N_{\!\!\textsc { h}}} ((1-u_1)S_{\!\!\textsc { h}} + R_{\!\!\textsc { hc}}) -((1+u_4)\alpha _{\!\!\textsc { hd}} +\varrho _{\!\!\textsc { h}} + \varphi _{\!\!\textsc { hd}}) I_{\!\!\textsc { hd}} - (1-u_3)\vartheta _{\!\!\textsc { 1}}\frac{\Lambda _{\!\!\textsc { hc}}(I_{\!\!\textsc { hc}} + I_{\!\!\textsc { dc}})}{N_{\!\!\textsc { h}}}I_{\!\!\textsc { hd}}\\&\quad + (1+u_5)\alpha _{\!\!\textsc { hc}}I_{\!\!\textsc { dc}},\\ \frac{{\mathrm{d}}R_{\!\!\textsc { hd}}}{{\mathrm{d}}t}&= (1+u_4)\alpha _{\!\!\textsc { hd}} I_{\!\!\textsc { hd}} - \varrho _{\!\!\textsc { h}} R_{\!\!\textsc { hd}} - \eta _{\!\!\textsc { hd}} R_{\!\!\textsc { hd}} - \frac{\Lambda _{\!\!\textsc { hc}}(I_{\!\!\textsc { hc}} + I_{\!\!\textsc { dc}})}{N_{\!\!\textsc { h}}} R_{\!\!\textsc { hd}},\\ \frac{{\mathrm{d}}I_{\!\!\textsc { hc}}}{{\mathrm{d}}t}&= \frac{\Lambda _{\!\!\textsc { hc}}(I_{\!\!\textsc { hc}} + I_{\!\!\textsc { dc}})}{N_{\!\!\textsc { h}}} ((1-u_2)S_{\!\!\textsc { h}} + R_{\!\!\textsc { hd}}) - ((1+u_5)\alpha _{\!\!\textsc { hc}} + \varrho _{\!\!\textsc { h}} + \varphi _{\!\!\textsc { hc}}) I_{\!\!\textsc { hc}} - (1-u_3) \vartheta _{\!\!\textsc { 2}}\frac{\Lambda _{\!\!\textsc { vd}}I_{\!\!\textsc { vd}}}{N_{\!\!\textsc { h}}}I_{\!\!\textsc { hc}}\\&\quad + (1+u_4)\alpha _{\!\!\textsc { hd}}I_{\!\!\textsc { dc}}, \\ \frac{{\mathrm{d}}R_{\!\!\textsc { hc}}}{{\mathrm{d}}t}&= (1+u_5)\alpha _{\!\!\textsc { hc}} I_{\!\!\textsc { hc}} - \varrho _{\!\!\textsc { h}} R_{\!\!\textsc { hc}} - \eta _{\!\!\textsc { hc}} R_{\!\!\textsc { hc}} - \frac{\Lambda _{\!\!\textsc { vd}}I_{\!\!\textsc { vd}}}{N_{\!\!\textsc { h}}} R_{\!\!\textsc { hc}},\\ \frac{{\mathrm{d}}I_{\!\!\textsc { dc}}}{{\mathrm{d}}t}&= (1-u_3)\vartheta _{\!\!\textsc { 1}}\frac{\Lambda _{\!\!\textsc { hc}}(I_{\!\!\textsc { hc}} + I_{\!\!\textsc { dc}})}{N_{\!\!\textsc { h}}}I_{\!\!\textsc { hd}} + (1-u_3)\vartheta _{\!\!\textsc { 2}}\frac{\Lambda _{\!\!\textsc { vd}}I_{\!\!\textsc { vd}}}{N_{\!\!\textsc { h}}}I_{\!\!\textsc { hc}} - (\varrho _{\!\!\textsc { h}} + \varphi _{\!\!\textsc { hd}} + \varphi _{\!\!\textsc { hc}} + (1+u_4)\alpha _{\!\!\textsc { hd}} \\&\quad + (1+u_5)\alpha _{\!\!\textsc { hc}})I_{\!\!\textsc { dc}}, \\ \frac{{\mathrm{d}}S_{\!\!\textsc { vd}}}{{\mathrm{d}}t}&= \omega _{\!\!\textsc { d}} - (1-u_1)\frac{\Lambda _{\!\!\textsc { hd}}(I_{\!\!\textsc { hd}} + I_{\!\!\textsc { dc}})}{N_{\!\!\textsc { h}}}S_{\!\!\textsc { vd}} - \varrho _{\!\!\textsc { v}}S_{\!\!\textsc { vd}}, \\ \frac{{\mathrm{d}}I_{\!\!\textsc { vd}}}{{\mathrm{d}}t}&= (1-u_1)\frac{\Lambda _{\!\!\textsc { hd}}(I_{\!\!\textsc { hd}} + I_{\!\!\textsc { dc}})}{N_{\!\!\textsc { h}}}S_{\!\!\textsc { vd}} - \varrho _{\!\!\textsc { v}} I_{\!\!\textsc { vd}}, \\ \end{aligned}\nonumber \\ \end{aligned}$$

(11)

subject to the initial conditions $S_{\!\!\textsc { h}}(0) = S_{\!\!\textsc { h}}^{\!\textsc { 0}}, I_{\!\!\textsc { hd}}(0)=I_{\!\!\textsc { hd}}^{\!\textsc { 0}}, R_{\!\!\textsc { hd}}(0)=R_{\!\!\textsc { hd}}^{\!\textsc { 0}}, I_{\!\!\textsc { hc}}(0)=I_{\!\!\textsc { hc}}^{\!\textsc { 0}}, R_{\!\!\textsc { hc}}(0) = R_{\!\!\textsc { hc}}^{\!\textsc { 0}}, I_{\!\!\textsc { dc}}(0) = I_{\!\!\textsc { dc}}^{\!\textsc { 0}}, S_{\!\!\textsc { vd}}(0) = S_{\!\!\textsc { vd}}^{\!\textsc { 0}}, I_{\!\!\textsc { vd}}(0) =I_{\!\!\textsc { vd}}^{\!\textsc { 0}}$.

The following objective function is considered.

$$\begin{aligned} J\big [u_{\!\!\textsc { 1}}, u_{\!\!\textsc { 2}}, u_{\!\!\textsc { 3}}, u_{\!\!\textsc { 4}}, u_{\!\!\textsc { 5}} \big ]= & {} \int _{0}^{T} \big [c_1 I_{\!\!\textsc { hd}}(t) + c_2 I_{\!\!\textsc { hc}}(t) + c_3 I_{\!\!\textsc { dc}}(t) + c_4 S_{\!\!\textsc { vd}}(t) + c_5 I_{\!\!\textsc { vd}}(t) + \frac{w_1}{2}u_{\!\!\textsc { 1}}^2 \nonumber \\&+\,\frac{w_2}{2}u_{\!\!\textsc { 2}}^2+ \frac{w_3}{2}u_{\!\!\textsc { 3}}^2+ \frac{w_4}{2}u_{\!\!\textsc { 4}}^2+ \frac{w_5}{2}u_{\!\!\textsc { 5}}^2 \big ]dt, \end{aligned}$$

(12)

where T is the final time. The total cost includes the cost of COVID-19 and dengue vaccinations, and other preventive measures and COVID-19 and dengue treatment for all infected individuals. As a result, the nonlinear cost functional is used. In this sequel, we apply the quadratic objective functional for measuring the cost of the control [20]. We seek to find an optimal control, $u_{\!\!\textsc { 1}}^*, u_{\!\!\textsc { 2}}^*, u_{\!\!\textsc { 3}}^*, u_{\!\!\textsc { 4}}^*, u_{\!\!\textsc { 5}}^*$, such that

$$\begin{aligned} J(u_{\!\!\textsc { 1}}^*, u_{\!\!\textsc { 2}}^*, u_{\!\!\textsc { 3}}^*, u_{\!\!\textsc { 4}}^*, u_{\!\!\textsc { 5}}^*)=min \{J(u_{\!\!\textsc { 1}}^*, u_{\!\!\textsc { 2}}^*, u_{\!\!\textsc { 3}}^*, u_{\!\!\textsc { 4}}^*, u_{\!\!\textsc { 5}}^*)|u_{\!\!\textsc { 1}}, u_{\!\!\textsc { 2}}, u_{\!\!\textsc { 3}}, u_{\!\!\textsc { 4}}, u_{\!\!\textsc { 5}}\in U\}, \end{aligned}$$

(13)

where $U= \{(u_{\!\!\textsc { 1}}^*, u_{\!\!\textsc { 2}}^*, u_{\!\!\textsc { 3}}^*, u_{\!\!\textsc { 4}}^*, u_{\!\!\textsc { 5}}^*)\}$, such that $u_{\!\!\textsc { 1}}^*, u_{\!\!\textsc { 2}}^*, u_{\!\!\textsc { 3}}^*$ are measurable with $0\le u_{\!\!\textsc { 1}}^* \le 0.9, 0\le u_{\!\!\textsc { 2}}^* \le 0.9, 0\le u_{\!\!\textsc { 3}}^* \le 0.9$, $0\le u_{\!\!\textsc { 4}}^* \le 1$, $0\le u_{\!\!\textsc { 5}}^* \le 1$ for $t\in [0,T]$ is the control set. The Hamiltonian is given by:

$$\begin{aligned} \begin{aligned} {\mathcal {Z}} \,=\,\,&c_1 I_{\!\!\textsc { hd}}(t) + c_2I_{\!\!\textsc { hc}}(t) + c_3 I_{\!\!\textsc { dc}}(t) +c_4 S_{\!\!\textsc { vd}}(t) + c_5 I_{\!\!\textsc { vd}}(t) + \frac{w_1}{2}u_{\!\!\textsc { 1}}^2+ \frac{w_2}{2}u_{\!\!\textsc { 2}}^2+ \frac{w_3}{2}u_{\!\!\textsc { 3}}^2 + \frac{w_4}{2}u_{\!\!\textsc { 4}}^2+ \frac{w_5}{2}u_{\!\!\textsc { 5}}^2\\&+ \lambda _1 \left( \omega _{\!\!\textsc { h}} - \left( (1-u_1)\frac{\Lambda _{\!\!\textsc { vd}}I_{\!\!\textsc { vd}}}{N_{\!\!\textsc { h}}} + (1-u_2)\frac{\Lambda _{\!\!\textsc { hc}}(I_{\!\!\textsc { hc}} + I_{\!\!\textsc { dc}})}{N_{\!\!\textsc { h}}} \right) S_{\!\!\textsc { h}}- \varrho _{\!\!\textsc { h}} S_{\!\!\textsc { h}} + \eta _{\!\!\textsc { hd}}R_{\!\!\textsc { hd}} + \eta _{\!\!\textsc { hc}}R_{\!\!\textsc { hc}}\right) \\&+ \lambda _2 \left( \frac{\Lambda _{\!\!\textsc { vd}}I_{\!\!\textsc { vd}}}{N_{\!\!\textsc { h}}} ((1-u_1)S_{\!\!\textsc { h}} + R_{\!\!\textsc { hc}}) -((1+u_4)\alpha _{\!\!\textsc { hd}} +\varrho _{\!\!\textsc { h}} + \varphi _{\!\!\textsc { hd}}) I_{\!\!\textsc { hd}} - (1-u_3)\vartheta _{\!\!\textsc { 1}}\frac{\Lambda _{\!\!\textsc { hc}}(I_{\!\!\textsc { hc}} + I_{\!\!\textsc { dc}})}{N_{\!\!\textsc { h}}}I_{\!\!\textsc { hd}} \right. \\&\left. + (1+u_5)\alpha _{\!\!\textsc { hc}}I_{\!\!\textsc { dc}} \right) +\lambda _3 \left( (1+u_4)\alpha _{\!\!\textsc { hd}} I_{\!\!\textsc { hd}} - \varrho _{\!\!\textsc { h}} R_{\!\!\textsc { hd}} - \eta _{\!\!\textsc { hd}} R_{\!\!\textsc { hd}} - \frac{\Lambda _{\!\!\textsc { hc}}(I_{\!\!\textsc { hc}} + I_{\!\!\textsc { dc}})}{N_{\!\!\textsc { h}}} R_{\!\!\textsc { hd}} \right) \\&+\lambda _4 \left( \frac{\Lambda _{\!\!\textsc { hc}}(I_{\!\!\textsc { hc}} + I_{\!\!\textsc { dc}})}{N_{\!\!\textsc { h}}} ((1-u_2)S_{\!\!\textsc { h}} + R_{\!\!\textsc { hd}}) - ((1+u_5)\alpha _{\!\!\textsc { hc}} + \varrho _{\!\!\textsc { h}} + \varphi _{\!\!\textsc { hc}}) I_{\!\!\textsc { hc}} - (1-u_3) \vartheta _{\!\!\textsc { 2}}\frac{\Lambda _{\!\!\textsc { vd}}I_{\!\!\textsc { vd}}}{N_{\!\!\textsc { h}}}I_{\!\!\textsc { hc}}\right. \\&\left. + (1+u_4)\alpha _{\!\!\textsc { hd}}I_{\!\!\textsc { dc}}\right) + \lambda _5 \left( (1+u_5)\alpha _{\!\!\textsc { hc}} I_{\!\!\textsc { hc}} - \varrho _{\!\!\textsc { h}} R_{\!\!\textsc { hc}} - \eta _{\!\!\textsc { hc}} R_{\!\!\textsc { hc}} - \frac{\Lambda _{\!\!\textsc { vd}}I_{\!\!\textsc { vd}}}{N_{\!\!\textsc { h}}} R_{\!\!\textsc { hc}} \right) \\&+ \lambda _6 \left( (1-u_3)\vartheta _{\!\!\textsc { 1}}\frac{\Lambda _{\!\!\textsc { hc}}(I_{\!\!\textsc { hc}} + I_{\!\!\textsc { dc}})}{N_{\!\!\textsc { h}}}I_{\!\!\textsc { hd}} + (1-u_3)\vartheta _{\!\!\textsc { 2}}\frac{\Lambda _{\!\!\textsc { vd}}I_{\!\!\textsc { vd}}}{N_{\!\!\textsc { h}}}I_{\!\!\textsc { hc}} - (\varrho _{\!\!\textsc { h}} + \varphi _{\!\!\textsc { hd}} + \varphi _{\!\!\textsc { hc}} + (1+u_4)\alpha _{\!\!\textsc { hd}} \right. \\&\left. + (1+u_5)\alpha _{\!\!\textsc { hc}})I_{\!\!\textsc { dc}} \right) + \lambda _7 \left( \omega _{\!\!\textsc { d}} - (1-u_1)\frac{\Lambda _{\!\!\textsc { hd}}(I_{\!\!\textsc { hd}} + I_{\!\!\textsc { dc}})}{N_{\!\!\textsc { h}}}S_{\!\!\textsc { vd}} - \varrho _{\!\!\textsc { v}}S_{\!\!\textsc { vd}} \right) \\&+ \lambda _8 \left( (1-u_1)\frac{\Lambda _{\!\!\textsc { hd}}(I_{\!\!\textsc { hd}} + I_{\!\!\textsc { dc}})}{N_{\!\!\textsc { h}}}S_{\!\!\textsc { vd}} - \varrho _{\!\!\textsc { v}} I_{\!\!\textsc { vd}} \right) .\\ \end{aligned}\nonumber \\ \end{aligned}$$

(14)

Theorem 5.1

Suppose the set $\{u_1 , u_2 , u_3, u_4, u_5\}$ minimizes J over U , then we have adjoint variables, $\lambda _1 , \lambda _2, . . . , \lambda _{8}$ (see Appendix for the expressions of $\frac{{\mathrm{d}} \lambda _i}{{\mathrm{d}}t}$) satisfying the adjoint equations

$$\begin{aligned} -\frac{\partial \lambda _i}{\partial t} = \frac{\partial {\mathcal {Z}}}{\partial i}, \end{aligned}$$

with

$$\begin{aligned} \lambda _i(t_f) = 0, \quad \quad \text {where}, \quad i= S_{\!\!\textsc { h}}, I_{\!\!\textsc { hd}}, R_{\!\!\textsc { hd}}, I_{\!\!\textsc { hc}}, R_{\!\!\textsc { hc}}, I_{\!\!\textsc { dc}}, S_{\!\!\textsc { vd}}, I_{\!\!\textsc { vd}}. \end{aligned}$$

(15)

Furthermore,

$$\begin{aligned} \begin{aligned} u_1^*&=\min \left\{ 1,\max \left( 0,\frac{\Lambda _{\!\!\textsc { hd}}(I_{\!\!\textsc { hd}} + I_{\!\!\textsc { dc}})S_{\!\!\textsc { v}} (\lambda _8 - \lambda _7) + I_{\!\!\textsc { vd}}\Lambda _{\!\!\textsc { vd}}S_{\!\!\textsc { h}} (\lambda _2 - \lambda _1)}{w_1N_{\!\!\textsc { h}}}\right) \right\} ,\\ u_2^*&=\min \left\{ 1,\max \left( 0,\frac{\Lambda _{\!\!\textsc { hc}}(I_{\!\!\textsc { hc}} + I_{\!\!\textsc { dc}})S_{\!\!\textsc { h}}(\lambda _4 - \lambda _1)}{w_2N_{\!\!\textsc { h}}}\right) \right\} ,\\ u_3^*&=\min \left\{ 1,\max \left( 0,\frac{I_{\!\!\textsc { hd}}\Lambda _{\!\!\textsc { hc}}(I_{\!\!\textsc { hc}} + I_{\!\!\textsc { dc}}) \vartheta _{\!\!\textsc { 1}} (\lambda _6-\lambda _2) + I_{\!\!\textsc { vd}}\Lambda _{\!\!\textsc { vd}} (I_{\!\!\textsc { hc}} + I_{\!\!\textsc { dc}}) \vartheta _{\!\!\textsc { 2}} (\lambda _6-\lambda _4)}{w_3N_{\!\!\textsc { h}}} \right) \right\} ,\\ u_4^*&=\min \left\{ 1,\max \left( 0,\frac{I_{\!\!\textsc { hd}}\alpha _{\!\!\textsc { hd}}(\lambda _2-\lambda _3) + I_{\!\!\textsc { dc}}\alpha _{\!\!\textsc { hd}}(\lambda _6-\lambda _4)}{w_4}\right) \right\} ,\\ u_5^*&=\min \left\{ 1,\max \left( 0,\frac{I_{\!\!\textsc { hc}}\alpha _{\!\!\textsc { hc}}(\lambda _4-\lambda _5) + I_{\!\!\textsc { dc}}\alpha _{\!\!\textsc { hc}}(\lambda _6-\lambda _2)}{w_5}\right) \right\} ,\\ \end{aligned}\nonumber \\ \end{aligned}$$

(16)

Proof of Theorem 5.1

Consider $U^*=(u_1^*,u_2^*, u_3^*, u_4^*, u_5^*)$ and $S_{\!\!\textsc { h}}^*, I_{\!\!\textsc { hd}}^*, R_{\!\!\textsc { hd}}^*, I_{\!\!\textsc { hc}}^*, R_{\!\!\textsc { hc}}^*, I_{\!\!\textsc { dc}}^*, S_{\!\!\textsc { v}}^*, I_{\!\!\textsc { vd}}^*$ being the associated solutions. Pontryagin’s Maximum Principle [39] is applied, such that there exist adjoint variables satisfying:

$$\begin{aligned} \begin{aligned} - \frac{{\mathrm{d}}\lambda _1}{{\mathrm{d}}t}&=\frac{ \partial {\mathcal {Z}}}{\partial S_{\!\!\textsc { h}}}, \quad \lambda _1(t_f)=0, \quad - \frac{{\mathrm{d}}\lambda _2}{{\mathrm{d}}t}=\frac{ \partial {\mathcal {Z}}}{\partial {I_{\!\!\textsc { hd}}}}, \quad \lambda _2(t_f)=0, \quad - \frac{{\mathrm{d}}\lambda _3}{{\mathrm{d}}t} =\frac{ \partial {\mathcal {Z}}}{\partial {R_{\!\!\textsc { hd}}}}, \quad \lambda _3(t_f)=0,\\ - \frac{{\mathrm{d}}\lambda _4}{{\mathrm{d}}t}&=\frac{ \partial {\mathcal {Z}}}{\partial {I_{\!\!\textsc { hc}}}}, \quad \lambda _4(t_f)=0, - \frac{{\mathrm{d}}\lambda _5}{{\mathrm{d}}t}=\frac{ \partial {\mathcal {Z}}}{\partial {R_{\!\!\textsc { hc}}}}, \quad \lambda _5(t_f)=0, \quad - \frac{{\mathrm{d}}\lambda _6}{{\mathrm{d}}t} =\frac{ \partial {\mathcal {Z}}}{\partial {I_{\!\!\textsc { dc}}}}, \quad \lambda _6(t_f)=0, \\ - \frac{{\mathrm{d}}\lambda _7}{{\mathrm{d}}t}&=\frac{ \partial {\mathcal {Z}}}{\partial {S_{\!\!\textsc { vd}}}}, \quad \lambda _7(t_f)=0, \quad - \frac{{\mathrm{d}}\lambda _8}{{\mathrm{d}}t} =\frac{ \partial {\mathcal {Z}}}{\partial {I_{\!\!\textsc { vd}}}}, \quad \lambda _8(t_f)=0,\\ \end{aligned}\nonumber \\ \end{aligned}$$

(17)

with $\lambda _1(t_f)=\lambda _2(t_f)=\lambda _3(t_f)=\lambda _4(t_f)=\lambda _5(t_f)=\lambda _6(t_f)=\lambda _7(t_f)=\lambda _8(t_f)=0$.

On the interior of the set, where $0<u_j<1$ for all ($j=1,2,3,4,5$), we have that

$$\begin{aligned} \begin{aligned} 0&=\frac{\partial {\mathcal {Z}}}{\partial u_1}=w_1N_{\!\!\textsc { h}}u_1^*-[I_{\!\!\textsc { hd}}\Lambda _{\!\!\textsc { hd}}S_{\!\!\textsc { v}} (\lambda _8 - \lambda _7) + I_{\!\!\textsc { vd}}\Lambda _{\!\!\textsc { vd}}S_{\!\!\textsc { h}} (\lambda _2 - \lambda _1)], \\ 0&=\frac{\partial {\mathcal {Z}}}{\partial u_2}=w_2N_{\!\!\textsc { h}}u_2^*- [I_{\!\!\textsc { hc}}\Lambda _{\!\!\textsc { hc}}S_{\!\!\textsc { h}}(\lambda _4 - \lambda _1)],\\ 0&=\frac{\partial {\mathcal {Z}}}{\partial u_3}=w_3N_{\!\!\textsc { h}}u_3^* - [I_{\!\!\textsc { hc}}I_{\!\!\textsc { hd}}\Lambda _{\!\!\textsc { hc}} \vartheta _{\!\!\textsc { 1}} (\lambda _6-\lambda _2) + I_{\!\!\textsc { hc}}I_{\!\!\textsc { vd}}\Lambda _{\!\!\textsc { vd}} \vartheta _{\!\!\textsc { 2}} (\lambda _6-\lambda _4)],\\ 0&=\frac{\partial {\mathcal {Z}}}{\partial u_3}=w_4u_4^* - [I_{\!\!\textsc { hd}}\alpha _{\!\!\textsc { hd}}(\lambda _2-\lambda _3) + I_{\!\!\textsc { dc}}\alpha _{\!\!\textsc { hd}}(\lambda _6-\lambda _4)],\\ 0&=\frac{\partial {\mathcal {Z}}}{\partial u_3}=w_5u_5^* - [I_{\!\!\textsc { hc}}\alpha _{\!\!\textsc { hc}}(\lambda _4-\lambda _5) + I_{\!\!\textsc { dc}}\alpha _{\!\!\textsc { hc}}(\lambda _6-\lambda _2)].\\ \end{aligned} \end{aligned}$$

(18)

Therefore,

$$\begin{aligned} \begin{aligned} u_1^*&=\frac{\Lambda _{\!\!\textsc { hd}}(I_{\!\!\textsc { hd}} + I_{\!\!\textsc { dc}})S_{\!\!\textsc { v}} (\lambda _8 - \lambda _7) + I_{\!\!\textsc { vd}}\Lambda _{\!\!\textsc { vd}}S_{\!\!\textsc { h}} (\lambda _2 - \lambda _1)}{w_1N_{\!\!\textsc { h}}},\\ u_2^*&=\frac{\Lambda _{\!\!\textsc { hc}}(I_{\!\!\textsc { hc}} + I_{\!\!\textsc { dc}})S_{\!\!\textsc { h}}(\lambda _4 - \lambda _1)}{w_2N_{\!\!\textsc { h}}},\\ u_3^*&=\frac{I_{\!\!\textsc { hd}}\Lambda _{\!\!\textsc { hc}} (I_{\!\!\textsc { hc}} + I_{\!\!\textsc { dc}}) \vartheta _{\!\!\textsc { 1}} (\lambda _6-\lambda _2) + I_{\!\!\textsc { vd}}\Lambda _{\!\!\textsc { vd}} (I_{\!\!\textsc { hc}} + I_{\!\!\textsc { dc}}) \vartheta _{\!\!\textsc { 2}} (\lambda _6-\lambda _4)}{w_3N_{\!\!\textsc { h}}},\\ u_4^*&=\frac{I_{\!\!\textsc { hd}}\alpha _{\!\!\textsc { hd}}(\lambda _2-\lambda _3) + I_{\!\!\textsc { dc}}\alpha _{\!\!\textsc { hd}}(\lambda _6-\lambda _4)}{w_4},\\ u_5^*&=\frac{I_{\!\!\textsc { hc}}\alpha _{\!\!\textsc { hc}}(\lambda _4-\lambda _5) + I_{\!\!\textsc { dc}}\alpha _{\!\!\textsc { hc}}(\lambda _6-\lambda _2)}{w_5}.\\ \end{aligned} \end{aligned}$$

(19)

$$\begin{aligned} u_1^*= & {} \min \left\{ 1,\max \left( 0,\frac{\Lambda _{\!\!\textsc { hd}}(I_{\!\!\textsc { hd}} + I_{\!\!\textsc { dc}})S_{\!\!\textsc { v}} (\lambda _8 - \lambda _7) + I_{\!\!\textsc { vd}}\Lambda _{\!\!\textsc { vd}}S_{\!\!\textsc { h}} (\lambda _2 - \lambda _1)}{w_1N_{\!\!\textsc { h}}}\right) \right\} ,\nonumber \\ u_2^*= & {} \min \left\{ 1,\max \left( 0,\frac{\Lambda _{\!\!\textsc { hc}}(I_{\!\!\textsc { hc}} + I_{\!\!\textsc { dc}})S_{\!\!\textsc { h}}(\lambda _4 - \lambda _1)}{w_2N_{\!\!\textsc { h}}}\right) \right\} ,\nonumber \\ u_3^*= & {} \min \left\{ 1,\max \left( 0,\frac{I_{\!\!\textsc { hd}}\Lambda _{\!\!\textsc { hc}}(I_{\!\!\textsc { hc}} + I_{\!\!\textsc { dc}}) \vartheta _{\!\!\textsc { 1}} (\lambda _6-\lambda _2) + I_{\!\!\textsc { vd}}\Lambda _{\!\!\textsc { vd}} (I_{\!\!\textsc { hc}} + I_{\!\!\textsc { dc}}) \vartheta _{\!\!\textsc { 2}} (\lambda _6-\lambda _4)}{w_3N_{\!\!\textsc { h}}} \right) \right\} ,\nonumber \\ u_4^*= & {} \min \left\{ 1,\max \left( 0,\frac{I_{\!\!\textsc { hd}}\alpha _{\!\!\textsc { hd}}(\lambda _2-\lambda _3) + I_{\!\!\textsc { dc}}\alpha _{\!\!\textsc { hd}}(\lambda _6-\lambda _4)}{w_4}\right) \right\} ,\nonumber \\ u_5^*= & {} \min \left\{ 1,\max \left( 0,\frac{I_{\!\!\textsc { hc}}\alpha _{\!\!\textsc { hc}}(\lambda _4-\lambda _5) + I_{\!\!\textsc { dc}}\alpha _{\!\!\textsc { hc}}(\lambda _6-\lambda _2)}{w_5}\right) \right\} , \end{aligned}$$

(20)

$\square$

5.1 Baseline values of the parameters and model fitting

The total population of Brazil is estimated at 212,559,409 [27], while the life expectancy in Brazil is estimated at 75.88 years. Hence, we set the recruitment rat of humans to $\frac{212,559,409}{75.88 \times 365}$ per day, whereas the natural death rate for humans is set at $\frac{1}{75.88 \times 365}$ per day. The dengue recovery rate is estimated at 0.15 per day [28, 29] so that we set $\alpha _{\!\!\textsc { hd}} = 0.15$ per day. There is no clinical evidence to tell us about the susceptibility of COVID-19 patients to dengue or vice versa. Hence, we set $\vartheta _{\!\!\textsc { 1}} = \vartheta _{\!\!\textsc { 2}}=1$. The other dengue-related parameters are presented in Table 1, together with the references. The initial conditions are set as follows: $S_{\!\!\textsc { h}}(0) = 200,000,000$. The total dengue cases in Brazil for 2021, as at August, 2021 are estimated at 671,732 [40]. Hence, we set $I_{\!\!\textsc { hd}} = 400,000$, $R_{\!\!\textsc { hd}}=4000$, $I_{\!\!\textsc { dc}}=200,000$. The total active COVID-19 cases as at February 1, 2021, were 942,878. Hence, we set $I_{\!\!\textsc { hc}}=942,878$. The total recoveries from COVID-19 as at that same day is 50,925. Thus, we set $R_{\!\!\textsc { hc}}(0)=50,925$.

The model fitting was performed using fmincon function in the Optimization Toolbox of MATLAB [41]. Using the data sets for cumulative confirmed daily COVID-19 cases and deaths for Brazil from February 1, 2021, to September 20, 2021 [27], the COVID-19-related parameters are estimated as follows: $\Lambda _{\!\!\textsc { hc}}$ = 0.1494, $\varphi _{\!\!\textsc { hc}}$ = 0.0047, $\alpha _{\!\!\textsc { hc}}$ = 0.36978. Figures 2 and 3 present the fitting of the model (3) to the cumulative confirmed daily COVID-19 cases and cumulative daily COVID-19 deaths for Brazil from 1 February, 2021 to 20 September, 2021. Both figures show that our model fits well to the two data sets obtained from [27]

Fig. 2

Model fitting to cumulative confirmed daily COVID-19 cases for Brazil from 1 February, 2021 to 20 September, 2021

Fig. 3

Model fitting to cumulative daily COVID-19 deaths for Brazil from 1 February, 2021 to 20 September, 2021

6 Numerical simulations

Simulations carried out on the control system (11), adjoint Eq. (17) and characterizations of the control (20) are run in MATLAB using the forward backward sweep by the Runge–Kutta method. The balancing factors are assumed as follows: $c_1=c_2=c_3=c_4=c_5=1$. Likewise, the quadratic cost functions $\frac{1}{2}w_1u_{\!\!\textsc { 1}}^2, \frac{1}{2}w_2u_{\!\!\textsc { 2}}^2, \frac{1}{2}w_3u_{\!\!\textsc { 3}}^2$, $\frac{1}{2}w_4u_{\!\!\textsc { 4}}^2$ and $\frac{1}{2}w_5u_{\!\!\textsc { 5}}^2$ are applied, over time, in order to compute the total cost for each strategy implemented. The weight constants are set as follows: $w_1=900, w_2=1500, w_3=2000, w_4=1000$ and $w_5=1200$. It is assumed here that, the cost of implementing the dengue prevention control (vaccination and use of treated bed nets and insecticides spray) is less compared to the cost of implementing the COVID-19 prevention control (vaccination, use of hand sanitizers and personal hygiene, use of face-masks in public places and use of personal protective equipments (PPE) by medical personnel). We also assumed that the cost of implementing control against co-infection with a second disease should be more than the cost of preventive control against respective diseases. It is also assumed that the cost of dengue treatment control should be less compared to the cost of COVID-19 treatment control.

6.1 Strategy A: Control against incident dengue infection ($u_1 \ne 0$)

Simulations of the optimal control system (11) when the strategy that prevents incident dengue infection ($u_1 \ne 0$) is administered are presented in Figs. 4, 5, 6 and 7 , respectively. It is observed that when this intervention strategy is implemented, for $\Lambda _{\!\!\textsc { vd}} =5.0, \Lambda _{\!\!\textsc { hd}} = 3.6$ and $\Lambda _{\!\!\textsc { hc}} = 0.5$, so that the reproduction number, $\displaystyle {\mathcal {R}}_{\!\!\textsc { 0cd}} =\max \left\{ {\mathcal {R}}_{\!\!\textsc { 0hc}}, {\mathcal {R}}_{\!\!\textsc { 0hd}} \right\} = 3.4598>1$, there is a significant reduction in the number of individuals infected with dengue (Fig. 4) (as expected). Interestingly, this dengue prevention strategy also averts about 870,000 new COVID-19 cases (as depicted by Fig. 5). Also, worthy of note is that, this strategy against incident dengue infection also has positive population level impact on number of individuals co-infected with dengue and COVID-19 (as shown in Fig. 6, where a total of 4,196,644 new co-infection cases were averted). Moreover, this strategy also averts 346,181 cases of vector infections, thereby significantly reducing the infectious vector population (Fig. 7). The control profile for this control strategy is presented in Fig. 8. The control profile shows that this strategy is at its peak for the first 15 days, drops and then rises again and remains at its peak for the remaining days of the simulation.

Fig. 4

Individuals infected with dengue when strategy A is implemented

Fig. 5

Individuals infected with COVID-19 when strategy A is implemented

Fig. 6

Individuals co-infected with dengue and COVID-19 infections when strategy1 A is implemented

Fig. 7

Infectious vectors infected with dengue when strategy A is implemented

Fig. 8

Control profile for strategy A

6.2 Strategy B: Control against incident COVID-19 infection ($u_2 \ne 0$)

Simulations of the model (11) when the strategy that prevents incident COVID-19 infection ($u_2 \ne 0$) is implemented are presented in Figs. 9, 10 and 11. It is observed that when this strategy is implemented, for $\Lambda _{\!\!\textsc { vd}} =5.0, \Lambda _{\!\!\textsc { hd}} = 3.6$ and $\Lambda _{\!\!\textsc { hc}} = 0.5$, so that the reproduction number, $\displaystyle {\mathcal {R}}_{\!\!\textsc { 0cd}} =\max \left\{ {\mathcal {R}}_{\!\!\textsc { 0hc}}, {\mathcal {R}}_{\!\!\textsc { 0hd}} \right\} = 3.4598>1$, there is a significant reduction in the number of individuals infected with COVID-19 infection (Fig. 9) (as expected). Also, it is imperative to note that, this control against incident COVID-19 infection also has positive impact on number of individuals co-infected with dengue and COVID-19 (Fig. 10, where a total of 4,204,905 new co-infection cases were prevented). In addition, this control also averts about 34,400 vector infections, thereby significantly bringing down the infectious vector population (Fig. 11). The control profile for this control strategy is depicted in Fig. 12. The control profile shows that this strategy is at its peak for the first 15 days, drops and then rises again, and remains at its peak for the remaining days of the simulation.

Fig. 9

Individuals infected with dengue when strategy B is implemented

Fig. 10

Individuals co-infected with dengue and COVID-19 infections when strategy B is implemented

Fig. 11

Infectious vectors infected with dengue when strategy B is implemented

Fig. 12

Control profile for strategy B

6.3 Strategy C: Control against co-infection with a second disease ($u_3 \ne 0$)

The optimal control simulations for the system (11) when the strategy that implements control against co-infection with a second disease ($u_3 \ne 0$) is administered are presented in Figs. 13 and 14 . It is revealed that when this strategy is implemented, for $\Lambda _{\!\!\textsc { vd}} =5.0, \Lambda _{\!\!\textsc { hd}} = 3.6$ and $\Lambda _{\!\!\textsc { hc}} = 0.5$, so that the reproduction number, $\displaystyle {\mathcal {R}}_{\!\!\textsc { 0cd}} =\max \left\{ {\mathcal {R}}_{\!\!\textsc { 0hc}}, {\mathcal {R}}_{\!\!\textsc { 0hd}} \right\} = 3.4598>1$, there is a great reduction in the co-infection new cases (as shown in Fig. 13), where a total of 3,155,000 co-infection new cases were averted. Also, this strategy also averts about 18,300 vector infections, thus significantly bringing down the infectious vector population (Fig. 14). The control profile for this control strategy is given in Fig. 15. The control profile shows that this strategy is at its peak for the entire simulation period.

Fig. 13

Individuals co-infected with dengue and COVID-19 infections when strategy C is implemented

Fig. 14

Infectious vectors infected with dengue when strategy C is implemented

Fig. 15

Control profile for strategy C

6.4 Strategy D: Dengue treatment control ($u_4 \ne 0$)

The control model (11) simulations when dengue treatment strategy alone ($u_4 \ne 0$) is administered are given in Figs. 16, 17 and 18 . It is observed that when dengue treatment strategy alone is implemented, for $\Lambda _{\!\!\textsc { vd}} =5.0, \Lambda _{\!\!\textsc { hd}} = 3.6$ and $\Lambda _{\!\!\textsc { hc}} = 0.5$, so that the reproduction number, $\displaystyle {\mathcal {R}}_{\!\!\textsc { 0cd}} =\max \left\{ {\mathcal {R}}_{\!\!\textsc { 0hc}}, {\mathcal {R}}_{\!\!\textsc { 0hd}} \right\} = 3.4598>1$, there is a gross reduction in the number of individuals infected with dengue (Fig. 16) (as expected). Also, worthy of note is that, this strategy equally has positive population level impact on co-infected individuals (as depicted by Fig. 17, where a total of 2,349,000 new co-infection cases were averted). In addition, this strategy also prevents 55,200 vector infections (as shown in Fig. 18). The control plot is presented in Fig. 19. The control profile shows that this strategy only attains its peak value and is most effective after about 120 days.

Fig. 16

Individuals infected with dengue when strategy D is implemented

Fig. 17

Individuals co-infected with dengue and COVID-19 infections when strategy D is implemented

Fig. 18

Infectious vectors infected with dengue when strategy D is implemented

Fig. 19

Control profile for strategy D

6.5 Strategy E: COVID-19 treatment control ($u_5 \ne 0$)

The control model (11) simulations when COVID-19 treatment control ($u_5 \ne 0$) is the only administered strategy are presented in Figs. 20, 21 and 22 . It is seen from the simulations, which were carried out for $\Lambda _{\!\!\textsc { vd}} =5.0, \Lambda _{\!\!\textsc { hd}} = 3.6$ and $\Lambda _{\!\!\textsc { hc}} = 0.5$, so that the reproduction number, $\displaystyle {\mathcal {R}}_{\!\!\textsc { 0cd}} =\max \left\{ {\mathcal {R}}_{\!\!\textsc { 0hc}}, {\mathcal {R}}_{\!\!\textsc { 0hd}} \right\} = 3.4598>1$, that there is a significant reduction in the number of individuals infected with COVID-19 (Fig. 20) (as expected). Interestingly, this strategy also has positive population level impact on co-infection new cases (as presented in Fig. 21, where a total of 2,544,000 new co-infection cases were averted). Nonetheless, this strategy also averts about 11,200 cases of vector infections, thus bringing down the infectious vector population (Fig. 22). The control plot against time, for this strategy, is given in Fig. 23. The control profile shows that this strategy attains its peak value for the initial days of the simulation, drops and then rises to its peak again after about 140 days.

Fig. 20

Individuals infected with dengue when strategy E is implemented

Fig. 21

Individuals co-infected with dengue and COVID-19 infections when strategy E is implemented

Fig. 22

Infectious vectors infected with dengue when strategy E is implemented

Fig. 23

Control profile for strategy E

6.6 Cost-effectiveness analysis

In this section, we seek to determine the intervention strategy which is most cost-effective in the fight against dengue and COVID-19 co-infections. In order to realize this, the methods used are: the average cost-effectiveness ratio (ACER) and the incremental cost-effectiveness ratio (ICER). ”The cost-effectiveness analysis is used to evaluate the health interventions related benefits so as to justify the costs of the strategies . This is obtained by comparing the differences among the health outcomes and costs of those interventions. ACER deals with a single intervention strategy and weighing the intervention against its baseline option. It is the ratio of the total cost of the intervention to the total number of infection averted by the intervention” [20]. The formula is as follows:

$$\begin{aligned} \text {ACER} = \dfrac{\text {Total cost produced by intervention}}{\text {Total number of infection averted}}. \end{aligned}$$

In a similar manner, ICER deals with the comparison of the differences in the costs and health benefits of two alternate competing interventions. ”It is the ratio of the change in costs of two alternative strategies to the change in the total number of infection averted by the two strategies”. The ICER formula is given by:

$$\begin{aligned} \text {ICER} = \dfrac{\text {Difference in costs between strategies}}{\text {Difference in health effects between strategies}}. \end{aligned}$$

Since our major objective is on how to reduce the co-infection of both dengue and COVID-19 in the population, the total co-infected cases averted and the total cost of the strategies applied shall be used for the cost-effectiveness analysis in this section. These are presented in Table 2. The total cases averted is obtained by calculating the difference in the population of co-infected individuals when control is applied and when control is not applied.

The incremental cost-effectiveness ratio (ICER) for strategies D (dengue treatment control $(u_4 \ne 0)$) and strategy E (COVID-19 treatment control $(u_5 \ne 0)$) are now evaluated.

$$\begin{aligned} \begin{aligned} \text {ICER (D)}&= \frac{1000}{2{,}349{,}000} = 0.0004257,\\ \text {ICER (E)}&= \frac{1200 - 1000}{2{,}544{,}000 - 2{,}349{,}000} = 0.001026.\\ \end{aligned} \end{aligned}$$

Comparing ICER (D) and ICER(E), it is observed that ICER (E) is greater than ICER (D), showing that strategy E is more costly and less effective compared to strategy D. Therefore, strategy E is removed from proceeding ICER computations, shown in Table 3. We shall now compare strategies D and C (Table 4).

Computing ICER for strategies D (dengue treatment control ($u_4 \ne 0$)) and C (control against co-infection with a second disease ($u_3 \ne 0$)), it is observed that ICER (C) is greater than ICER (D), showing that strategy C is more costly and less effective compared to strategy D. Hence, strategy C is removed from subsequent ICER computations. We shall now compare strategies D and A (Table 5).

$$\begin{aligned} \begin{aligned} \text {ICER (D)}&= \frac{1000}{2{,}349{,}000} = 0.0004257,\\ \text {ICER (C)}&= \frac{2000 - 1000}{3{,}155{,}000 - 2{,}349{,}000} = 0.001241.\\ \end{aligned} \end{aligned}$$

Computing ICER for strategies D (dengue treatment control ($u_4 \ne 0$)) and A (control against incident dengue infection ($u_1 \ne 0$)), it is observed that ICER (D) is greater than ICER (A), showing that strategy D is more costly and less effective compared to strategy A. Thus, strategy A is used in subsequent ICER computations, while strategy D is removed. We shall now compare strategies A and B (Table 6).

$$\begin{aligned} \begin{aligned} \text {ICER (D)}&= \frac{1000}{2{,}349{,}000} = 0.0004257,\\ \text {ICER (A)}&= \frac{900 - 1000}{4{,}196{,}644 - 2{,}349{,}000} = -0.00005412.\\ \end{aligned} \end{aligned}$$

Table 2 Increasing order of the total co-infection cases averted using the control strategies

Table 3 ICER computations for strategies D and E

Table 4 ICER computations for strategies D and C

Table 5 ICER computations for strategies D and A

Computing ICER for strategies A (control against incident dengue infection ($u_1 \ne 0$)) and B (control against incident COVID-19 infection ($u_1 \ne 0$)), it is observed that ICER (B) is greater than ICER (A), showing that strategy B is more costly and less effective compared to strategy A. Hence, strategy A is the most cost-effective in controlling dengue and COVID-19 co-infections.

$$\begin{aligned} \begin{aligned} \text {ICER (A)}&= \frac{900}{4{,}196{,}644} = 0.0002145,\\ \text {ICER (B)}&= \frac{1500 - 900}{4{,}204{,}905 - 4{,}196{,}644} = 0.07263.\\ \end{aligned} \end{aligned}$$

Table 6 ICER computations for strategies A and B

7 Conclusion

We have formulated and analyzed a mathematical model for the co-infection and COVID-19 and dengue transmission dynamics, with optimal control and cost-effectiveness analysis. The sub-models are shown to be locally asymptotically stable when the respective reproduction numbers are below unity. Using available data sets, the model is fitted to the cumulative confirmed daily COVID-19 cases and deaths for Brazil (a country with high co-endemicity of both diseases) from February 1, 2021 to September 20, 2021. The fitting was done using the fmincon function in the Optimization Toolbox of MATLAB. Parameters denoting the COVID-19 contact rate, death rate and loss of infection acquired immunity to COVID-19 were estimated using the two data sets. The appropriate conditions for the existence of optimal control and the optimality system for the co-infection model are established using the Pontryagin’s Principle. Different control strategies were considered and simulated for the model, which include: controls against incident dengue and COVID-19 infections, control against co-infection with a second disease and treatment controls for both dengue and COVID-19. Highlights of the simulation results show that:

• i. dengue prevention strategy could avert as much as 870,000 new COVID-19 infections (as depicted in Fig. 5);

• ii. dengue only control strategy or COVID-19 only control strategy significantly reduces new co-infection cases. These are presented in Figs. 6 and 10 ;

• iii. the strategy that implements control against incident COVID-19 infection (vaccination, use of hand sanitizers and personal hygiene, use of face-masks in public places and use of personal protective equipments (PPE) by medical personnel) averts the highest number of co-infection cases than any of the administered control strategies. This is depicted by Fig. 10, where a total of 4,204,905 new co-infection cases were averted.

• iv. the strategy implementing control against incident dengue infection is the most cost-effective in controlling dengue and COVID-19 co-infections as presented in Sect. 6.6).

As COVID-19 pandemic threatens the delivery of dengue services, the lockdown could have further impacted the continuation of dengue control programs, and there is an urgent need for rapid and effective responses to avoid dengue outbreaks.

Our model was created based on the focus of COVID-19 and dengue co-infection only. Thus, we did not investigate the impact of multiple COVID-19 stains and waves on the dynamics of the two diseases. Also, the emergence of COVID-19 and dengue co-infection warrants further investigations at the individual level (within host dynamics). While this is to the best of our knowledge seemingly the first study of the co-interaction of COVID-19 and dengue, more studies should be devoted to the mathematical (agent based, within/intra-host dynamics) and epidemiological dynamics of this co-infection. Due to the uncertainty of several aspects and characteristics of both diseases, we realize the difficulty in finding estimates for certain parameters in this study. Specific parameters difficult to estimate are the modification parameters accounting for susceptibility of dengue infected individuals to COVID-19 and the modification parameter accounting for susceptibility of COVID-19 infected individuals to dengue. Because mathematical models are symbolic representations of biological systems, by construction, they inherit the loss of information which could potentially make the prediction of model outcomes imprecise. Therefore, exploring sensitivity analysis of the model variables and parameters to changes in the assumptions made regarding the characteristics of the disease is viable.

It is important to note from Fig. 15 that the control $u_3$ for the co-infection is optimally applied for the entire simulation period. On the contrary, the other four controls $u_1, u_2, u_4$ and $u_5$ all start at a high lever, but deepens after a short time before returning to the optimum. This striking result may be due to the initial inadequate treatment when the disease emerged, and vaccine hesitancy/denial at the onset of the vaccination campaign. It has been reported how hesitancy/denial is challenging the vaccination campaigns with the number of those infected increasing relative to the percentage of population that avoids getting vaccinated [42]. This is an area that warrants further investigation as well as the impact of multiple waves of COVID-19 on disease co-interaction.

Appendix: Adjoint functions of optimality system

$$\begin{aligned} \lambda _1^{'}&= {\lambda _2}\, \left( \frac{I_{\!\!\textsc { v}}\, {\beta }_{\!\!\textsc { v}}\, \left( {u_1} - 1\right) }{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) } + \frac{I_{\!\!\textsc { v}}\, {\beta }_{\!\!\textsc { v}}\, \left( R_{\!\!\textsc { c}} - S_{\!\!\textsc { h}}\, \left( {u_1} - 1\right) \right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2} \right. \\&\quad \left. + \frac{I_{\!\!\textsc { d}}\, {\beta }_{\!\!\textsc { c}}\, {\vartheta }_{\!\!\textsc { 1}}\, \left( {u_3} - 1\right) \, \left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { dc}}\right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2}\right) \\&\quad - {\lambda _6}\, \left( \frac{I_{\!\!\textsc { c}}\, I_{\!\!\textsc { v}}\, {\beta }_{\!\!\textsc { v}}\, {\vartheta }_{\!\!\textsc { 2}}\, \left( {u_3} - 1\right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2} + \frac{I_{\!\!\textsc { d}}\, {\beta }_{\!\!\textsc { c}}\, {\vartheta }_{\!\!\textsc { 1}}\, \left( {u_3} - 1\right) \, \left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { dc}}\right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2}\right) \\&\quad + {\lambda _1}\, \left( {\varrho }_{\!\!\textsc { h}} + S_{\!\!\textsc { h}}\, \left( \frac{{\beta }_{\!\!\textsc { c}}\, \left( {u_2} - 1\right) \, \left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { dc}}\right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2} + \frac{I_{\!\!\textsc { v}}\, {\beta }_{\!\!\textsc { v}}\, \left( {u_1} - 1\right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2}\right) \right. \\&\quad \left. - \frac{{\beta }_{\!\!\textsc { c}}\, \left( {u_2} - 1\right) \, \left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { dc}}\right) }{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) } - \frac{I_{\!\!\textsc { v}}\, {\beta }_{\!\!\textsc { v}}\, \left( {u_1} - 1\right) }{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }\right) \\&\quad + {\lambda _4}\, \left( \frac{{\beta }_{\!\!\textsc { c}}\, \left( {u_2} - 1\right) \, \left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { dc}}\right) }{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) } \right. \\&\quad \left. + \frac{{\beta }_{\!\!\textsc { c}}\, \left( R_{\!\!\textsc { d}} - S_{\!\!\textsc { h}}\, \left( {u_2} - 1\right) \right) \, \left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { dc}}\right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2} + \frac{I_{\!\!\textsc { c}}\, I_{\!\!\textsc { v}}\, {\beta }_{\!\!\textsc { v}}\, {\vartheta }_{\!\!\textsc { 2}}\, \left( {u_3} - 1\right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2}\right) \\&\quad - \frac{R_{\!\!\textsc { d}}\, {\beta }_{\!\!\textsc { c}}\, {\lambda _3}\, \left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { dc}}\right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2} \\&\quad - \frac{I_{\!\!\textsc { v}}\, R_{\!\!\textsc { c}}\, {\beta }_{\!\!\textsc { v}}\, {\lambda _5}}{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2} + \frac{S_{\!\!\textsc { v}}\, {\beta }_{\!\!\textsc { d}}\, {\lambda _7}\, \left( {u_1} - 1\right) \, \left( I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}}\right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2}\\&\quad - \frac{S_{\!\!\textsc { v}}\, {\beta }_{\!\!\textsc { d}}\, {\lambda _8}\, \left( {u_1} - 1\right) \, \left( I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}}\right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} \quad + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2},\\ \lambda _2^{'}&= {\lambda _8}\, \left( \frac{S_{\!\!\textsc { v}}\, {\beta }_{\!\!\textsc { d}}\, \left( {u_1} - 1\right) }{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) } - \frac{S_{\!\!\textsc { v}}\, {\beta }_{\!\!\textsc { d}}\, \left( {u_1} - 1\right) \, \left( I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}}\right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2}\right) \\&\quad - {\lambda _7}\, \left( \frac{S_{\!\!\textsc { v}}\, {\beta }_{\!\!\textsc { d}}\, \left( {u_1} - 1\right) }{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }\right. \\&\quad \left. - \frac{S_{\!\!\textsc { v}}\, {\beta }_{\!\!\textsc { d}}\, \left( {u_1} - 1\right) \, \left( I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}}\right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2}\right) - {c_1} \\&\quad - {\lambda _3}\, \left( {\alpha }_{\!\!\textsc { d}}\, \left( {u_4} + 1\right) + \frac{R_{\!\!\textsc { d}}\, {\beta }_{\!\!\textsc { c}}\, \left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { dc}}\right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2}\right) \\ \end{aligned}$$

$$\begin{aligned}&\quad - {\lambda _6}\, \left( \frac{I_{\!\!\textsc { c}}\, I_{\!\!\textsc { v}}\, {\beta }_{\!\!\textsc { v}}\, {\vartheta }_{\!\!\textsc { 2}}\, \left( {u_3} - 1\right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2} - \frac{{\beta }_{\!\!\textsc { c}}\, {\vartheta }_{\!\!\textsc { 1}}\, \left( {u_3} - 1\right) \, \left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { dc}}\right) }{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) } \right. \\&\quad \left. + \frac{I_{\!\!\textsc { d}}\, {\beta }_{\!\!\textsc { c}}\, {\vartheta }_{\!\!\textsc { 1}}\, \left( {u_3} - 1\right) \, \left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { dc}}\right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2}\right) \\&\quad + {\lambda _4}\, \left( \frac{{\beta }_{\!\!\textsc { c}}\, \left( R_{\!\!\textsc { d}} - S_{\!\!\textsc { h}}\, \left( {u_2} - 1\right) \right) \, \left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { dc}}\right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2} + \frac{I_{\!\!\textsc { c}}\, I_{\!\!\textsc { v}}\, {\beta }_{\!\!\textsc { v}}\, {\vartheta }_{\!\!\textsc { 2}}\, \left( {u_3} - 1\right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2}\right) + {\lambda _2}\, \left( {\varphi }_{\!\!\textsc { d}} + {\varrho }_{\!\!\textsc { h}} + {\alpha }_{\!\!\textsc { d}}\, \left( {u_4} + 1\right) \right. \\&\quad \left. + \frac{I_{\!\!\textsc { v}}\, {\beta }_{\!\!\textsc { v}}\, \left( R_{\!\!\textsc { c}} - S_{\!\!\textsc { h}}\, \left( {u_1} - 1\right) \right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2} - \frac{{\beta }_{\!\!\textsc { c}}\, {\vartheta }_{\!\!\textsc { 1}}\, \left( {u_3} - 1\right) \, \left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { dc}}\right) }{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) } + \frac{I_{\!\!\textsc { d}}\, {\beta }_{\!\!\textsc { c}}\, {\vartheta }_{\!\!\textsc { 1}}\, \left( {u_3} - 1\right) \, \left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { dc}}\right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2}\right) \\&\quad + S_{\!\!\textsc { h}}\, {\lambda _1}\, \left( \frac{{\beta }_{\!\!\textsc { c}}\, \left( {u_2} - 1\right) \, \left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { dc}}\right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2} + \frac{I_{\!\!\textsc { v}}\, {\beta }_{\!\!\textsc { v}}\, \left( {u_1} - 1\right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2}\right) \\&\quad - \frac{I_{\!\!\textsc { v}}\, R_{\!\!\textsc { c}}\, {\beta }_{\!\!\textsc { v}}\, {\lambda _5}}{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2},\\ \lambda _3^{'}&= {\lambda _4}\, \left( \frac{{\beta }_{\!\!\textsc { c}}\, \left( R_{\!\!\textsc { d}} - S_{\!\!\textsc { h}}\, \left( {u_2} - 1\right) \right) \, \left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { dc}}\right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2} - \frac{{\beta }_{\!\!\textsc { c}}\, \left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { dc}}\right) }{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) } + \frac{I_{\!\!\textsc { c}}\, I_{\!\!\textsc { v}}\, {\beta }_{\!\!\textsc { v}}\, {\vartheta }_{\!\!\textsc { 2}}\, \left( {u_3} - 1\right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2}\right) \\&\quad - {\lambda _1}\, \left( {\eta }_{\!\!\textsc { d}} - S_{\!\!\textsc { h}}\, \left( \frac{{\beta }_{\!\!\textsc { c}}\, \left( {u_2} - 1\right) \, \left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { dc}}\right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2} + \frac{I_{\!\!\textsc { v}}\, {\beta }_{\!\!\textsc { v}}\, \left( {u_1} - 1\right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2}\right) \right) \\&\quad - {\lambda _6}\, \left( \frac{I_{\!\!\textsc { c}}\, I_{\!\!\textsc { v}}\, {\beta }_{\!\!\textsc { v}}\, {\vartheta }_{\!\!\textsc { 2}}\, \left( {u_3} - 1\right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2} + \frac{I_{\!\!\textsc { d}}\, {\beta }_{\!\!\textsc { c}}\, {\vartheta }_{\!\!\textsc { 1}}\, \left( {u_3} - 1\right) \, \left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { dc}}\right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2}\right) \\&\quad + {\lambda _2}\, \left( \frac{I_{\!\!\textsc { v}}\, {\beta }_{\!\!\textsc { v}}\, \left( R_{\!\!\textsc { c}} - S_{\!\!\textsc { h}}\, \left( {u_1} - 1\right) \right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2} + \frac{I_{\!\!\textsc { d}}\, {\beta }_{\!\!\textsc { c}}\, {\vartheta }_{\!\!\textsc { 1}}\, \left( {u_3} - 1\right) \, \left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { dc}}\right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2}\right) \\&\quad + {\lambda _3}\, \left( {\eta }_{\!\!\textsc { d}} + {\varrho }_{\!\!\textsc { h}} + \frac{{\beta }_{\!\!\textsc { c}}\, \left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { dc}}\right) }{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) } - \frac{R_{\!\!\textsc { d}}\, {\beta }_{\!\!\textsc { c}}\, \left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { dc}}\right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2}\right) \\&\quad - \frac{I_{\!\!\textsc { v}}\, R_{\!\!\textsc { c}}\, {\beta }_{\!\!\textsc { v}}\, {\lambda _5}}{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2} + \frac{S_{\!\!\textsc { v}}\, {\beta }_{\!\!\textsc { d}}\, {\lambda _7}\, \left( {u_1} - 1\right) \, \left( I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}}\right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2} - \frac{S_{\!\!\textsc { v}}\, {\beta }_{\!\!\textsc { d}}\, {\lambda _8}\, \left( {u_1} - 1\right) \, \left( I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}}\right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2},\\ \end{aligned}$$

$$\begin{aligned} \begin{aligned} \lambda _4^{'}&= {\lambda _6}\, \left( \frac{I_{\!\!\textsc { d}}\, {\beta }_{\!\!\textsc { c}}\, {\vartheta }_{\!\!\textsc { 1}}\, \left( {u_3} - 1\right) }{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) } + \frac{I_{\!\!\textsc { v}}\, {\beta }_{\!\!\textsc { v}}\, {\vartheta }_{\!\!\textsc { 2}}\, \left( {u_3} - 1\right) }{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) } - \frac{I_{\!\!\textsc { c}}\, I_{\!\!\textsc { v}}\, {\beta }_{\!\!\textsc { v}}\, {\vartheta }_{\!\!\textsc { 2}}\, \left( {u_3} - 1\right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2} \right. \\&\left. \quad - \frac{I_{\!\!\textsc { d}}\, {\beta }_{\!\!\textsc { c}}\, {\vartheta }_{\!\!\textsc { 1}}\, \left( {u_3} - 1\right) \, \left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { dc}}\right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2}\right) - {c_2} + {\lambda _3}\, \left( \frac{R_{\!\!\textsc { d}}\, {\beta }_{\!\!\textsc { c}}}{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) } \right. \\&\left. \qquad \qquad - \frac{R_{\!\!\textsc { d}}\, {\beta }_{\!\!\textsc { c}}\, \left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { dc}}\right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2}\right) \\&\quad - {\lambda _5}\, \left( {\alpha }_{\!\!\textsc { c}}\, \left( {u_5} + 1\right) + \frac{I_{\!\!\textsc { v}}\, R_{\!\!\textsc { c}}\, {\beta }_{\!\!\textsc { v}}}{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2}\right) + {\lambda _2}\, \left( \frac{I_{\!\!\textsc { v}}\, {\beta }_{\!\!\textsc { v}}\, \left( R_{\!\!\textsc { c}} - S_{\!\!\textsc { h}}\, \left( {u_1} - 1\right) \right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2} \right. \\&\quad \left. - \frac{I_{\!\!\textsc { d}}\, {\beta }_{\!\!\textsc { c}}\, {\vartheta }_{\!\!\textsc { 1}}\, \left( {u_3} - 1\right) }{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) } + \frac{I_{\!\!\textsc { d}}\, {\beta }_{\!\!\textsc { c}}\, {\vartheta }_{\!\!\textsc { 1}}\, \left( {u_3} - 1\right) \, \left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { dc}}\right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2}\right) + {\lambda _4}\, \left( {\varphi }_{\!\!\textsc { c}} + {\varrho }_{\!\!\textsc { h}} + {\alpha }_{\!\!\textsc { c}}\, \left( {u_5} + 1\right) \right. \\&\quad \left. - \frac{{\beta }_{\!\!\textsc { c}}\, \left( R_{\!\!\textsc { d}} - S_{\!\!\textsc { h}}\, \left( {u_2} - 1\right) \right) }{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) } + \frac{{\beta }_{\!\!\textsc { c}}\, \left( R_{\!\!\textsc { d}} - S_{\!\!\textsc { h}}\, \left( {u_2} - 1\right) \right) \, \left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { dc}}\right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2} - \frac{I_{\!\!\textsc { v}}\, {\beta }_{\!\!\textsc { v}}\, {\vartheta }_{\!\!\textsc { 2}}\, \left( {u_3} - 1\right) }{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) } \right. \\&\quad \left. + \frac{I_{\!\!\textsc { c}}\, I_{\!\!\textsc { v}}\, {\beta }_{\!\!\textsc { v}}\, {\vartheta }_{\!\!\textsc { 2}}\, \left( {u_3} - 1\right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2}\right) + S_{\!\!\textsc { h}}\, {\lambda _1}\, \left( \frac{{\beta }_{\!\!\textsc { c}}\, \left( {u_2} - 1\right) \, \left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { dc}}\right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2}\right. \\&\left. \quad - \frac{{\beta }_{\!\!\textsc { c}}\, \left( {u_2} - 1\right) }{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) } \right. \\&\quad \left. + \frac{I_{\!\!\textsc { v}}\, {\beta }_{\!\!\textsc { v}}\, \left( {u_1} - 1\right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2}\right) + \frac{S_{\!\!\textsc { v}}\, {\beta }_{\!\!\textsc { d}}\, {\lambda _7}\, \left( {u_1} - 1\right) \, \left( I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}}\right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2} - \frac{S_{\!\!\textsc { v}}\, {\beta }_{\!\!\textsc { d}}\, {\lambda _8}\, \left( {u_1} - 1\right) \, \left( I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}}\right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2},\\ \end{aligned} \end{aligned}$$

$$\begin{aligned} \lambda _5^{'}&= {\lambda _2}\, \left( \frac{I_{\!\!\textsc { v}}\, {\beta }_{\!\!\textsc { v}}\, \left( R_{\!\!\textsc { c}} - S_{\!\!\textsc { h}}\, \left( {u_1} - 1\right) \right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2} - \frac{I_{\!\!\textsc { v}}\, {\beta }_{\!\!\textsc { v}}}{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }\right. \\&\quad \left. + \frac{I_{\!\!\textsc { d}}\, {\beta }_{\!\!\textsc { c}}\, {\vartheta }_{\!\!\textsc { 1}}\, \left( {u_3} - 1\right) \, \left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { dc}}\right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2}\right) \\&\quad - {\lambda _6}\, \left( \frac{I_{\!\!\textsc { c}}\, I_{\!\!\textsc { v}}\, {\beta }_{\!\!\textsc { v}}\, {\vartheta }_{\!\!\textsc { 2}}\, \left( {u_3} - 1\right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2} + \frac{I_{\!\!\textsc { d}}\, {\beta }_{\!\!\textsc { c}}\, {\vartheta }_{\!\!\textsc { 1}}\, \left( {u_3} - 1\right) \, \left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { dc}}\right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2}\right) \\&\quad - {\lambda _1}\, \left( {\eta }_{\!\!\textsc { c}} - S_{\!\!\textsc { h}}\, \left( \frac{{\beta }_{\!\!\textsc { c}}\, \left( {u_2} - 1\right) \, \left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { dc}}\right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2} + \frac{I_{\!\!\textsc { v}}\, {\beta }_{\!\!\textsc { v}}\, \left( {u_1} - 1\right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2}\right) \right) \\&\quad + {\lambda _5}\, \left( {\eta }_{\!\!\textsc { c}} + {\varrho }_{\!\!\textsc { h}} + \frac{I_{\!\!\textsc { v}}\, {\beta }_{\!\!\textsc { v}}}{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) } - \frac{I_{\!\!\textsc { v}}\, R_{\!\!\textsc { c}}\, {\beta }_{\!\!\textsc { v}}}{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2}\right) \\&\quad + {\lambda _4}\, \left( \frac{{\beta }_{\!\!\textsc { c}}\, \left( R_{\!\!\textsc { d}} - S_{\!\!\textsc { h}}\, \left( {u_2} - 1\right) \right) \, \left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { dc}}\right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2} + \frac{I_{\!\!\textsc { c}}\, I_{\!\!\textsc { v}}\, {\beta }_{\!\!\textsc { v}}\, {\vartheta }_{\!\!\textsc { 2}}\, \left( {u_3} - 1\right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2}\right) \\&\quad - \frac{R_{\!\!\textsc { d}}\, {\beta }_{\!\!\textsc { c}}\, {\lambda _3}\, \left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { dc}}\right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2} \\&\quad + \frac{S_{\!\!\textsc { v}}\, {\beta }_{\!\!\textsc { d}}\, {\lambda _7}\, \left( {u_1} - 1\right) \, \left( I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}}\right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2} - \frac{S_{\!\!\textsc { v}}\, {\beta }_{\!\!\textsc { d}}\, {\lambda _8}\, \left( {u_1} - 1\right) \, \left( I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}}\right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2},\\ \end{aligned}$$

$$\begin{aligned} \lambda _6^{'}&= {\lambda _8}\, \left( \frac{S_{\!\!\textsc { v}}\, {\beta }_{\!\!\textsc { d}}\, \left( {u_1} - 1\right) }{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} +\quad I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) } - \frac{S_{\!\!\textsc { v}}\, {\beta }_{\!\!\textsc { d}}\, \left( {u_1} - 1\right) \, \left( I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}}\right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2}\right) \\&\quad - {\lambda _4}\, \left( {\alpha }_{\!\!\textsc { d}}\, \left( {u_4} + 1\right) + \frac{{\beta }_{\!\!\textsc { c}}\, \left( R_{\!\!\textsc { d}} - S_{\!\!\textsc { h}}\, \left( {u_2} - 1\right) \right) }{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) } - \frac{{\beta }_{\!\!\textsc { c}}\, \left( R_{\!\!\textsc { d}} - S_{\!\!\textsc { h}}\, \left( {u_2} - 1\right) \right) \, \left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { dc}}\right) }{{\left( I_{\!\!\textsc { c}} +\quad I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2} \right. \\&\left. \quad - \frac{I_{\!\!\textsc { c}}\, I_{\!\!\textsc { v}}\, {\beta }_{\!\!\textsc { v}}\, {\vartheta }_{\!\!\textsc { 2}}\, \left( {u_3} - 1\right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} +\quad I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2}\right) \\&\quad - {\lambda _7}\, \left( \frac{S_{\!\!\textsc { v}}\, {\beta }_{\!\!\textsc { d}}\, \left( {u_1} - 1\right) }{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) } - \frac{S_{\!\!\textsc { v}}\, {\beta }_{\!\!\textsc { d}}\, \left( {u_1} - 1\right) \, \left( I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}}\right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2}\right) - {c_3} \\&\quad - {\lambda _2}\, \left( {\alpha }_{\!\!\textsc { c}}\, \left( {u_5} + 1\right) - \frac{I_{\!\!\textsc { v}}\, {\beta }_{\!\!\textsc { v}}\, \left( R_{\!\!\textsc { c}} - S_{\!\!\textsc { h}}\, \left( {u_1} - 1\right) \right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2} + \frac{I_{\!\!\textsc { d}}\, {\beta }_{\!\!\textsc { c}}\, {\vartheta }_{\!\!\textsc { 1}}\, \left( {u_3} - 1\right) }{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }\right. \\&\quad \left. - \frac{I_{\!\!\textsc { d}}\, {\beta }_{\!\!\textsc { c}}\, {\vartheta }_{\!\!\textsc { 1}}\, \left( {u_3} - 1\right) \, \left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { dc}}\right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2}\right) \\&\quad + {\lambda _3}\, \left( \frac{R_{\!\!\textsc { d}}\, {\beta }_{\!\!\textsc { c}}}{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) } - \frac{R_{\!\!\textsc { d}}\, {\beta }_{\!\!\textsc { c}}\, \left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { dc}}\right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2}\right) \\&\quad + {\lambda _6}\, \left( {\varphi }_{\!\!\textsc { c}} + {\varphi }_{\!\!\textsc { d}} + {\varrho }_{\!\!\textsc { h}} + {\alpha }_{\!\!\textsc { c}}\, \left( {u_5} + 1\right) + {\alpha }_{\!\!\textsc { d}}\, \left( {u_4} + 1\right) \right. \\&\quad \left. + \frac{I_{\!\!\textsc { d}}\, {\beta }_{\!\!\textsc { c}}\, {\vartheta }_{\!\!\textsc { 1}}\, \left( {u_3} - 1\right) }{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) } - \frac{I_{\!\!\textsc { c}}\, I_{\!\!\textsc { v}}\, {\beta }_{\!\!\textsc { v}}\, {\vartheta }_{\!\!\textsc { 2}}\, \left( {u_3} - 1\right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2}\right. \\&\quad \left. - \frac{I_{\!\!\textsc { d}}\, {\beta }_{\!\!\textsc { c}}\, {\vartheta }_{\!\!\textsc { 1}}\, \left( {u_3} - 1\right) \, \left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { dc}}\right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2}\right) \\&\quad + S_{\!\!\textsc { h}}\, {\lambda _1}\, \left( \frac{{\beta }_{\!\!\textsc { c}}\, \left( {u_2} - 1\right) \, \left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { dc}}\right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2} - \frac{{\beta }_{\!\!\textsc { c}}\, \left( {u_2} - 1\right) }{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) } \right. \\&\quad \left. + \frac{I_{\!\!\textsc { v}}\, {\beta }_{\!\!\textsc { v}}\, \left( {u_1} - 1\right) }{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2}\right) - \frac{I_{\!\!\textsc { v}}\, R_{\!\!\textsc { c}}\, {\beta }_{\!\!\textsc { v}}\, {\lambda _5}}{{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }^2},\\ \lambda _7^{'}&= {\lambda _7}\, \left( {\varrho }_{\!\!\textsc { v}} - \frac{{\beta }_{\!\!\textsc { d}}\, \left( {u_1} - 1\right) \, \left( I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}}\right) }{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }\right) - {c_4} + \frac{{\beta }_{\!\!\textsc { d}}\, {\lambda _8}\, \left( {u_1} - 1\right) \, \left( I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}}\right) }{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }\\ \end{aligned}$$

$$\begin{aligned} \lambda _8^{'}&= {\lambda _8}\, {\varrho }_{\!\!\textsc { v}} - {c_5} + \frac{R_{\!\!\textsc { c}}\, {\beta }_{\!\!\textsc { v}}\, {\lambda _5}}{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) } - \frac{{\beta }_{\!\!\textsc { v}}\, {\lambda _2}\, \left( R_{\!\!\textsc { c}} - S_{\!\!\textsc { h}}\, \left( {u_1} - 1\right) \right) }{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) } \\&\quad - \frac{S_{\!\!\textsc { h}}\, {\beta }_{\!\!\textsc { v}}\, {\lambda _1}\, \left( {u_1} - 1\right) }{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) } - \frac{I_{\!\!\textsc { c}}\, {\beta }_{v}\, {\lambda _4}\, {\vartheta }_{2}\, \left( {u_3} - 1\right) }{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) } + \frac{I_{\!\!\textsc { c}}\, {\beta }_{\!\!\textsc { v}}\, {\lambda _6}\, {\vartheta }_{\!\!\textsc { 2}}\, \left( {u_3} - 1\right) }{\left( I_{\!\!\textsc { c}} + I_{\!\!\textsc { d}} + I_{\!\!\textsc { dc}} + R_{\!\!\textsc { c}} + R_{\!\!\textsc { d}} + S_{\!\!\textsc { h}}\right) }, \end{aligned}$$