COVID-19 and Chikungunya: an optimal control model with consideration of social and environmental factors.

Chikungunya is one of the Aedes aegypti diseases that mosquito transmits to humans and that are common in tropical countries like Yemen. In this work, we formulated a novel dynamic mathematical model framework, which integrates COVID-19 and Chikungunya outbreaks. The proposed model is governed by a system of dynamic ordinary differential equations (ODEs). Particle swarm optimization was employed to solve the parameters estimation problem of the outbreaks of COVID-19 and Chikungunya in Yemen (March 1, 2020, to May 30, 2020). Besides, a bi-objective optimal control model was formulated, which minimizes the number of affected individuals and minimizes the total cost associated with the intervention strategies. The bi-objective optimal control was also solved using PSO. Five preventive measures were considered to curb the environmental and social factors that trigger the emergence of these viruses. Several strategies were simulated to evaluate the best possible strategy under the conditions and available resources in Yemen. The results obtained confirm that the strategy, which provides resources to prevent the transmission of Chikungunya and provides sufficient resources for testing, applying average social distancing, and quarrying the affected individuals, has a significant effect on flattening the epidemic curves and is the most suitable strategy in Yemen.

Introduction

Yemen is a country located in the southwest of Asia on 2000 km coastline, and it is considered one of the tropical countries. Therefore, infectious and seasonal diseases are highly prevalent there, such as Aedes aegypti, Aedes albopictus, dengue, Chikungunya, malaria, cholera, and other infectious diseases. Recently, the Chikungunya virus (CHIKV) has spread in some areas, especially in Al-Hadidah, Lahj and Aden. Spread of the Chikungunya virus coincided with the COVID-19 pandemic, causing many deaths. It is well known in epidemiology that Chikungunya is an infectious disease caused by a virus of the Togaviridae family and transmitted to humans by either the genus Aedes aegypti or Aedes albopictus mosquitoes. COVID-19 is an infectious disease caused by the virus SARS-CoV-2 of the Coronaviridae family and transmitted to humans by infectious droplets of another human. Hence, most of the studies that addressed epidemiological optimal control models focused on break of the vector-borne using some intervention strategies.

Lately, several studies have been proposed to explore optimal control models for infectious diseases including Chikungunya diseases. Ruiz-Moreno et al. (2012) combined stochastic model, which is related to season-based mosquito vector dynamics with an epidemiological model, to investigate the potential risk of Chikungunya entering the US. Moulay et al. (2012) formulated an optimal control problem for the Chikungunya model, which is controlled by three time-dynamic variables: prevention and treatment and vector breeding sites destruction. Yakob and Clements (2013) addressed the Chikungunya outbreak in Réunion 2006, and Monte Carlo simulation is used in the sensitivity analysis study. Liu and Stechlinski (2015) discussed the spread of CHIKV outbreak in Réunion with respect to some time-varying parameters that are related to the rate of breeding of mosquitoes in the rainy and dry climate and the rate of dynamic contact between the mosquitoes and the human. Agusto et al. (2016) developed a transmission model with focusing on the influence of age on Chikungunya transmission. The population is classified into three age-structured: Juveniles, adults and senior subcategories. Zhu et al. (2018) employed the input–output polynomials with QR decomposition to estimate the parameters of a CHIKV transmission model. Liu et al. (2020) investigated the impact of temperature and rainfall on the CHIKV model. Dodero-Rojas et al. (2020) formulated a compartmental mathematical model, which was implemented in the city of Rio de Janeiro in the period 2017–2019 on the Chikungunya epidemic. The basic reproduction number was estimated, and a potential outbreak of Mayaro virus was predicted. Moreover, a number of different intervention strategies were simulated with the aim of reducing the number of infected individuals. The optimal control model was proposed for the Chikungunya epidemic model in Gonzalez-Parra et al. (2020) when Pontryagin’s maximum principle was utilized to solve the optimal controls and to find the optimal final time. This study was also carried out on real data in Colombia 2015. The optimal control policies of the Chikungunya epidemic model were presented through three time-varying variables to control the spread of the CHIKV between humans and vectors. The first time-varying variable related to educational awareness and personal protection through bed nets, wearing full body clothing, and avoiding water stagnation. The second time-varying variable relied on the impact of the treatment of infected individuals. The third time-varying variable related to spraying insecticides in order to reduce mosquito population to reduce the proportion of the infected individuals. (Ali et al. 2020) presented a clinical study on the Chikungunya outbreak in Yemen (March 2020–May 2020).

On the other hand, this year, coinciding with the outbreak of the COVID-19, many studies have emerged that addressed the epidemiological model of COVID-19. (Abbasi et al. 2020) proposed a dynamic of SQEIAR Susceptible, Quarantined, Exposed, Infected, Asymptomatic, and Recovered individuals model considering two time-varying variables to control the outbreak: quarantine and treatment. Besides, they used Pontryagin’s maximum principle to solve the optimal controls model. Kouidere et al. (2020) considered three actions to control the COVID-19 outbreak, i.e. prevention, quarantine, and treatment. Also, Pontryagin’s maximum principle utilized and analyzed the model. The Susceptible-Exposed-Asymptomatic-Infected-Removed (SEAIR) model was expanded to include the perished class due to infection with COVID-19 in Tsay et al. (2020). Social distancing, quarantining, and availability of testing kits were used as time-dynamic variables for dynamic optimal control. The aim of the optimization problem was to minimize social and economic costs provided that the size of the epidemic remains under a given epidemic peak value. The optimal control model was solved directly by Python by considering the discretization of time domain into one finite element per day. Morato et al. (2020a, b) proposed epidemic models of COVID-19 with respect to the real data from Brazil. Social distancing was investigated as the time-varying parameter for controlling the optimization model. Ullah and Altaf (2020) formulated a mathematical model to explore COVID-19 in Pakistan, presenting a mathematical analytic study. The parameters were estimated using a least square fitting method. Moreover, two actions were investigated for controlling the epidemic, namely, quarantine, and hospitalization. Also, some articles addressed the optimal distribution of vaccines, such as Narayanamoorthy et al. (2021), Hezam et al. (2021b) or discussed the humanitarian response plan for high-priority countries as an optimal distribution model (Hezam 2021a).

There are a number of epidemiological models in the literature that try to describe the co-infection models. Samat and Ma’Arof (2014) presented a stochastic susceptible–infective–recovered (human); susceptible–infective (vector) SIR-SI model to estimate the relative risk of Chikungunya and dengue infection in Malaysia. They determined the high-risk and low-risk area of Dengue and Chikungunya, which can be employed for the prevention and control of both outbreaks. Kumar et al. (2019) studied the Chikungunya, dengue and zika outbreaks in Mexico 2015–2016. They estimated the parameters for each virus. Besides, they presented some strategies based on isolation and self-protection factors to control the outbreaks. Aldila and Agustin (2018) proposed a transmission model of dengue-Chikungunya coinfection and presented the mathematical analysis and numerical simulation. In Isea and Lonngren (2016), transmission dynamics models were developed of Dengue, Chikungunya and Zika, including the possibility of a dual-infection of any two infectious diseases in the same population. Musa et al. (2020) proposed a co-infection transmission dynamics model for Chikungunya and dengue. The model was implemented on real data from India and a study of sensitivity analysis for model parameters was presented. An interesting study was presented in Jindal (2020), in which the authors investigated the effect of the lockdowns due to COVID-19 on outbreaks of mosquito-borne diseases, including dengue, Chikungunya and malaria. They found that the risk and severity of the outbreak increased with lockdowns. Hezam et al. (2021a) addressed the COVID-19 and Cholera coinfection in Yemen 2020. Four time-dynamic variables were investigated to show its impact on epidemics curbing: social distancing, quarantine, test kits availability, rate of the population able to access to pure water. Also, some studies discussed the interactive COVID-19 with other outbreaks like dengue (Lam et al. 2020), HIV (Doungmo Goufo et al. 2020), Ebola (Zhang and Jain 2020), or the impact of COVID-19 pandemic on TP patients (Marimuthu et al. 2020), or on unemployment problem (Hezam 2021b).

Limited studies have been carried out with respect to solving the epidemiological models using metaheuristic algorithms. Yan and Zou (2008) used genetic algorithm to find sub-optimal solution of the optimal control in SARS epidemics. Florentino et al. (2014) employed genetic algorithm to solve multi-objective optimal control models derived from a dynamical model of dengue virus epidemic. Florentino et al. (2018) also used genetic algorithm to solve the optimal control model that aims to reduce Aedes mosquitoes through preventing breeding or spreading insecticides. Chaikham and Sawangtong (2018) used differential evolution to solve the optimal control model applied to curb the Zika virus outbreak. Rahmalia and Herlambang (2018) applied artificial bee colony algorithm to find the optimal weights of objective function in optimal control derived from Influenza epidemic models. In Akman et al. (2018), particle swarm optimization was employed to solve the parameter estimation problem in ordinary differential equations. Putra and Mu (2019) used particle swarm optimization to estimate the parameters, then Euler method was used to solve Susceptible, Infected and Resistant (SIR) models. Kmet and Kmetova (2019) introduced two approaches for solving the optimal control model, i.e. the direct method and the indirect method. The direct methods were based on collaborating Bernstein–Bézier parametrization for control variables, in addition to invasive weed optimization algorithm and particle swarm optimization for solving the optimal control model, while Pontryagin’s maximum principle and the necessary conditions of optimal control model are the indirect methods. Windarto et al. (2020) employed particle swarm optimization to estimate the parameters of dengue transmission model. Lobato et al. (2020) applied stochastic fractal search algorithm to solve both parameters estimation problem and multi-objective optimal control model of COVID-19 outbreak in China. Mahmoodabadi (2020) used PSO to find the unknown variables in the epidemic model, comparing it to other classical methods. Yousefpour et al. (2020) used a genetic algorithm to solve the multi-objective optimal control model of COVID-19 outbreak. Besides, some control strategies were discussed and analyzed, which can contribute to curbing the COVID-19 epidemic. Similarly, in Okuonghae and Omame (2020), social distancing, use of face mask and case detection of COVID-19 outbreak in Nigeria were the investigated actions. Salgotra et al. (2020) used Genetic Evolutionary Programming (GEP) to predict confirmed cases and death cases of COVID-19 in 15 worst affected countries in the world for the period (February 1st, 2020, to last of May 2020. He et al. (2020) employed PSO to estimate the parameters of the SEIR model based on real data from Hubei province. A similar study, Abdallah and Nafea (2021), used PSO to estimate the parameters of the SEIQRD model and predict the COVID-19 outbreak in Italy.

However, up to the time of writing this paper, no study has formulated a mathematical model of the coinfection COVID-19 and Chikungunya. The contributions of this work are as follows. Firstly, the present study formulates a new dynamic transmission model of the coinfection COVID-19 and Chikungunya. Secondly, PSO will be used to estimate the parameters based on historical data for both infections in Yemen from March 1, 2020, to May 30, 2020, and also to predict the trajectories of both outbreaks for one year. Thirdly, we present an optimal control model that includes two objective functions: minimizing the expected cumulative number of infected people for COVID-19 and Chikungunya, and minimizing the total cost associated with the intervention strategies. The optimal solutions of the optimal control model will be found also by PSO. Fourthly, we investigate the sensitivity of the time-dynamic inputs of the optimal control by evaluating a novel set of control strategies. The focus is on time-dynamic inputs related to lockdown, social distancing, test kits numbers, personal protection and insecticides and larvicides to fight mosquitoes. These five actions are used to identify the optimal strategies that contribute to infection mitigation. Lastly, all the strategies are carried out for COVID-19 and Chikungunya outbreaks in Yemen to determine the suitable optimal strategy in Yemen under the conditions and available resources.

The remainder of this work is organized as follows. In Sect. 2, we discuss the co-infection model formulation. PSO is presented in Sect. 3, while the estimated parameters are discussed in Sect. 4. In Sect. 5, we provide the optimal control model and its analysis. In Sect. 6, we discuss several strategies. Finally, in Sect. 7, we summarize and conclude the work.

Model formulation

To describe the co-infection dynamics of COVID-19 and Chikungunya, the dynamical compartment model includes both human and mosquito population. The human population will be divided into sixteen different epidemiological classes, while mosquito vector will include three classes. Descriptions of the present classes with associate parameters are defined in Tables 1 and 2, respectively.

Table 1

Description and initial values of model variables

Variables	Description	Initial conditions	Source
N_H	Total population size of humans	29, 825, 964	(Worldometers 2020)
S_H	Number of susceptible humans to both COVID-19 and CHIKV	N_H	Assumed
E_1	Number of humans exposed to COVID-19	0,N_H×10^-6	(Tsay et al. 2020)
A_1	Number of asymptomatic COVID-19 individuals	0	Assumed
I_1	Number of COVID-19- infected individuals	0	Assumed
R_1	Number of recovered from COVID-19 individuals	0	Assumed
P_1	Number of perished with COVID-19 individuals	0	Assumed
I_R2	Number of infected to CHIKV after recovery from COVID-19	0	Assumed
E_2	Number of humans exposed to CHIKV	0	Assumed
A_2	Number of asymptomatic CHIKV individuals	0	Assumed
I_2	Number of CHIKV infected individuals	0	Assumed
R_2	Number of recovered from CHIKV individuals	0	Assumed
P_2	Number of perished with CHIKV individuals	0	Assumed
I_1R	Number of infected to COVID-19 after recovery from CHIKV	0	Assumed
I_12	Number of dually infected with both COVID-19 and CHIKV	0	Assumed
R_12	Number of recovered from both COVID-19 and CHIKV	0	Assumed
P_12	Number of perished by both COVID-19 and CHIKV	0	Assumed
N_m	The total population of mosquitoes	10^3,5×10^4	(Sanchez et al. 2018)
X	Number of susceptible mosquitoes to CHIKV	10^3,5×10^4	(Sanchez et al. 2018)
Y	Number of mosquitoes exposed to CHIKV	0	(Sanchez et al. 2018)
Z	Number of CHIKV-infected mosquitoes	1	(Sanchez et al. 2018)

Table 2

Parameters description and their values

Parameters	Description	Value (range)	Reference
I1peak	The peak limit of the infected individuals with COVID-19	10^4,28×10^6	Assumed
I2peak	The peak limit of the infected individuals with CHIKV	10^4,28×10^6	Assumed
Z^peak	The peak limit of the infected mosquitoes	10^4,10^5	Assumed
u_1t	Time-dynamic function to measure the social distancing rate	0.05,0.5	(Tsay et al. 2020)
u_2t	Time-dynamic function to measure the quarantining rate	0.01,0.3	(Tsay et al. 2020)
u_3t	Time-dynamic function to measure the testing rate	0.1,0.3	(Tsay et al. 2020)
u_4t	Time-dynamic function to measure the personal protection rate	0.2,1	Assumed
u_5t	Time-dynamic function to measure the ratio of used the insecticides and larvicides for reducing the mosquitoes population	0.2,1	Assumed
α_c	Rate of exposed to CHIKV	0.0209	Estimated
α_1	Transmission rate of COVID-19	0.09	Estimated
α_2	Transmission rate of CHIKV	0.0780	Estimated
α_12	Rate of dually infected with both COVID-19 and CHIKV	0	Assumed
α_3	Rate of mosquitoes exposed to CHIKV	0.3574	Estimated
φ_1	Rate of infected to COVID-19 Asymptomatically-Symptomatically	0.55	Estimated
φ_2	Rate of infected to CHIKV Asymptomatically-Symptomatically	0.5	Estimated
δ_1	Rate of infected to CHIKV of COVID-19 patient	0	Assumed
δ_2	Rate of infected to COVID-19 of CHIKV patient	0	Assumed
β_1,β_2	Recovery rate from COVID-19 and CHIKV respectively	0.0964, 0.98	Estimated
ω_1,ω_2,ω_12	Death rate from COVID-19, CHIKV and both respectively	0.0974, 0.0127, 0.0889	Estimated
μ_1,μ_2,μ_12	Susceptible for COVID-19, CHIKV and both again	0.0889,0.09,0.07	Estimated
λ_1	Infected rate by CHIKV after recovered from COVID-19	0	Assumed
λ_2	Infected rate by COVID-19 after recovered from CHIKV	0	Assumed
λ_3	Rate of mosquitoes infected to CHIKV	0.5	Estimated
υ_1	Recovery rate of COVID-19 and CHIKV sequentially	0	Assumed
υ_2	Recovery rate of CHIKV and COVID-19 sequentially	0	Assumed
υ_12	Recovery rate of both COVID-19 and CHIKV simultaneously	0	Assumed

The classes of human population are and correspond to the number of individuals in the sixteen epidemiological classes at time . The total human population at time t, denoted by , is given by:

Moreover, the classes of mosquito vector population are and where the total mosquitoes’ population at time is given by:

The structure schematic diagram of the compartmental COVID-19- Chikungunya co-infection model is illustrated in Fig. 1.

Fig. 1

Schematic diagram of the compartmental Chikungunya and COVID-19 co-infection model

The proposed co-infection model is described by the following dynamic differential equations system.

Equation (1) defines all susceptible humans’ populations for both COVID-19 and Chikungunya infections. We will consider that the total population in Yemen, including those individuals recovered from any infectious, are susceptible to infection with COVID-19 or Chikungunya infections, or both, in varying proportions. Four times-variant variables are suggested to reduce the susceptible individuals: social distancing, quarantining, personal protection and rate of pesticides that are used to kill mosquitoes. and which correspond to the social measures taken during the course of the COVID-19 pandemic. refers to the control rate of social distancing of the exposed individuals. refers to the control rate of the quarantine of the infected individuals. On the other hand, and correspond to the measures taken during the course of the Chikungunya pandemic. The time-dependent input measures the personal protection from CHIKV through education and awareness from the risk of CHIKV and how to protect yourself and your family using mosquitoes’ nets, wearing long clothes and not exposing the body to mosquitoes’ bites. refers to the control rate of using the insecticides and larvicides to fight mosquitoes. Equation (2) describes the individuals exposed to COVID-19. This epidemiological class includes a fraction of susceptible individuals minus infected individuals with COVID-19, irrespective of the fact that they have symptoms or not. Equation (3) describes the individuals infected with COVID-19 that are asymptomatic or unconfirmed. The individuals for this epidemiological class come from the exposed class at a rate of . This class will be left, either by confirmation of COVID-19 infection through a screening test, by direct recovery, by direct perished due to COVID-19 infection, or by co-infection at a rate of . This epidemiological class contributes to spreading infection in a latent way, so our aim in this work is to reduce the number of the individuals of this class by increasing the level of test kits for the largest possible rate of the exposed individuals. The screening level is measured by a time-dependent parameter . Equation (4) describes the confirmed infected individuals who have been tested. The individuals of this class come directly from the exposed class at a rate of or by confirming the asymptomatic individuals at a rate of . This epidemiological class will be left either for recovery at a rate of , for death at a rate of , or for co-infection at a rate of . Equation (5) describes the individuals who have recovered from COVID-19 with or without symptoms. This class is increased by the recovery of the individuals affected with COVID-19 of both classes at a rate of . Also, this class is decreased by infection with Chikungunya at a rate or again infection by COVID-19 at a rate . Equation (6) describes the individuals perished due to COVID-19 at a rate of . Equation (7) represents the individuals who have recovered from COVID-19 but are still susceptible to Chikungunya. This class is increased by recovering the individuals infected with COVID-19 and still susceptible to CHIKV at a rate of . It is decreased by recovery from both diseases at a rate of . Equation (8) characterizes the individuals infected with both COVID-19 and Chikungunya diseases, either consecutively or simultaneously. The inputs of this class come from either the infection of susceptible individuals or if a patient with one infection becomes infected of the other infection at rates respectively. It is decreased by the individuals’ recovered at a rate of or individuals’ death at a rate of Equation (9) describes the individuals recovered from both COVID-19 and Chikungunya diseases, either consecutively or simultaneously. It is increased when infected individuals recover from Chikungunya, COVID-19, or both at rates of , respectively. It is decreased by natural death or by being susceptible to both infections at a rate of Equation (10) describes the individuals perished due to both COVID-19 and Chikungunya diseases at a rate of .

Equation (11) defines the individuals exposed to Chikungunya. This epidemiological class includes a fraction of susceptible individuals minus infected individuals with Chikungunya, with or without symptoms. The time-dynamic variables and are being entered to control the number of the exposed individuals to Chikungunya. Equation (12) describes the individuals infected with Chikungunya who are asymptomatic or unconfirmed. The individuals for this epidemiological class come from the exposed class at a rate of . Leaving this class is carried out with four possibilities, either by confirming the individuals infected with CHIKV by testing, by direct recovery, by direct death due to Chikungunya, or by co-infection. Equation (13) describes the individuals infected with Chikungunya. The individuals of this class come directly from the exposed class at a rate of . Leaving this class also carried out three possibilities, either by recovery at a rate of , or death at a rate of , or for co-infection at a rate of . Equation (14) defines the individuals recovered from Chikungunya at a rate of . This epidemiological class is decreased by returning to being susceptible to coinfection or COVID-19 at rates of , respectively. Equation (15) describes individuals’ deaths due to Chikungunya at a rate of . Equation (16) describes the individuals recovered from Chikungunya but are still susceptible to COVID-19 at a rate of .

On the other hand, the last three equations describe mosquito vectors. Equation (17) determines all the susceptible mosquito populations. Equation (18) describes the exposed mosquito. This epidemiological class increases when the mosquito bites a person infected with Chikungunya at a rate of . Equation (19) describes all the infected mosquito populations. This class is the most dangerous class in the spread of the virus among humans. However, time-dynamic variables and are used to control of the number of the mosquitoes in these classes. controls the number of the individuals infected with Chikungunya and thus reduces the number of the infected mosquitoes. is used to eliminate mosquito colonies or larvae and eggs as a whole by using pesticides as well as not to create a fertile environment for mosquitoes to breed. Overall, the parameters in this work are divided into three categories. The first category of parameters includes time-varying inputs that reflect levels of measures taken during the course of the COVID-19 and Chikungunya pandemics. Social distancing, quarantining, COVID-19 test kits, personal protection from Chikungunya using mosquito nets, and wearing long-sleeved shirts and pants as well as insecticides and larvicides to fight mosquitoes are the five parameters in this category. The second category of parameters includes values that need to be estimated from real data collected from Yemen such as transmission rate, recovered rate, and death rate. The third category of parameters includes values taken from the literature.

Particle swarm optimization

Particle swarm optimization is a metaheuristic algorithm based on the principles of behavior of flocks of birds, fish schooling, or insect swarms for finding food positions. This algorithm was developed by Kennedy and Eberhart (1995). The optimal solution of the given function is the swarm position that can be calculated using the following equation:where is the particle position, and is the velocity of the particle that are given in the following equation:where is the inertia weight, are coefficients that indicate the swarm’s ability level to the cognitive of personal and social successes,, is the personal best solution of swarm position, and is the global best solution for all swarms and for all iterations as well. Table 3 report the parameters of PSO that are used in this work. Besides, the PSO pseudo-code is expressed in detail by the following algorithm:

Table 3

PSO configuration parameters

Parameter	NP	Max_iter	ω	c_1	c_2
Value	1000	200	0.75	1.5	2

Parameter estimation

In this section, the parameters estimation problem is solved by using PSO. The real data on COVID-19 and Chikungunya outbreaks in Yemen during the period March 1, 2020 to May 30, 2020 are used to solve the parameters estimation problem. The COVID-19 data were taken from the Center for Systems Science and Engineering (CSSE) at Johns Hopkins University (https://github.com/CSSEGISandData/COVID-19). The Chikungunya data were obtained from health officials in Yemen. Since the Chikungunya data in Yemen were recorded weekly, we considered the data for both diseases weekly, not daily.

The mean squared error (MSE) between the historical data and the predicted data is minimized using the PSO, where the Eq. (22) is used as the fitness function in PSO.where represents normalization weights and represent the number of each epidemiological class and its predictions, respectively. Table 1 reports the mean and standard deviation (Std.) of the estimated parameters for Yemen during the period (March 2020–May 2020).

For brevity, only one of the experiments, which were carried out during solving the parameters estimation problem, is discussed and shown in Fig. 2. The first panel displays simulations for all human classes of historical data and the data calculated using the system of Eqs. (1)–(19). The second panel shows the estimated trajectories of time-dynamic functions, and the rest of the panels illustrate the infected, death, and recovered for both the COVID-19 and CHIKV, while the last panel shows the trajectory of the mosquito vector obtained from simulating the system (1)–(19). It’s clear from the panel of infected with COVID-19 that the end of the epidemic is not coming soon, while the CHIKV epidemic curve will come back down faster as shown in the panel of infected with CHIKV. Also, from the same panel, the curve fitting does not exactly match the historical data, but the predicted data are still logical. We can see from the panels of the death and recovered individuals for both infections the cumulative numbers continue to increase over time. Besides, the panel of the mosquito populations continues to decrease by the passage of time, which explains the effectiveness of pesticides and the different seasons as well. On the other hand, we can clearly see in panels 3 to 9 that the curves obtained using PSO are more fitness than those obtained using Pyomo optimization modelling. This confirms that the proposed algorithm is effective.

Fig. 2

The result of simulation for all human classes, the estimated trajectories of the time-dynamic functions, and analysis and week-wise prediction using the proposed system. The real data represented by the black dotted, the red dotted lines represent the fitting line from the proposed system, the predicted data using PSO represented by the red solid lines, and the predicted data using Pyomo optimization modelling represented by the blue solid lines (color figure online)

Optimal control problem

The basic principle of the proposed optimal control model is to curb epidemic outbreaks by breaking the link between the epidemic vector with individuals through time-dynamic interventions, so that the number of infected individuals is reduced as well as the detection of infected individuals without symptoms by increasing test kits and then isolating and treating affected individuals. Moreover, minimizing the cost associated with optimal control is another goal for decision makers. Five control functions, i.e. , and , are considered in this work. The first control function represents the control rate of the social distancing. This variable can be applied to the whole society by raising awareness and obliging individuals to wear masks, personal protection, hygiene and spacing between individuals at least a meter distance, and this, in turn, contributes greatly to breaking the chain of transmission between humans for COVID-19 directly and also indirectly contributes to limiting the transmission of the Chikungunya. By increasing the control rate of this variable, the number of people who are at risk of infection decreases. The second control function indicates quarantine and insolation rates. It is easy to apply isolation of individuals who are confirmed to be infected, but it is difficult to implement quarantine on a large scale, especially in Yemen, where most of the population live on daily wages. Hence, comes the importance of the third variable , which is to increase the number of tests. It is possible to discover infected individuals without symptoms and work to isolate and treat them, which contributes to curbing the spread of COVID-19. The fourth variable refers to the rate of personal protection through protective clothing and the rate of awareness and education of the risks of infectious diseases by mosquitoes including CHIKV. Increasing the rate contributes to reducing the infected individuals with CHIKV. The fifth variable indicates the importance of purity of the surrounding environment, the seriousness of the swamps, and the fertile places for mosquito breeding. Increasing the use of pesticides contributes to killing mosquitoes and reduces the transmission rate of CHIKV to humans.

Two objectives are proposed in this work. The first one is to minimize the infected individuals with each infection or both. The second one is to minimize the cost associated with optimal control strategies.

(1–19)

The objective function in Eq. (23) minimizes the number of the individuals infected with COVID-19, CHIKV, and both. The objective function in Eq. (24) represents the cost associated with the optimal control strategies. The dynamic differential Eqs. (1)–(19) are also included as constraints of the optimal control model. Constraint (25) defines the maximum limit of the infected individuals for each infection and the peak value of the infected mosquitoes as well. The constraint bounds (26)–(30) define the ranges of the time-dynamic parameters, where and are weight coefficients. The optimal control model is assumed during the full-time horizon , where is 52 weeks.

The summarized steps of the present work are given as the following:

There are two approaches to solve the above optimal control model, namely, indirect methods such as mathematical analysis and direct methods such as PSO. In this work, we use the PSO to solve the optimal control model. Before we solve it using PSO, we briefly illustrate the mathematical analysis to prove the theorem of the solution.

In mathematical (indirect) methods, Pontryagin’s maximum principle is used to solve the optimal control model, through driving the Hamiltonian function and defining the necessary conditions for the optimal control. The optimal solution of the control model is in whichwhere and , , , , and

The Hamiltonian H is defined as:where is an adjoint variables, which corresponds to each that represents the right-hand side of the system Eqs. (1)–(19).where , and are adjoint variables or co-state variables.

Theorem

Given the optimal control and solutions , and of the corresponding state system (1)–(19) that minimize over There exist adjoint variables , and such thatwhere , and with transversality conditionsandwhere

Proof.

See Appendix A.

Numerical simulations

In this section, we investigated and analyzed the impact of the five dynamic control measures, the peak limits, and the coefficient weights of control measures in the objective function, which will determine the best strategy to curb the rapid spread of both Chikungunya and COVID-19 in Yemen in 2020. The initial conditions that are used in all the strategies of the present model are mentioned in Table 1. Besides, the parameters that are needed for estimation are reported in Table 4, and the other parameters are reported in Table 2. We assume the weights A, B, C, D, and E are equal one all simulations except when discussing the impact of their changes. We also used the same limit peaks for all simulations except when discussing the impact of limit peaks constraints. The weighted method was employed to convert the bi-objective optimization model into a single-objective optimization model. In addition, the dynamic differential equations system (1–19) was discretized using the orthogonal collocation method with considering the time domain which consists of weekly finite elements. Also, we used the penalty function to deal with the constraints. Then, the new optimization model was be solved using PSO. It is worth mentioning here that the proposed algorithm was implemented in Python 3.7. Plus, the figures were omitted for each case, and only the figures that are related to the comparisons were shown, for the sake of brevity. The strategies are discussed in this work as the following.

Table 4

The estimated parameter values using PSO

Parameters	α_c	α_1	α_2	α_3	φ_1	φ_2	λ_3
Mean	0.0209	0.0900	0.0780	0.3574	0.55	0.5000	0.5000
Std	0.0630	0.1449	0.1839	0.0701	0.2415	0.0000	0.0000

Strategy 1: In this strategy, we discussed three levels of controls, which are max control, partial control and no control (baseline). In the max control level, the range of time-dynamic variables were , , , , , while were , , , , in the partial control level. Here, the range of the inequalities for the and means that a fraction of the population who implement these rules lies between the values of the left and right sides of the inequalities. For , the rate of the availability of screening tests for the individuals exposed with COVID-19 is located between the boundaries of inequality. Also, the rate of using the insecticide lies between the values of the left and right sides of the inequality. We also set the coefficients of the time-dynamic variables equals zeros in the objective function in the baseline (without control) level. Figure 3 illustrates the compression of the max control, partial control, and baseline levels for both infected, death, the individuals recovered from COVID-19, infected, death, individuals recovered from CHIKV and infected mosques vector throughout the simulated time. We can see that max control comes in the first rank. It is more effective to reduce the number of infected individuals and to decrease the time-span of both epidemics, then partial control comes in the second rank, and no-control (without control) level comes in the last rank. It is normal to expect this ranking due to the impacts of interventions through applying social distancing, quarantine, detection of latent cases by providing test kits, personal protection and use of insecticides, which in turn reduces the number of individuals with both diseases. On the other hand, the cost associated with max optimal control was higher than the cost associated with other optimal controls.

Fig. 3

Comparison of the simulation results of the control levels of strategy 1

Strategy 2: we concentrated on the effect of the limit of the peaks in the constraint (25). Two cases were investigated, i.e. small and large limit peaks. For infected humans, the small limit of peaks was while the large limit of the peaks equals the total of population in Yemen . The small and large limit peaks for the infected mosquitos were and respectively. We also used the same range of time-dynamic variables in all cases. As it is shown in Fig. 4, the right-hand side values of limit constraints affect all epidemiological classes. The number of individuals was less when the limit constraints were small. We conclude that the number of individuals can be reduced by decreasing the right-hand side values of limit constraints and vice versa, i.e., we can increase the number of individuals for all classes by neglecting limit constraints. However, the cost associated with small limit peaks was higher than the cost associated with large limit peaks.

Fig. 4

Comparison of the simulation results of two cases of strategy 2

Strategy 3: In this strategy, we focused on the effect of the combination of the coefficients of intervention variables in the second objective function, and whether or not they will be reflected in the number of individuals in the epidemiological classes. We investigated that in two cases. In the first case, all the coefficients were equal one , and in the another case, , , , , and Figure 5 shows the simulation of the two cases throughout the time horizon. It is clear that the first case is superior to the other case in all epidemiological classes. It is also known that the coefficients of the intervention variables are the relative cost associated with these interventions. Therefore, the more coefficient values are, the higher the cost value is. Thus, as expected, the cost of the first case (with coefficients are one) is higher than the cost of the other case (with coefficients are less than one).

Fig. 5

Comparison of the simulation results of two cases of strategy 3

Figure 6 shows the trade-off curve of both functions that can be obtained when using different weights,

Fig. 6

Trade-off curve of both functions for all cases

It can be concluded that we need some interventions to curb the spread of the epidemics. These interventions depend on the awareness of the community as well as its economic capacity and response to implement them. It should be borne in mind that some interventions cannot be applied such as imposing quarantine and total closure, especially in poor societies that depend on daily income. Hence, we recommend that partial control is the closest and most appropriate solution that can be applicable in the Yemen case, as it will be more harmonious and compatible with other economic, social, cultural factors outside the framework of this study.

Conclusions

In this study, an optimal control model with strategies of COVID-19 and CHIKV co-infection was investigated. First, the state-space model was formulated, and then the parameters were estimated for the outbreaks in Yemen for the period (March 2020–May 2020). Critical epidemiological interventions were identified, and their impact on the optimal control model was investigated. The bi-objective optimal control model aims to minimize the infected individuals and the cost associated with the corresponding control. PSO algorithm was used to solve the parameter estimation problem and the optimal control model as well. Numerical simulations clearly show the importance of the present optimal control model for controlling epidemics, which may give the government highlights to select the suitable strategies that curb epidemics while simultaneously preventing economic collapse. As a future suggestion, the introduction of uncertainty to some parameters of the model can be considered.