Parameter Estimation of Compartmental Epidemiological Model Using Harmony Search Algorithm and Its Variants

Epidemiological models play a vital role in understanding the spread and severity of a pandemic of infectious disease, such as the COVID-19 global pandemic. The mathematical modeling of infectious diseases in the form of compartmental models are often employed in studying the probable outbreak growth. Such models heavily rely on a good estimation of the epidemiological parameters for simulating the outbreak trajectory. In this paper, the parameter estimation is formulated as an optimization problem and a metaheuristic algorithm is applied, namely Harmony Search (HS), in order to obtain the optimized epidemiological parameters. The application of HS in epidemiological modeling is demonstrated by implementing ten variants of HS algorithm on five COVID-19 data sets that were calibrated with the prototypical Susceptible-Infectious-Removed (SIR) compartmental model. Computational experiments indicated the ability of HS to be successfully applied to epidemiological modeling and as an efficacious estimator for the model parameters. In essence, HS is proposed as a potential alternative estimation tool for parameters of interest in compartmental epidemiological models.


Introduction
Epidemiological models play a vital role in understanding the spread and severity of a pandemic (or epidemic) of infectious diseases [1]. During an outbreak of an infectious disease, it is crucial to simulate the potential outbreak growth for planning the outbreak control measures in order to provide useful insights into measurable outcome of existing interventions, predictions of subsequent growth, risk estimations, and guiding alternative interventions [2][3][4]. Epidemiological constraints, such as delays in symptom appearance (due to incubation period) and positive test confirmation (due to limited testing and detection resources), may limit the real-time use of epidemiological models [5,6]. In order to overcome such constraints, mathematical modeling of infectious diseases was employed in epidemiology, as recognized by WHO [7] and proven to be effective [8,9]. Compartmental modeling as a class of mathematical modeling was widely applied to infectious diseases modeling [10]. The Susceptible-Infectious-Removed (SIR) model is the first compartmental modeling approach for simulating the probable outbreak trajectory [11]. Besides the standard model, various extensions of SIR were developed in the recent past, mostly by including additional compartments, such as the Susceptible-Exposed-Infectious-Removed (SEIR) model. During the current COVID-19 global pandemic crisis, many studies (e.g., [12][13][14][15][16][17][18]) applied the SIR model and its extensions to analyze the dynamics of the disease. A recent review on SIR family of models used for studying, predicting, and managing COVID-19 is presented in [19].

The Susceptible-Infectious-Removed (SIR) Model
Despite the development of various extensions of the SIR model, the standard SIR model remains the preferred first approach for analyzing the spreading of an infectious disease (especially in the beginning or first phase) and it is reasonably predictive [11]. The SIR model splits a given population N into three compartments (non-intersecting classes) at any given time t (measured in days) namely, (i) susceptible (not yet infectious and disease free individuals at t) denoted S t , (ii) infectious (confirmed or isolated individuals) denoted I t , and (iii) removed (no longer infectious or dead) denoted R t . The number of individuals in each compartment vary over time. In general, the dynamics is described with a large number of susceptible individuals at the beginning, since the entire population that is not infected is considered to be susceptible, while infectious individuals remain low at the beginning of a pandemic. At subsequent times, the number of infectious individuals will increase, the number of removed individuals will begin to gradually increase, and the number of susceptible individuals will decrease. Finally, towards the end of a pandemic, the number of removed individuals will increase, the number of infectious individuals will decrease gradually, while the number of susceptible individuals will remain the lowest. The disease dynamics according to SIR model can be visualized, as in Figure 1. The variation rate over time in each compartment is modeled using a system of non-linear ordinary differential equations (ODEs), (1) Figure 1. Infectious disease dynamics according to SIR model. Adapted from [11].
The primary assumption of SIR model is that the population is closed and fixed, and it is the sum of individuals in all of t compartments i.e., N = S + I + R for all t. The epidemiological parameters of interest are, i) the transmission rate β and ii) the recovery rate γ. Accordingly, the average transmission (from an infectious individual to a susceptible individual) period is 1/β days and the average infectious period is 1/γ days. Estimation of the SIR parameters is a critical task, since there is no closed-form analytical solution to the SIR model. Numerical approximation methods, such the Runge-Kutta methods, are often employed in order to solve the SIR model with estimated parameters. Thus, the quality of the simulation heavily relies on the estimates of the epidemiological parameters. Apart from that, a good estimate of the parameters is crucial in assessing the transmissibility of an infectious diseases in real-time through the effective reproduction number R t . One of the approaches for inferring R t is by using compartmental epidemiological models, in which R t is treated at the deflation of the basic reproduction number R 0 [20], which can be estimated using the relationship among the epidemiological parameters [21] that is given by:

Determining the SIR Parameters
The surging interest in infectious diseases modeling due to the COVID-19 pandemic has led to a different approach of compartmental modeling from the usual mathematical modeling strategy to statistical modeling strategy [22]. In the latter strategy, the model parameters are estimated instead of being specified by or adapted from certain subjective prior information as in the former strategy (usually from previous studies, previous pandemic or values estimated by WHO). The estimation of the SIR model parameters is essentially an optimization problem attempting to find a model that best fits the data. The estimated parameters are then used along with any numerical approximation methods to obtain the simulated SIR compartments' trajectories. A least squares loss function in terms of the simulated and observed trajectories, such as the Sum of Squared Error (SSE), is usually applied to quantify the discrepancy that arises from the simulation. Hence, the objective of the optimization is to minimize the loss function by estimating the parameters that lead to the best fit curve through any standard optimization techniques.
Bayesian approach is becoming the popular optimization and estimation tool, as evidenced by the studies concerning COVID-19 pandemic. This approach is commonly employed by calibrating the available data while using Markov Chain Monte Carlo (MCMC) method with Metropolis-Hastings (MH) algorithm sampling as implemented in [6,[12][13][14]18,23] to obtain posterior estimates and credible intervals of the epidemiological parameters. Although the estimation of the model parameters is an obvious optimization problem, the metaheuristics family of optimization techniques received very little attention in epidemiological modeling. The first metaheuristic algorithm used in estimating the parameters in ODEs in general, and infectious diseases specifically, is the Genetic Algorithm (GA) (e.g., [24][25][26]). Recently, the Particle Swarm Optimization (PSO) algorithm is implemented in order to estimate the parameters in ODEs governing the SIR model variants, as presented in [27]. As for COVID-19 pandemic, very few studies applied metaheuristic algorithms for estimating the epidemiological model parameters, such as GA in [28], PSO in [29,30], Stochastic Fractal Search in [31], Marine Predators Algorithm in [32], and Flower Pollination Algorithm with Salp Swarm Algorithm in [33]. In this regard, we are interested in applying a metaheuristic algorithm, namely the Harmony Search and its variants, to the optimization problem of estimating the SIR model parameters.
This paper is organized, as follows. Section 2 details the Harmony Search algorithms that were used in this study. Section 3 presents the experimental setup of estimating the epidemiological parameters of SIR model while using Harmony Search algorithms. Section 4 provides the simulation results with detailed discussions. Finally, Section 5 provides the conclusion and possible future works of this study.

Harmony Search Algorithm and Selected Variants
The Harmony Search (HS) algorithm [34] is a well-known population-based metaheuristic algorithm. The optimization process in HS is a mimicry of the underpinning principles of jazz music orchestra, where musicians attain a pleasant harmony through several improvisation steps. HS has been successfully applied to a wide variety of real-world optimization problems, such as system reliability, robot path planning, renewable energy systems, hyper-parameter tuning of deep neural networks, intelligent manufacturing, and credit scoring (see [35,36]); university timetabling, structural design, water distribution, and supply chain management (see [37,38]); and, music composition, Sudoku puzzle solving, tour planning, web page clustering, vehicle routing, dam scheduling, groundwater modeling, soil stability analysis, ecological conservation, heat exchanger design, transportation energy modeling, satellite heat pipe design, medical physics, medical imaging, RNA structure prediction, and image segmentation (see [39], among others). Besides that, the implementation of HS in various parameter estimation studies indicated the potentiality of HS as an effective parameter estimation tool. Some of the notable parameter estimation problems that applied HS include parameter estimation of the nonlinear Muskingum model [40,41], parameter estimation in vapor-liquid equilibrium modeling [42], parameter estimation in an electrochemical lithium-ion battery model [43], parameter identification of synthetic gene networks [44], and design storm estimation from probability distribution models [45]. In addition, HS was also successfully employed in human activity pattern modeling, such as disease spread and disaster response [46]. Hereinafter, the term HS represents the family of HS variants, while the standard HS is denoted SHS. The five primary steps in SHS, as outlined in [47], are as follows: 1. Initialization of HS parameters and the objective function. The control parameters are Harmony Memory Size HMS, Harmony Memory Considering Rate HMCR, Pitch-Adjusting Rate PAR, bandwidth (now known as fret width) BW, and maximum improvisations MaxImp. If f (·) is the objective function with n decision variables x = (x 1 , ..., x n ) in the range (LB i , UB i ), then the continuous optimization problem can be formulated as follows: 2. Initialization of Harmony Memory (HM). HM is a HMS × n dimensional matrix that consists of randomly generated harmonies (candidate solutions) from the uniform distribution U(0, 1) within the ranges of the decision variables. In general, it is more convenient to represent HM as an augmented matrix of order HMS × (n + 1), as follows: where each row of the HM represents the solution vector x j , (j = 1, 2, ..., HMS) in the first n−columns, followed by the fitness that is generated from the solution vector, f (x j ). It is also common to have HM as a sorted matrix in the ascending order of the fitness value (the last column of HM). 3. Improvisation. Improvisation is performed to generate a new harmony by exploring and exploiting the search space. Thus, a new harmony is randomly selected from the HM with a probability of HMCR or it is randomly generated outside of HM with probability 1 − HMCR. If a new harmony is obtained from HM, then the harmony may be improvised by adjusting the harmony with neighborhood values that are based on BW with probability of PAR or remain as is with probability 1 − PAR. Note that HMCR is inversely proportional to the explorative power of different search spaces, while PAR is directly proportional to the exploitation power of local search space. 4. Update HM. New harmony from Step 3 is evaluated with the objective function to obtain the new fitness. If the new fitness is lower than the worst fitness, then the worst solution in HM will be replaced with the new harmony. 5. Termination. Repeat Steps 3 and 4 until MaxImp has been reached or other termination criteria are satisfied.
The complete details of SHS are provided as in Algorithm 1. SHS is designed with fewer mathematical operations and it is relatively easy to code, but easily applicable to a wide variety of optimization problems. The advantages of HS are discussed in [37]. Over the two decades since the inception of HS, many HS variants have been developed to date. A majority of the variants modifies the improvisation procedures either by internal modification or hybridization with other heuristics. Some of the recent comprehensive reviews on HS variants are presented in [37,39,[47][48][49]. In this study, only the internally modified HS variants were considered and selected on a minimal parameter setting requirement basis. for each i ∈ [1, n] do 6: if U(0, 1) ≤ HMCR then 7: x i = x j i where j ∼ U(1, HMS) % memory consideration 8: if U(0, 1) ≤ PAR then 9: x i = x i + (2 × U(0, 1) − 1) × BW % pitch adjustment 10: end if 11: else 12:

Improved Harmony Search (IHS)
The Improved Harmony Search (IHS) [50] is the prototypical HS variant. While still requiring to fine-tune the HMCR parameter, the parameters PAR and BW are made dynamic in this variant with the introduction of the PARmax (BWmax) and PARmin (BWmin) iterative parameters. Although IHS has shown better performance than SHS, this variant increased the burdensome process of setting suitable values for four parameters instead of just two in SHS [51]. Nevertheless, IHS is considered to be a breakthrough that paved the way for the development of various HS variants to date. PAR and BW are adjusted at each iteration while using: The computational procedure of IHS is provided, as in Algorithm 2. for each i ∈ [1, n] do 6: if U(0, 1) ≤ HMCR then 7: if U(0, 1) ≤ PAR t then 9: end if 11: else 12: A year later, a new variant that was inspired by the swarm intelligence concept of the PSO was introduced by the principal developer of IHS. This variant is known as Global Harmony Search (GHS) [52] and aims to mimic the best harmony in the HM. The parameter PAR is adapted from IHS, while BW is removed. Thus, the pitch adjustment step in SHS is replaced with a random selection of best harmony of any decision variable from the HM, as follows: In general, GHS is claimed to perform better than SHS and IHS, especially in highdimensional optimization problems. However, [53] asserted that GHS has flaws that will cause premature convergence and the name of this variant is also said to be misleading. The most serious flaw, as noted by [51], is the frequent generation of infeasible new harmonies, whenever the upper and lower bounds of each decision variable are not identical in the given optimization problem. Hence, GHS is not considered in this paper.

Novel Global Harmony Search (NGHS)
The Novel Global Harmony Search (NGHS) [54] adapts the swarm intelligence concept of the PSO algorithm in the improvisation step of the SHS. This approach enables the new harmony to mimic the global-best harmony in the HM. Thus, the HMCR and PAR parameters are removed and the improvisation is only dependent on the best and worst harmonies in HM. The random generation of harmony is, in fact, analogous to SHS, and the only difference is that, instead of randomly selecting a harmony with the probability of 1 − HMCR, NGHS performs genetic mutation with the probability of Pm, based on the ideas from evolutionary algorithms. NGHS is presented, as in Algorithm 3.

Self-Adaptive Global Best Harmony Search (SGHS)
The Self-Adaptive Global Best Harmony Search (SGHS) [55] aims to improve the GHS [52] in terms of avoiding getting trapped at local optima. In this approach, HMCR and PAR are dynamically adjusted to a suitable range after a number of iterations by tracking their previous values that allowed for the replacement of new harmony in HM. Further, HMCR and PAR are assumed to be normally distributed where HMCR ∼ N(HMCRm, 0.01), HMCR ∈ [0.9, 1.0] and PAR ∼ N(PARm, 0.05), PAR ∈ [0.0, 1.0]. The initial values of HMCRm and PARm are set at 0.98 and 0.9, respectively. Subsequently, SGHS begins with HMCR and PAR values being generated from the Normal distribution. During each iteration, the values of HMCR and PAR that corresponds to a replacement of new harmony in HM is recorded until a number of solutions are generated within the specified learning period LP. Once LP is reached, the recorded HMCR and PAR values in previous iterations are averaged to obtain new HMCRm and PARm to be used in upcoming iterations. This process is repeated until the termination criterion is satisfied. As for BW, the values are dynamically adapted, as follows: SGHS is outlined, as in Algorithm 4.

Intelligent Tuned Harmony Search (ITHS)
Based on the idea of sub-population approach for optimization (such as [56]) and decision-making theory, the Intelligent Tuned Harmony Search (ITHS) [57] attempts to intelligently control the exploration and exploitation in HS based on consciousness or previous experience. This approach begins by assigning the x best as the leader of the whole population. The leader divides the HM into two groups (sub-populations), say Group I and Group II, in order to achieve a good balance between exploration and exploitation. Group I consists of the harmonies with fitness less than or equal to average fitness and Group II vice-versa. In that sense, Group I will undergo both exploration and exploitation stages, while Group II will only undergo exploration stage. ITHS uses the same adaptation of dynamic PAR as the Self-Adaptive Harmony Search (SAHS) [53] that is given by: Algorithm 5 presents the computational steps of ITHS.

Novel Self-Adaptive Harmony Search (NSHS)
The Novel Self-Adaptive Harmony Search (NSHS) is a HS variant that was developed by [51], being inspired by the defects that the creator found in SHS and other variants, namely IHS [50], GHS [52], SAHS [53], Dynamic Local Harmony Search (DLHS) [56], and SGHS [55]. In NSHS, the HMCR parameter is constructed based on the dimension of the optimization problem to be solved, With reference to Equation (10), HMCR is set to be directly proportional to n, in order to use the HM more frequently, and it lies in the interval (0.5, 1). The parameter PAR is removed. Furthermore, a dynamic fine-tuned BW is introduced and it depends on the standard deviation S of the objective function, f std = S( f (x j )), ∀j = 1, 2, ..., HMS. BW diminishes in stages according to the iteration number t, while increasing with a larger range of decision variables. The improvisation step in NSHS generates a new harmony within the narrow range of [x worst i , x best i ] based on conditions of HMCR and f std. Algorithm 6 outlines the computational procedure of NSHS.

Global Dynamic Harmony Search (GDHS)
Based on IHS [50], the Global Dynamic Harmony Search (GDHS) [58] further improves the improvisation step of HS with dynamic parameters, as well as dynamic upper and lower bounds of the decision variables. The iterative values of HMCR and PAR are made to be both decreasing and increasing in the search of global optima, as given by: For BW, the dynamic adjustment is adapted from IHS, but with a few modifications, as follows: and, based on Equation (6) from IHS, the equation reduces to, Next, a correction coefficient, coe f , is introduced at each iteration by: where j is the index of the selected harmony in the memory consideration step. Finally, the dynamic lower and upper bounds are obtained for the random selection step. Algorithm 7 provides the computational steps of GDHS. if U(0, 1) ≤ HMCR t 7: compute coe f t using (15) 10:  The Parameter Adaptive Harmony Search (PAHS) [59] focused on the modification of the improvisation step of IHS [50]. The dynamic values of HMCR, PAR, and BW are generated during each iteration to ensure the global optima is achieved. The authors explored four different combinations of dynamic HMCR and PAR iterative values i.e., (i) linear HMCR and PAR; (ii) exponential HMCR and linear PAR; (iii) linear HMCR and exponential PAR; and, (iv) exponential HMCR and PAR. Through computational experiments, it was concluded that linear HMCR and exponential PAR yields the best performance. In this way, HMCR gradually increases, while PAR exponentially decreases with respect to the iterations. HMCR and PAR at each iteration are computed using: whereas, BW is adapted from IHS as it is. PAHS further aggravates the difficulty of finding suitable values as there are six parameters to be set now, rather than only four in IHS. PAHS is detailed in Algorithm 8.

Algorithm 8: Parameter Adaptive Harmony Search
1: Set HMS, HMCR using (16), PAR using (17), BW using (6) if U(0, 1) ≤ HMCR t then 7: if U(0, 1) ≤ PAR t then 9: end if 11: else 12: The Enhanced Self-Adaptive Global Best Harmony Search (ESHS) [60] is one of the recent HS variants that retains the simplicity and distinctive framework of SHS. In order to eliminate the troublesome parameter fine-tuning process in SHS, a new parameter settingfree strategy is proposed without requiring any extra statistic and external archive. HMCR is dynamically obtained at each iteration as a random normal number, while PAR is given by: and BW is defined by: ESHS employs the Gaussian mutation technique in contrary to the uniform randomization in SHS. Gaussian mutation is claimed to be more efficient in exploring the global optimum solution as compared to uniform randomization. Thus, the random selection in ESHS is performed while using: with the probability of 1 − HMCR. ESHS is detailed in Algorithm 9.

14:
if (U(0, 1) ≤ PAR) then 15: x i = x i + U(0, 1) × BW 16: Despite being successfully applied in various fields, to the best of our knowledge, HS algorithms have yet to be applied in epidemiology, particularly in epidemiological modeling. Thus, in this study, ten variants of HS algorithm is proposed to be applied in order to estimate the epidemiological parameters of interest in the prototypical compartmental epidemiological SIR model and compare the estimation performance of each algorithm.

Estimating the Epidemiological Parameters of SIR Model
HS algorithms are applied to the SIR parameter estimation problem while using the cumulative infectious cases (total cases) of the COVID-19 pandemic as the use-case. The optimized (final) values of parameters β and γ are estimated by calibrating the available COVID-19 data and the SIR model with the generated harmonies (candidate solutions) from HM.

COVID-19 Data Sets
The authors obtain the time series of cumulative infectious cases per day for five countries, namely the United States of America (USA), France (FR), South Korea (SK), Ireland (IR), and Singapore (SG), by web scraping the figures in https://www.worldometers.info/coronavirus (which gathers data from various reliable sources including European Centre for Disease Prevention & Control and Johns Hopkins University & Medicine Coronavirus Resource Center). The data collected are for a period of 240 days, beginning from the first day of the outbreak in each country. Next, the first 220 days data is used for calibration to obtain the estimates of the parameters. The remaining 20 days of data are used for validation against the projection of simulation produced while using the HS optimized parameter values.

SIR Model Setup
The SIR compartments are initialized (at time t = 0) with initial conditions I 0 and R 0 , according to the actual number of infectious and removed cases, respectively, on the first day of the outbreak in each country. Meanwhile, S 0 is set as the remaining individuals in the population N who are yet to be infected and removed, i.e., S 0 = N − I 0 − R 0 . Table 1 provides the modeling initialization.

SIR Parameters Estimation as Optimization Problem
The estimation of SIR parameters is formulated as an optimization problem with decision variables x = {x 1 , x 2 }, where β = x 1 and γ = x 2 . The upper and lower bounds of the decision variables are set according to the complete possible ranges of β and γ, which are essentially the epidemiological dynamics' rates and, hence, they share the same bounds, i.e., β, γ ∈ (0, 1). The objective function is formulated, as follows: 1.
Let C T = ∑ t=T t=0 I t be the observed cumulative infectious cases on day t = T, (t, T ∈ [0, 219]), where I t is the observed infectious cases on day t.

2.
Let C T = ∑ t=T t=0 I t be the simulated cumulative infectious cases (rounded to the corresponding integers) on day t = T, (t, T ∈ [0, 219]), where I t is the SIR model simulated infectious cases on (rounded to the corresponding integers) day t while using x ∈ HM. 3.
Subsequently, an objective function f (·) that minimizes the SSE between C T and C T can be formulated as: A set of ten independent runs is performed for each HS algorithm while using each of the data sets. The optimization steps are detailed, as follows: 1.
The common control parameters shared among all ten algorithms are set to be identical, HMS = 30 and MaxImp = 50, 000, following the recommendation in [47]. 2.

3.
Optimization is performed for ten independent runs for each data set while using the identical HM and f using each of the algorithms (Algorithm 1 to Algorithm 10), as described in Section 2. For instance, the first run for the USA data set will be performed while using the same HM 1 and f 1 for all ten HS algorithms. The specific control parameters that are shared among some algorithms are also set to be identical as displayed in Table 2. 4. Repeat Step 3 until the runs are completed for all five data sets. The combination of parameters' values that yields the best fitness (lowest SSE) are designated as the optimized parameters.

5.
The Subsequently, {x 1 ,x 2 } are designated as the overall optimized SIR parameters with respect to each algorithm, thus x optimized = {x 1 ,x 2 }. Note that the identical common and specific control parameters, as well as the identical HM and f , are used in this study to ensure a fair comparison among the HS algorithms.

Evaluation of HS Estimation and Performance Comparison
The accuracy of the SIR parameter estimates is deduced from the fitness values (SSE), where the algorithm with lowest fitness is designated as the best performing estimator and vice-versa. The optimized epidemiological parameters x optimized are supplied to the SIR model once again to produce a projection of simulation for a period of 20 subsequent days from the end date of calibration period in order to evaluate the predictive capability of SIR model while using the parameters that were estimated from HS algorithms. Let C T be the observed cumulative infectious cases on the projection period and C T be the projected SIR model simulated cumulative infectious cases using x optimized , then the accuracy of estimation is evaluated by computing the Root Mean Squared Error (RMSE) between C T and C T in the period of projection using: The accuracy of the SIR model's projected simulation is decided based on the RMSE value, where a lower RMSE value indicates a better prediction and vice-versa. Eventually, the predictive capability indicates the parameter estimation accuracy of the HS variants. The corresponding RMSE for each algorithm within the same data set will be statistically compared. The comparison is performed while using the Friedman test to see whether there is an overall significant difference in the estimations produced by each of the algorithms. Furthermore, if the performance of SHS is found to be comparable with the rest of the HS variants, then the post-hoc procedure of Wilcoxon signed-rank tests is conductedin order tto determine whether there is any statistically significant difference between the estimation performance of SHS and the rest of the HS variants individually.

SIR Simulation Experiments and Discussion
In this section, SIR simulation experiments are performed to illustrate the ability of HS algorithms as an efficacious estimator of SIR parameters. All of the variants of HS algorithm were coded in MATLAB R2017b on a laptop computer with 2.50 GHz Intel i7-4710HQ CPU with 32 GB of RAM. The discussion is discretely presented for each data set used in this study.  Table 3 presents the optimized epidemiological parameters and corresponding fitness values (average from ten independent runs) for the USA data set. Figure 2 displays the visualization of the 220 days simulation together with the 20 days projected simulation using the optimized parameters. Visual inspection for USA data set is not informative enough as some of the lines representing the algorithms overlapped, which suggested that the estimates are very close to each other, except for NGHS, IBGHS, SHS, ITHS, ESHS, PAHS, and IHS that are visible. Simulations of each algorithms for approximately the first 100 days were indeed close to each other as well as to the observed cumulative cases. Beginning from the hundredth day of the outbreak, the simulations of each algorithm showed differences while they started deviating from the observed values, except for PAHS and IHS, which are still intact with the observed values. The projected simulation for subsequent 20 days were consistent in terms of the pattern in the calibration period. It is observed from the fitness values (SSE) in Table 3 that IHS appears to be the best performing estimator for USA data set, while NGHS is the least performing estimator.  For the USA data set, it is observed that the simulation from each HS algorithms are not very far off the observed cumulative cases, even in the projection period. Most of the simulations are similar as far as the parameters' values are concerned and they were able to approximately resemble the observed trend. Simulations of NGHS and IBGHS are the farthest deviation from the actual values. The similar behavior between these two algorithms may be due to the use of genetic mutation in the improvisation step that sets them apart from the rest of the algorithms. Table 4 provides the optimized epidemiological parameters and the corresponding fitness values (the average from 10 independent runs) for the FR data set. The visualization of the 220 days simulation together with the 20 days projected simulation using the optimized parameters are depicted in Figure 3. The simulations can be visually inspected as the lines representing the algorithms do not overlap each other except for SHS, which overlaps with GDHS. During the calibration period the simulation of IHS was similar but not so close to ESHS and PAHS. However, in the projection period, IHS's simulation seems to be closer to ESHS and PAHS. The other groups of algorithms that produced similar simulations are ITHS, NSHS, SGHS and SHS, GDHS, IBGHS. The simulations started diverging from each after day 60 of the outbreak, and gradually digressed from the observed cumulative values all the way up to the projection period. The projected simulation for subsequent 20 days depicted a further deviance of the simulations from the observed values. It is observed from the fitness values (SSE) in Table 4, that IBGHS appears to best performing estimator for FR data set, while NGHS is the least performing estimator, just as the case of USA data set.  For the FR data set, observe that the simulations are slightly far off the observed cumulative cases. The simulations were not able to accurately mimic the observed trend. NGHS's simulation was apart all the way, whereas we can infer that the simulations of IHS, ESHS, and PAHS are similar; ITHS, NSHS, and SGHS are alike; and finally, SGHS, GDHS, SHS, and IBGHS are close to each other. NGHS's simulation is again the farthest and different from the rest, which is probably due to the use of genetic mutation, which sets it apart from the rest. Genetic mutation is also used in IBGHS, but the standard PAR setting in IBGHS could have contributed to the high similarity in behavior as SHS than NGHS. Hence, for this particular data set, the PAR parameter was more influential then genetic mutation, as compared to the USA data set. Table 5 presents the optimized epidemiological parameters and the corresponding fitness values (average from 10 independent runs) for the SK data set. Figure 4 displays the visualization of the 220 days simulation, together with the 20 days projected simulation using the optimized parameters. The simulations can be well visualized as the lines representing each algorithm is distinct, except for SHS, which overlaps with ESHS. The simulations started departing from the observed cumulative cases as early as the fiftieth day of the outbreak. The band of simulations also started deviating around the same time and was gradually separated up to the projection period. All of the simulations were distinguishable, except in the case of IHS and NSHS, which were close to each other and SHS and ESHS that were similar. The projected simulation for subsequent 20 days indicated a larger deviation of the simulations from the observed values. We observe from the fitness values (SSE) in Table 5 that ITHS appears to be the best performing estimator for SK data set, while NGHS is the least performing estimator, just as the case of USA and FR data sets.   For the SK data set, observe that the simulations are quite far off the observed cumulative cases, with NGHS being the farthest and ITHS being the nearest. Yet, ITHS's simulation is still far off the observed cumulative cases comparatively. The simulations were not able to resemble the observed trend well. The simulation for SK data set is clearly not as good as USA or FR data sets. the simulations of all the algorithms were distinct except for the overlapping SHS and ESHS. NGHS's simulation is again very different from the rest, probably due to use of genetic mutation, which sets it apart from the rest. As far as this data set is concerned, it is also interesting to observe ESHS that requires zero parameter setting and uses Gaussian mutation in the place of random generation of harmony, behaves somewhat similar to SHS. Table 6 provides the optimized epidemiological parameters and the corresponding fitness values (average from 10 independent runs) for the IR data set. Figure 5 depicts the visualization of the 220 days simulation, together with the 20 days projected simulation using the optimized parameters. The visual inspection of the simulations is informative, although the lines representing most of the algorithms are quite close to each other, with only SGHS and PAHS overlapping. The simulations drifted away from the observed cumulative cases as early as before the fiftieth day of the outbreak itself. However, the simulations of SGHS, PAHS, and NSHS managed to stay intact with the observed values, even during the projection period. On the other hand, the simulations of IHS, IBGHS, and SHS were still close to the observed values and they became closer in the projection period. ITHS remained steadily far from the observed values, while ESHS, NGHS, and GDHS maintained a constant separation from the observed values. The projected simulation for subsequent 20 days indicated a smaller deviation of the simulations from the observed values, except for ITHS. It is observed from the fitness values (SSE) in Table 6 that NSHS appears to be the best performing estimator for IR data set, while ITHS emerged as the least performing estimator. For the IR data set, it is observed that the simulations are reasonably close to the observed cumulative cases, except for ITHS. This shows that most of the simulations managed to mimic the observed trend well. Excluding ITHS, we can group the algorithms with similar simulations as (i) ESHS, NGHS, and GDHS; (ii) IHS, IBGHS, and SHS; and, (iii) SGHS, PAHS, and NSHS. Note that the difference in the simulation of ITHS was better in the previous three data sets, but it appeared to be the worst in this data set. The combination of algorithms that produced similar simulations is also different from combinations in previous data set, notably the FR data set.  Table 7 provides the optimized epidemiological parameters and the corresponding fitness values (average from 10 independent runs) for the SG data set. Figure 6 displays the visualization of the 220 days simulation, together with the 20 days projected simulation using the optimized parameters. The simulations are distinguishable, thus the visualization is informative, except for a slight overlap between ITHS and IBGHS. The simulations band begin to diverge from the observed cumulative cases approximately around the hundredth day of the outbreak, except for the simulations of ESHS, PAHS, and IHS, which diverged gradually. While other simulations deviated from the rest, ESHS, PAHS, and IHS were close to each other until day 200 of the outbreak. Only after that, the simulations diverged from each other in a small amount up to the projection period. All of the simulations were far from the observed values, resembling a similarity with the case of simulations in the SK data set. The projected simulation for subsequent 20 days indicated an even larger deviation of the simulations from the observed values. It is observed from the fitness values (SSE) presented in Table 7 that IHS appears to best performing estimator for SG data set, while ITHS is the least performing estimator, similar to IR data set.

IR Data Set (Ireland)
For the SG data set, it is observed that the simulations are quite far off the actual cumulative cases, with ITHS being the farthest and IHS being the nearest. Nevertheless, IHS's simulation is still far off the observed cumulative cases. The simulations were not able to resemble the observed trend well and it is undoubtedly not as good as for USA, FR, and IR data sets. The quality of simulations are similar to the SK data set. The simulations of ITHS and IBGHS are really far off the observed values. We note that the simulations of ITHS and IBGHS are similar, NSHS and SGHS are close to each other, GDHS and NGHS are approximately close, ESHS and PAHS are close to each other, while SHS and IHS were not similar or close to the rest.

Performance Comparison
Following the SIR simulation experiments that were performed on the five data sets, it is noted that the performance of each algorithm (based on the fitness values (SSE)) varies across the data sets. Based on the parameter estimates for each data set presented in Tables 3-7, the estimates produced by each algorithm are fairly similar, which indicates that the optimization performed is consistent due to the underlying nature of HS, regardless of the type of HS variant. Table 8 lists the best performing algorithm for each data set. Apparently, there is no one clear winner algorithm for this particular application of HS.

Data Set
Best Algorithm   USA  IHS  FR  IBGHS  SK  ITHS  IR  NSHS  SG  IHS We obtain the RMSE values from the 20 days projected SIR simulation while using Equation (24) in order to statistically compare the performance of each algorithm within each of the data sets. The Friedman test is conducted to determine whether there is any statistically significant differences among the estimation performances (in terms of RMSE) of each algorithm. Table 9 displays the RMSE values. The low RMSE values within each data set indicates that the predictive capability of SIR model while using the HS optimized parameters are fairly good. Eventually, it attests that the parameter estimation accuracy of HS variants is satisfactory. The Friedman test elicited no statistically significant difference in the estimation performance of the HS algorithms at significance level of 0.05 (χ 2 (9) = 11.749, p = 0.228). It is noteworthy to observe that the estimation performance of SHS is comparable with the rest of the HS variants and fairly consistent across the data sets. The simulations of SHS are also reasonably close to the observed cumulative cases in FR, SK, IR, SG, and USA data sets (in the closest order). The post-hoc analysis for SHS is performed using the Wilcoxon signed-rank tests in order to identify whether there is any statistically significant difference between the estimation performance of SHS and the rest of the HS variants individually. Table 10 displays the test results. The Wilcoxon signed-rank tests that were conducted with a Bonferroni correction applied at the resulting significance level of 0.05 9 = 0.006 elicited no statistically significant difference in the estimation performance of SHS when compared to the rest of the HS variants individually. Indeed, the insignificant difference supports that the performance of SHS is comparable. Therefore, although SHS may not be the best performing algorithm for any of the data sets, the consistency and statistical tests' results elucidate that SHS is competent enough to be a potential efficacious estimator for the epidemiological parameters of SIR model. A slight manual fine-tuning of the control parameters may do the job of increasing the estimation accuracy of SHS. In essence, the primary advantages of applying HS to estimate the epidemiological parameters of compartmental models are as follows:

1.
No initial values for the epidemiological parameters are required. One does not need to adapt the values of epidemiological parameters from previous studies, so as to alleviate any bias in the estimation process.

2.
No specified upper and lower bounds for the epidemiological parameters are required to suit the data sets. The burdensome process of finding an appropriate range of the parameters for each data set are evaded by using the complete range of the parameters, regardless of the data set. This may increase the computational time, but it can be traded off with better computing resources.

3.
No in-depth information about the infectious disease is necessary. HS optimization can well be applied to other infectious disease modeling without extra specific information about the disease.

Conclusions and Future Work
The application of HS is a novel approach in the field of epidemiology, particularly in epidemiological modeling. In this study, HS was implemented in order to estimate the epidemiological parameters of the prototypical compartmental SIR model as an optimization problem. Ten variants of HS algorithm were applied on five data sets to simulate the trajectory of COVID-19 cumulative infectious cases. The computational experiments demonstrated the ability of HS to be successfully applied to epidemiological modeling and as an efficacious estimator for the model parameters. As such, HS is proposed as a potential alternative estimation tool for the epidemiological parameters. An interesting insight from this study is that SHS is competent enough and it exhibited comparable performance with the rest of the HS variants in this particular application of HS optimization. For future work, the application of HS can be expanded to parameter estimations in advanced compartmental epidemiological models (e.g.,: SEIR model) and to the modeling of other existing infectious diseases (e.g.,: H1N1) or potential novel infectious diseases in the future.  Data Availability Statement: Publicly available datasets were analyzed in this study. This data can be found here: https://www.worldometers.info/coronavirus.

Conflicts of Interest:
The authors declare no conflict of interest.