Reassessment of contact restrictions and testing campaigns against COVID-19 via spatio-temporal modeling.

Since the earliest outbreak of COVID-19, the disease continues to obstruct life normalcy in many parts of the world. The present work proposes a mathematical framework to improve non-pharmaceutical interventions during the new normal before vaccination settles herd immunity. The considered approach is built from the viewpoint of decision makers in developing countries where resources to tackle the disease from both a medical and an economic perspective are scarce. Spatial auto-correlation analysis via global Moran's index and Moran's scatter is presented to help modulate decisions on hierarchical-based priority for healthcare capacity and interventions (including possible vaccination), finding a route for the corresponding deployment as well as landmarks for appropriate border controls. These clustering tools are applied to sample data from Sri Lanka to classify the 26 Regional Director of Health Services (RDHS) divisions into four clusters by introducing convenient classification criteria. A metapopulation model is then used to evaluate the intra- and inter-cluster contact restrictions as well as testing campaigns under the absence of confounding factors. Furthermore, we investigate the role of the basic reproduction number to determine the long-term trend of the regressing solution around disease-free and endemic equilibria. This includes an analytical bifurcation study around the basic reproduction number using Brouwer Degree Theory and asymptotic expansions as well as related numerical investigations based on path-following techniques. We also introduce the notion of average policy effect to assess the effectivity of contact restrictions and testing campaigns based on the proposed model's transient behavior within a fixed time window of interest.

Introduction

COVID-19 outbreaks have been curtailing socio-economic activities around the globe. Over 150 million total confirmed cases had been reported by Apr 29, 2021, and the number of deaths exceeded 3 million by Apr 16, 2021, reflecting the burden of the pandemic [1]. This unprecedented health crisis has shown how far time and spatial propagation of incidence matter to each individual on a micro-scale and subsequently to a country on a macro-scale. Toward the ultimate herd immunity, several vaccines have been introduced; however, their efficacy must be scrutinized amidst virus mutations [2, 3]. World Health Organization sets a minimum efficacy of 50% with a preferable threshold of 70% [4]. Although many of the vaccines are well above these efficacy levels, effectiveness in the field might be different due to the variations in affordability, public compliance, healthcare planning, etc. [5, 6]. Moreover, equitable access to vaccines, in particular for developing countries, is also a challenging task [7]. Therefore, all the non-pharmaceutical interventions (NPIs) by means of contact restrictions (physical distancing, wearing face masks, washing hands, crowd clearance, workplace clearance, school closure, lockdown, public curfew, mobility restriction), and testing campaigns (including contact tracing) must be maintained until vaccination programs take substantial control over the further spread [8-10]. Many developing countries are still subject to financial restrictions against the import of vaccines [7], and at the same time, NPIs give a variable impact due to wavering laws and public compliance that mostly weigh upon socio-economic reasons [11]. As far as the spatial aspect is concerned, these NPIs should be implemented considering disease and societal impact according to international, national, and regional epidemiological situations [12]. Research on the actual performance of NPIs in developing countries is limited, and thus related government decisions usually are over- or underestimated [13, 14]. It further creates a dilemma on what is more important between intra-regional and inter-regional contact restrictions, in particular for reopening the economy [15].

As vaccines with yet unknown success rates toward herd immunity are not even equally affordable across different economic classes, the only alternatives are enforcing laws and reshaping public awareness toward upholding NPIs. In relation to the spatial aspect, we start our investigation with the following questions: This work puts up not only the prioritization of healthcare capacity and NPIs among spatial units but also the route for them in a more robust way than incidence-driven approaches. Endeavor to this has been known from the field of spatial mapping, namely to group spatial units into meaningful clusters.

On what sense may the decision maker appropriately perform the prioritization of healthcare capacity (e.g., hospital beds, ICU units, testing capacity, monitoring quarantine, including limited vaccines) among all spatial units in a country?

Under limited data of confounding factors, how can the decision maker value and reassess the flow of epidemics as well as the impeding NPIs?

In Sec. 2, we adopt global Moran’s index and Moran’s scatter to measure the timely spatial pattern of COVID-19 incidence in a country as well as to set the grouping. Particularly in developing countries, prioritizing high-risk areas or hotspots is driven by careful utilization of healthcare capacity [16]. The two aforementioned tools stand out among simplistic case mappings for their power to localize and group hotspots. Accordingly, priority for intra-cluster NPIs remains the same within a cluster but sequential between clusters. This strategy is important for developing countries like Sri Lanka that has not yet been covered by a holistic spatial analysis of this caliber. In addition to prioritization and route, the clustering study can bear the locations for placing border controls, which in this case are those in the main inter-cluster mobility streams. There remains, however, one caveat from these tools. That is, they are not able to parameterize the ongoing government decisions in terms of numbers and thus fail to impart how sensitive the incidence is against changes in those decisions.

Focusing on Sri Lanka, in Sec. 3 we propose a metapopulation model for Moran’s clusters determined from available panel COVID-19 incidence data. The preference of the dynamic model over functional regression models stems from integrable mechanistic processes behind COVID-19 infection and that no spatio-temporal data of confounding factors were found. A complexity reduction is proposed based on the unavailability of related field data, resulting in a simple model but rational enough such that contact restrictions and testing campaigns are mediated. Sects. 3.2–3.4 are then devoted to studying the likelihood if the incidence persists for a long time. To this, the model solution is compared with certain equilibria in the local sense whereby the basic reproduction number and effective reproduction numbers play the key role within. The model fitting as in Sec. 4 will provide a proxy to not only the approximate reproduction numbers but also non-observable dynamics, including contact matrix and the ongoing government decision on testing campaigns.

Finally, Sec. 5 extends the bifurcation analysis numerically using a path-following technique for the case where according to the fitting, the clusters are not strongly connected. In addition to this, the performance of the government decisions on contact restrictions and testing campaigns during the observations is reassessed via maximal average policy effect, which measures the average number of individuals per 1,000,000 inhabitants that could have been saved from COVID-19 infection on the virtue of better interventions. Scenarios to cost-to-benefit ratio are also presented alongside.

Spatio-temporal analysis

This section is devoted to answering question 1 in the context of Sri Lanka. Particularly under consideration is prioritization of healthcare capacity and NPIs as well as classification of the 26 Regional Director of Health Services (RDHS) divisions into Moran’s clusters.

Study area and observation period

Sri Lanka is a South Asian island country situated in the Indian Ocean between latitudes 555’ and 950’ N and between longitudes 7931’ and 8153’ E. Sri Lanka has a population of about 21.9 million [17]. From an administrative perspective, the country is divided into 9 provinces that cater to 25 districts. In health administration, there are 26 Regional Director of Health Services (RDHS) divisions that mainly coincide with administrative districts, except the district Ampara that is covered by two RDHS divisions. The primary units of health administration are called Medical Officer of Health (MOH) areas. There are 356 MOH areas wherein the health surveillance activities are carried out [18]. Over 100,000 total confirmed cases and 600 deaths had been reported in Sri Lanka by Apr 24 and by Apr 13, 2021, respectively [19]. The public has been asked to follow health guidelines such as wearing masks, washing hands, and keeping one-meter distance since the early stage of the outbreak [20]. All the confirmed cases are directed to hospitals, and close contacts in addition to overseas returnees are requested to be quarantined [20]. The data used in this study are the daily new cases recorded by the Epidemiology Unit, Ministry of Health of Sri Lanka, spanning over the period from Nov 14 until Mar 31, 2021 [21]. Recording data in RDHS level began from Nov 14, which lies within the post-curfew period after major superspreading events (apparel factory cluster [22] and fish market cluster [23]). Earlier to that, the data had been listed only according to the clusters arisen from superspreading events and quarantine centers. This is due to that only several clusters were significant rather than a community level spread up to the end of Oct 2020 [21]. The RDHS-wise normalized daily new cases (per 1,000,000 inhabitants) are illustrated in Fig. 1. Note that no major mobility restrictions had been imposed within the observation period.

Fig. 1

Daily COVID-19 new cases in RDHS divisions per 1,000,000 inhabitants (a–b) and their timely average including spatial average and maximum (c). The RDHS divisions are COL (Colombo), GAM (Gampaha), PUT (Puttalam), KAL (Kalutara), ANU (Anuradhapura), KAN (Kandy), KUR (Kurunegala), POL (Polonnaruwa), JAF (Jaffna), RAT (Ratnapura), KEG (Kegalle), MON (Moneragala), KLM (Kalmunai), MAL (Matale), GAL (Galle), AMP (Ampara), BAD (Badulla), MAT (Matara), BAT (Batticaloa), HAM (Hambantota), VAV (Vavuniya), TRI (Trincomalee), NE (Nuwaraeliya), KIL (Kilinochchi), MAN (Mannar), and MUL (Mullaitivu)

Global Moran’s index and Moran’s scatter

For spatial auto-correlation, interconnectivity between spatial units indexed by i and j is usually represented by a spatial weight matrix . These weights can be designed according to shared boundaries of spatial units or distance between centers. The usual adjacency matrix can be an option, but a distance measure may better articulate connectivity since adjacency only captures interaction among neighbors. In our case, the distances among RDHS divisions are based on placing appropriate centers , which are taken from averaging those from MOH areas , namely , weighted by their population . The centers consist of the latitude and longitude of the most attractive points, for example a city center, main administrative/commercial building, transport hub, main junction, etc. It then followsNow that the distances are computable by the standard Haversine formula, we take the power functional form [24] of the weightThe exponential decay parameter serves to scale the influence of the distance while the threshold distance cuts the inessential interconnectivity. It is important to note that sufficiently large d values help make W irreducible, i.e., all regions become strongly connected.

Suppose that time is frozen and the mean normalized cases over the period shown in Fig. 1 for RDHS divisions are reported as with mean . Taking , the global Moran’s index [25] with a row stochastic matrix W as in (2) is given byThe global Moran’s index basically is the Rayleigh quotient of evaluated at Z, which brings the spatial autocovariance standardized by the variance of the data. The interpretation of the index usually comes in connection with the so-called Moran’s scatter where . It is quite apparent that the latter compares every spatial unit’s self-incidence magnitude against the mean with the weighted magnitudes from its corresponding neighbors as spatial lags of the unit. The four Moran’s clusters are then the cluster Q1 (first quadrant in 2-dimensional Euclidean space) referring to a set of spatial units of high incidence surrounded by their spatial lags of high incidence (high–high, hotspots), the cluster Q2 (second quadrant) for spatial units of low incidence surrounded by their spatial lags of high incidence (low–high), the cluster Q3 (third quadrant) for those of low incidence surrounded by their spatial lags of low incidence (low–low, coldspots), and the cluster Q4 (fourth quadrant) for those of high incidence surrounded by their spatial lags of low incidence (high–low). We obtain two facts accordingly. First, the regressing line of that passes through the origin has the slope . Second, if and denote the minimum and maximum eigenvalue of the symmetric matrix , then the standard Rayleigh–Ritz (min–max) theorem (see e.g., [26] or [27]) suggests that . This gives somewhat the tightest range due to . Since the diagonal entries of W are 0, evaluating the Rayleigh quotient at any of vectors in the standard basis of , namely , yields . If , then more points are aligned with the regressing line of that slope, making Q1 and Q3 full of points leaving out Q2 and Q4 scarce. A locally clustered spatial pattern is then observed. If and in case , then points are more concentrated in Q2 and Q4, indicating a locally dispersed spatial pattern. In between, under , there is no relation between self-incidence magnitudes and those from their spatial lags, leading to a random spatial pattern. We shall comment that the case may be observed in many cases where but generally not always true. Besides assuring the upper bound 1 to any of the aforementioned bounds of , sufficient conditions for this may include: W is symmetric (doubly stochastic) such that , W and commute in which case [28], and W is diagonalizable since then .

As far as Sri Lankan data are concerned, a technical question arises: which values of and d in the weight matrix are suitable for the data? We answer this question by computing the smallest absolute elasticity indices of Moran’s index on the average new cases. Now suppose that is decreased to a certain percentage from its current value, i.e., , where . In this way, represents the total percentage post perturbation and the percentage of increment. Taking this definition of is more technically sound for a comparison among parameters as they may live in disparate scales. In response, also changes from its initial data in the same fashionFor “fair” treatment, one usually designates , which is sufficiently small. Therefore, the first-order terms from the previous expressions take the role in determining which parameter, to which is more sensitive. We then say is more sensitive to the increase of than d in the regimeIn the literature, e.g., [29, 30], these two compared expressions in (4) are called the (first-order) elasticity indices. There is one issue, namely the non-smoothness of the index with respect to d limits the definition to its approximation; see Fig. 2a–d. Apparently, Moran’s index is highly sensitive to d in case is relatively small () but insensitive to d as . By m and , the elasticity indices are roughly zero, meaning that Moran’s index changes only very slightly under the variation of in a neighborhood of these values. Additionally, plotting Moran’s scatter on a daily basis (Fig. 2f) gives maximal concurrence percentages across RDHS divisions that agree with Moran’s scatter on the average data (Fig. 2e). To the latter, we obtain (p value ) indicating a locally clustered spatial pattern for the average data.

Fig. 2

Approximates of the elasticity indices using (a) and (b) and using (c) and (d) . Moran’s scatter plot using m and for the average daily new cases is presented in (e) while (f) gives the timely Moran’s scatter

Accordingly, we classify the RDHS divisions as follows: cluster Q1 (COL, GAM, and KAL); cluster Q2 (KEG, MAL, GAL, and NE); cluster Q3 (PUT, KUR, ANU, POL, JAF, MON, AMP, BAD, MAT, BAT, HAM, VAV, TRI, KIL, and MUL); cluster Q4 (KAN, RAT, KLM, and MAN). From the application point of view, the cluster Q1 amasses all the hotspots. Ameliorating the burdens of infection follows from putting a first-level priority on healthcare capacity and possible vaccinations as well as providing more strict border controls that would reduce mobility from and to its spatial lags, i.e., neighbors in the sense of the weight matrix. Intra-cluster border controls cannot change the situation much, but interventions can be realized through the applications of NPIs including public curfew and testing campaigns. We argue that the intra-cluster prioritization as well as deployment route for NPIs can be left to the decision maker, which can rely on the available resources. The cluster Q4 requires not only a second-level priority on healthcare capacity but also mobility restrictions to its spatial lags, otherwise the epidemics outwardly diffuses. Meanwhile, the cluster Q2 may receive a third-level priority as well as isolation from its spatial lags such that it does not attract epidemics. Lives in the cluster Q3 can be the easiest in terms of mobility restrictions as long as reasonable hygiene practices and physical distancing are upheld. Border controls can now be localized to any point that shares the borders between clusters, which could be an intersection point on the main road or an intersecting railway station. For example, no border controls are required between KUR and ANU, but between KUR and COL. If the authorities follow a clustering based on administrative provinces, the border between KUR and ANU should be controlled as they belong to separated provinces. Thus, our analysis suggests more scientific clustering that may overrule general administrative choices.

Modeling-based reassessment of NPIs

Now that the 26 RDHS divisions are classified into the four Moran’s clusters, this section is devoted to modeling the incidence dynamics on the clusters and parameterizing ongoing government decisions on contact restrictions and testing campaigns toward answering 1. Here is the idea: Once the essential performance measures for contact restrictions and testing campaigns are gained through the model fitting, we can optimize the model toward specific goals whereby different magnitudes of the measures are tested. We focus on three goals, namely minimizing the basic reproduction number, maximizing the average policy effect, and minimizing the associated policy cost. All forms of NPIs can be reassessed toward these goals.

The nature of standard metapopulation models suggests that, in contrast to the kinematic models, the whereabouts of every single individual are no longer concerned. Among first generations of metapopulation model, two-patch models were proposed for their accessibility to sophisticated analytical investigation on the disease endemicity via the basic reproduction number [31-34]. These studies shared similar results: the disease-free equilibrium is globally asymptotically stable if , and all the state variables are uniformly persistent if leading to the existence of an endemic equilibrium, which was proven to be globally asymptotically stable. The SIS model in [34] stands out among the cited models as it incorporates infections during travels. General n-patch SIR-type models considering mobile humans with memory over their origin zones admit short visits to other zones that allow them to infect other humans or to acquire infections ex-situ [35-37]. The notion of transit time becomes the key determinant to the latter. Models without memory were proposed in [38, 39] where in [38], a more generalized population growth function was used, taking into account the relationship between and the disease extinction and persistence. Metapopulation models for COVID-19 have also appeared recently. Citron et. al. [40] consider metapopulation versions of an SIR, an SIS, and a Ross-Macdonald model integrating Eulerian movement (direct out-flux) and Lagrangian movement model (net out-flux and influx). The two movement models were analyzed to synthesize conditions under which one model can be superior against the other with respect to epidemiological outcomes. A model including transit time and infection due to exposed cases was proposed in [41] With known network attributes of the tested case, the study was brought to determine the infection rates and the ratio between asymptomatic and symptomatic cases. Recently, metapopulation models including vaccinated compartments [42] and age structure [43] were proposed and validated using field data. A memory-less migration (or diffusion) model including human mobility [44] was used for modeling daily confirmed cases on a network of 343 cities in China.

In this study, our model is concerned with the COVID-19 epidemics that naturally include undetected and deceased cases. We use the concept memory in the model, but unlike in [35-37], the number of humans from cluster i that are in cluster j and thus at which cluster the contacts happen, are not displayed. In contrast to [42], we combine the infection rate and the matrix representing the fraction of total daily time for i-residents to be in j-region into what we called a contact matrix.

Metapopulation model for COVID-19 epidemics and biological assumptions

We divide the regional population () into five subpopulations based on their health status: susceptible (healthy but vulnerable to infection), detected active cases comprising some portions of asymptomatic and symptomatic (hospitalized) cases, undetected cases (dark figures, mostly asymptomatic), recovered , and deceased cases . Due to a small incubation duration, we count the intermediate exposed (pre-symptomatic) cases to the susceptible subpopulation to simplify the model presentation. Net population growth due to imports, migrations, natural births, and deaths is assumed to be negligible during the observation period, inducing constant total cluster population . The point of departure from our modeling is concerned with

In this basic model, denotes the infection rate that determines the likelihood of a susceptible person from ith region to meet with an infected person from jth region. In the standard SIR models for airborne diseases, the infection rates depend on many factors including sneezing rate, probability of sneezing during encounters [45], infectiousness measure (viral load, case index) [46, 47], effectiveness measure determining how probable an average susceptible person contracts infection (health condition, age) [45], human mobility for bearing corrections of the possible number of encounters [45], influence of media reports on public awareness [48], and possibly weather factors that enhance aerosol transmissions [49-51]. As the exposed cases were gathered in , the infection rates also give another correction as the individuals cannot both be infected and spread the viruses. After contracting infection, the remaining time known as viral shedding period determines the average duration from the onset of symptoms until the cessation of viral shedding (when viruses can no longer be released from an infected person), indicating the end of the infectiousness period [52, 53]. The parameters and m denote the duration of temporary immunity and the fatality rate from the detected cases. We impose a strong assumption that during the limited observations, the entire infected cases are timely distributed into the detected and undetected cases with the average proportions and , respectively. To accommodate different transmission scales from detected and undetected cases, we have used the parameter . Finally, denotes the natural birth or death rate.

Our next task is to simplify the model even further. Due to unknown dark figures , several ideas and estimates have been appearing in the literature, see e.g., [54]. Ours is based on the assumption that the dark figures proportionate the detected cases to a certain constant, i.e., where for all clusters i and time. By the range of p, we impose that most cases are detected. As we specify , Eqs. (5b) and (5c) apparently become equivalent. This choice justifies the idea that the constant ratio between undetected and detected cases requires constant detection rate (in the sense of averaging) and that the detection rate also holds . Apart from this, we bring forward the non-observability assumption to due to data credibility. Our study designates as to proportionate to a certain constant from time to time, namely for a new parameter and all i. It is straightforward to see that . In the next model presentation, we would like to use the re-scaled variables and with as well as the following notations , , , , as the contact matrix, as the diagonal matrix whose main diagonal is S, and as a matrix or a vector whose entries are 1. We acquire the final model with an initial value . This model portrays the situation where all infected cases are distributed to the detected classes in case , i.e., when the quality of the testing campaigns is extremely good. Moreover, as much as half of the infected cases will be distributed to the detected cases when no essential tests are done, i.e., when . In Sri Lanka, PCR and antigen tests are carried out on a random and targeted basis [55]. However, limited financial allocations may curtail arbitrary increase in testing capacity [56]. Another factor for a large p is the compromised public compliance to tracing technology that tolerates the effectiveness [57]. As a result, lack of tests retards the process of unraveling possibly infected close contacts and thus hotspot identification, which eventually delays blocking the routes of transmission [58].

Basic reproduction number

We study the basic reproduction number to determine the local behavior of model solution around the disease-free (DFE) and endemic equilibrium (EE) for fleeting observations. Therefore, given that optimal parameters are subject to data availability, the predictive power of this behavioral analysis is limited to a short-range prediction window. Let be a point of interest that is compared to the solution of the model represented by (6). For simplicity, we assume that is either DFE or EE such that . The error measure then follows providing that z is close enough to 0 or x is close to . A compelling property of such a linearized system is that the short-term trend of the model solution around 0 can be predicted by the local (even global) stability. The basic reproduction number will then be used to parameterize a condition for the maximal real part of the eigenvalues of Jacobian matrix , which eventually determines the local stability of z around 0.

When , we obtainwhere denotes the identity matrix. The next generation matrix as well as the basic reproduction number can now be formulated asrespectively. Here, denotes the spectral radius of G. According to Berman and Plemmons [59], becomes a nonsingular M–matrix if and only if or . The fact that being an eigenvalue of G is equivalent to being the corresponding eigenvalue of (with the non-changing eigenvectors), Perron–Frobenius Theorem on simplicity and dominance of also guarantees that holds true if and only if all other eigenvalues of (or ) lie in the open disk of center and radius (or center and radius ). Additionally, we acquire another fact that (or ) becomes singular if and only if , in which case a zero eigenvalue occurs. Due to the equivalence relations, the final case happens if and only if there exists at least one eigenvalue of with a positive real part. We obtain the following summary: z is attracted to 0 or DFE becomes locally attractive to the solution x of (6) if and it becomes locally repelling to x if .

Existence and attractiveness of an endemic equilibrium

Computing an endemic equilibrium (EE) from model (6) also returns complexity on its own. The second subsystem (6b) gives the equilibrium equationfor all j. Substituting this to the first subsystem (6a) together with also multiplying every i-th entry with gives usfor all i. The preceding system of equations folds under multiplication by intowhere G denotes the next generation matrix as in (8). This is a multidimensional quadratic equation whose solutions cannot be derived explicitly. In order to guarantee the existence of EE, we first see if the point is indeed a branch point of the quadratic equation. Since is a quadratic function defined on some open subset , it is verifiable that . Let be a point such that , where denotes the boundary of . The point q is said to be regular if either for all points return invertible . Otherwise, q is called critical.

Adopting definitions from [60, 61], the map defined aswith being regular such that , denotes the Brouwer degree of in with respect to a reference point q. Another convention narrows the singular value down to with the neighborhood of 0 is chosen in such a way that is isolated. In this case, the mapdefines the index of at the isolated singular value . According to the last two references, is a branching point providing that changes values around .

In case , the fact that the multiplication between complex conjugate numbers return a positive number, gives us . We can always impose continuous perturbation on parameters such that s(0) solves and is equivalent to for and some . In case , the eigenvalue of with the largest real part apparently returns and the other eigenvalues lie in the open disk of center and radius 1. We can appoint the eigenvalue of with the largest negative real part and of algebraic multiplicity , and define . The function with G(0) corresponding to is holomorphic in and continuous in . We can appoint such that is the only root in the closed disk . There must now exist such thatholds for all and . According to Rouché’s Theorem [62], has roots in of counting multiplicities when . The same continuity argument can be used to show that all the remaining eigenvalues can never have largest negative real part which exceeds . In summary, as for a new it holds that 1 is not an eigenvalue of G and slightly increases from 1 such thatfor all eigenvalues of G. This returns two results: (i) becomes non-singular such that serves as an isolated singular value of in its neighborhood due to there and (ii) . The index of at the singular value now reads asfor . This confirms that that is indeed a branching point.

The next task is to verify the positivity of the local branch. We took the asymptotic expansion of from 1 [63, 64], i.e., the coefficient of in the quadratic equation (11) via the direct relation between and :such that the branch I takes the expansionSubstituting the preceding expressions to the quadratic equation (11) returnsZeroing the first-order term () gives usThis means that is the eigenvector of G associated with , whose existence and positivity are guaranteed by Perron–Frobenius Theorem. The latter also guarantees the existence and positivity of the left eigenvector associated with . Now, multiplying the second-order term () with from the left gives uswherebyby . Moreover, substituting from (6b) to (6a) leads us to the following summaryThese parametric expressions suggest that as increases from 0, increases from 1 and a unique local branch I takes off from 0 with the initial direction with respect to whereby the susceptible state decreases from c simultaneously for all clusters. Finally, one yieldsThis indicates the existence of a continuum of endemic equilibria in the neighborhood of and in the direction of increasing . For , let us write one endemic equilibrium EE as with the expression given in (20). The Jacobian matrix evaluated at EE takes the formThe matrix in the leading order has eigenvalues of algebraic multiplicity 4 and where are the eigenvalues of G. Due to simplicity and dominance of , all the eigenvalues locate in the open disk of center and radius where only is in the origin. We can use Rouché’s Theorem one more time with the function to show that all eigenvalues of , except the one that corresponds to , stay in the open left-half plane in for a sufficiently small .

The fate of this last eigenvalue can be analyzed as follows. The eigenvalue of associates with the right and left eigenvector (by and of ):respectively. Using Taylor expansion on a simple eigenvalue of a perturbed matrix [65], we obtain the eigenvalue of the Jacobian matrix that corresponds to :by substituting from (20) and taking Taylor expansion over . This shows the existence of where implies all eigenvalues of having negative real part. Combining with (20) and (21), we acquire a forward bifurcation of the model system (6) at . This means that the local branch of EEs becomes locally attractive as .

Effective reproduction numbers

Providing epidemics are going on, , we have from (26):Observe that if and only if the local instantaneous reproduction numbersfor all i. Practically speaking, if the active cases I determines the endemicity levels, speaks about the epidemics progression. The fact that gives the inflow of new cases in all clusters, gives the expected number of new cases from at the timestamp t attributed to the entire infected individuals from () per viral shedding period . Therefore, represents the expected number of new cases from attributed to the normalized infected individuals () per viral shedding period . For a realistic approximation, we took the inflow of new cases from the data while the active cases, which serve as the denominators, will be taken from the fitting. Such a definition of instantaneous reproduction number has been used by Fraser in [66], except where the active cases (as the denominator) were taken from weighted new cases in the past n days for a fixed n. The weights were later known as serial intervals [67, 68], estimating the distribution of delays taken from the onset of symptoms until hospital admission (i.e., when the data of new cases are usually recorded). Under the two facts that (1) the instantaneous reproduction number is, by the definition, too fluctuating and (2) infected individuals can already infect susceptible individuals from the onset of symptoms, Fraser also introduced some moving average such that the ‘real’ new cases at a certain timestamp should come from the new cases ‘recorded by hospitals’ in the future timestamps (up to n) weighted by the serial intervals, while the active cases come from the sum of those corresponding to the used timestamps, where again, at each timestamp the active cases are weighted sum of new cases in the past n days. Inspired by such refinement with, however, lack of serial interval data, ours becomesfor some averaging window size . The forward moving average thus allows the serial intervals to be of uniform distribution around days. In the numerical computations, we designate .

Data assimilation

The basic aim of parameter estimation is to find agreement between model solution for weekly new cases at all time points with known data subject to identifiability of unknown parameters . We assume that the fitting would be subject to time-invariant i.i.d. error of the weekly covariance (for all the four clusters) and the prior was set to be Gaussian. The latter means that the parameter estimation will be based on minimizing the Mahalanobis distance between the model solution and empirical data. For simplicity, no correlation was imposed for the cluster-wise error such that . The non-degenerate joint likelihood function for one time point is then given byAssuming timely i.i.d. measurement, the joint likelihood for the entire observations is then given byby taking that serves to simplify the representation of the likelihood function [72]. Our study designates the variance terms as the mean of the data throughout the observations so as to avoid a blow-up in the likelihood function.

As the parameter dimension is much smaller than the data size, the standard asymptotic confidence interval [73] has been suggested to delineate the parameter uncertainty [74, 75]. The formula of the confidence interval for each optimal parameter takes the formThe operator denotes the inverse of the Hessian, while the notation denotes the quantile of the distribution with the degree of freedom df. The degree of freedom can be chosen between two that further determines the type of confidence interval: gives pointwise asymptotic confidence interval (PACI) that works on the individual parameter, gives a simultaneous asymptotic confidence interval (SACI) that works jointly for all the parameters.

In the present study, the Hessian matrix in (25) will be approximated up to the second order using the queen-type stencil. Due to disparate scales of the parameters, the step size will be made dependent on the parameter’s order of magnitude, i.e., for a uniformly small . Our study uses . After all, the fitting will be done in MATLAB using the toolbox fmincon accompanied by interior-point as the core optimization solver. The fitting result together with the effective local reproduction numbers can be seen in Fig. 3. Meanwhile, we keep at the fitted values, we vary from their confidence interval to have a shaded region around the fitting curves. Due to model simplicity (no time-dependent parameters), we can only expect to see an almost stationary model solution to fit the almost variance-stationary dataset, also subject to the constraint on of the four clusters: . The fitted parameter values can be seen in Table 1.

Fig. 3

Fitting results from the deterministic SI model. The vertical axes give the numbers of weekly new cases per 1,000,000 inhabitants and effective local reproduction numbers from all clusters. The fluctuating curves represent the data and the shaded region around the fitted model solution determine the variation of from their confidence interval

Table 1

Parameters of the SI model (26). All zero -values were due to rounding numbers smaller than 1e–07. This is intentional against floating-point error in the numerical continuation, while at the same time, almost no visible difference in the model response in comparison to that using positive values was observed

Parameter	Description	Unit	Range	Value (ε PACI)	Ref.
Known:
1/μ~	Average human lifespan	[d]	70-80·365	76.9·365	[69]
1/γ~	Viral shedding period	[d]	1-30	20	[70]
γ	μ~+γ~	[d^-1]
c	Scaling factor	–		10^6
Optimized:
ϱ	Transmission scale from undetected to susceptible	–	1-20	11.6373 (0.6727)
p	Case detection ratio	–	0-35/65	0.4698 (0.0148)	[71]
η	Rescaled loss-of-immunity rate	[d^-1]	1/3·365	4.5666e–04 (2.0482e–05)
β_11	Infection rate between S_1 and I_1	[d^-1]	0-3/c	0.0040 (2.4694e–05)
β_12	Infection rate between S_1 and I_2	[d^-1]	0-3/c	0.0308 (5.7954e–04)
β_13	Infection rate between S_1 and I_3	[d^-1]	0-3/c	0 (0)
β_14	Infection rate between S_1 and I_4	[d^-1]	0-3/c	0 (0)
β_21	Infection rate between S_2 and I_1	[d^-1]	0-3/c	0.0032 (1.3493e–04)
β_22	Infection rate between S_2 and I_2	[d^-1]	0-3/c	0 (0)
β_23	Infection rate between S_2 and I_3	[d^-1]	0-3/c	0 (0)
β_24	infection rate between S_2 and I_4	[d^-1]	0-3/c	0 (0)
β_31	Infection rate between S_3 and I_1	[d^-1]	0-3/c	0.0012 (4.3181e–05)
β_32	Infection rate between S_3 and I_2	[d^-1]	0-3/c	0 (0)
β_33	Infection rate between S_3 and I_3	[d^-1]	0-3/c	0.0126 (3.0470e–04)
β_34	Infection rate between S_3 and I_4	[d^-1]	0-3/c	0 (0)
β_41	Infection rate between S_4 and I_1	[d^-1]	0-3/c	0.0046 (5.8647e–05)
β_42	Infection rate between S_4 and I_2	[d^-1]	0-3/c	0 (0)
β_43	Infection rate between S_4 and I_3	[d^-1]	0-3/c	0.0039 (4.4670e–05)
β_44	Infection rate between S_4 and I_4	[d^-1]	0-3/c	0 (0)
S_10	Initial condition for S_1	–	0.01c-c	8.9227e+05 (2.3154e+04)
S_20	Initial condition for S_2	–	0.01c-c	9.3000e+05 (4.7650e+04)
S_30	Initial condition for S_3	–	0.01c-c	6.5110e+05 (1.6834e+04)
S_40	Initial condition for S_4	–	0.01c-c	8.1692e+05 (3.7530e+04)
I_10	Initial condition for I_1	–	0-0.1c	569.9877 (18.4876)
I_20	Initial condition for I_2	–	0-0.1c	569.9791 (10.1970)
I_30	Initial condition for I_3	–	0-0.1c	37.4413 (0.8390)
I_40	Initial condition for I_4	–	0-0.1c	569.9829 (17.9763)
Adjusted:
μ	μ~(1+p)	[d^-1]		5.2365e–05
T_Ref	Reference time for transient analysis	[d]		136
p_Ref	Reference value for average policy effect	–		p
ω_Ref	Reference value for average policy effect	–		1

Numerical study of the COVID-19 model via path-following techniques

In this section, our main goal is to investigate the dynamical response of the model as certain selected parameters are varied. To evaluate the impact of reassessment on government policy against COVID-19 posterior to fitting, we shall introduce one more control parameter that hereafter is referred to as the contact restriction factor. This parameter will serve to decrease the intra- and inter-cluster contacts as so far portrayed by the fitted values of . From the application point of view, this parameter can be realized by enhancing NPIs and all possible interventions that likely reduce the contact between susceptible and infected persons. Furthermore, the parameter p (case detection ratio) will be interpreted as a factor determining the quality of COVID-19 testing campaigns in such a way that p close to zero represents an effective testing policy, while a large p indicates that a great number of infections are not detected and therefore are able to spread the disease at higher infections rates (according to the factor in the SI model (26)). Consequently, the reassessment yields a small modification in the model asThe numerical investigation will be based on the parameter fitting obtained in the previous section. There, the pair represents the susceptible and infected population in the cluster i. In this way, our study will focus on the effect of the main disease control parameters on the model behavior including the basic reproduction numberin such a way that a fixed combination of will be interpreted as a specific disease control policy determined by the decision makers. The numerical study will be carried out using the path-following software COCO (Computational Continuation Core [76]). This is an analysis and development platform for the numerical treatment of continuation problems using MATLAB. A remarkable feature of COCO is its set of toolboxes that covers, to a large extent, the functionality of available continuation packages, such as AUTO [77] and MATCONT [78]. In particular, we will make extensive use of the COCO-toolbox ep, which encompasses a set of numerical routines for the bifurcation analysis of parameter-dependent families of equilibria in smooth dynamical systems.

Monitor and cost functions

In this investigation, one of the main goals is to study the effectiveness of the disease control policy to reduce the number of COVID-19 cases in the proposed biological scenario, and for this purpose suitable performance measures will be considered in our numerical implementation. Let us assume thatis a bounded reference solution of system (26) computed for the parameter values and initial conditions given in Table 1, with being a reference final time and . In this setting, we define the performance measure given bywhere and is the p-value given in Table 1. In the above expression, S(t), I(t), , stand for a solution of system (26) computed for the parameter values and initial conditions given in Table 1, but for arbitrary . From a practical point of view, the quantity (hereafter referred to as the average policy effect) represents the average COVID-19 cases that could have been free from infection on a daily basis by applying a particular disease control policy , in comparison to the reference solution case defined above. In connection to this definition, we introduce the associated policy cost given bywhere is a coefficient that characterizes the cost distribution among contact restrictions and testing campaigns. As can be seen from (29), a strict mobility reduction () implies a high policy cost, representing the bad impact on the economy and other negative effects associated with the mobility reduction. Similarly, a widely spread and effective COVID-19 testing campaign () also produces very high costs, due to the personals required for implementation, expenditure on test kits and other logistics, media advertisement, organization, etc. In our investigation, the value will be assigned, which portrays a realistic distribution between the two cost terms in (29) according to our numerical simulations. Nevertheless, we give such a higher contribution from contact restrictions based on bad economic impact in Sri Lanka due to job and earning losses associated with mobility restriction and crowd clearance, which additionally force the government to spend much on welfare activities targeting low-income citizens [79]. Therefore, the cost function given in (29) takes not only the view of government spending but also the economic recession in the whole country into account.

Numerical investigation of the modified COVID-19 model

With the mathematical framework introduced in the earlier section, we can now move on to the numerical study of the modified COVID-19 model (26) using parameter values and initial conditions given in Table 1. Observe that the contact matrix is no longer irreducible. As a result, the initial direction of the continuum of endemic equilibria as in (20) is only guaranteed to be nonnegative according to Perron–Frobenius Theorem. A preliminary system response can be seen in Fig. 4. The picture shows time series for the active cases and susceptible population , corresponding to Moran’s clusters Qi where . As can be seen in the figure, for the selected parameter values the system shows a damped oscillatory behavior that settles down after a long time to an endemic equilibrium, i.e., a steady state where the COVID-19 infection is present in all clusters. This equilibrium state will then be used as starting point for our numerical investigation based on path-following techniques.

Fig. 4

Dynamical response of the COVID-19 model (26), computed for the parameter values and initial conditions given in Table 1. The picture shows time series for the infected () and susceptible population ()

Fig. 5

One-parameter continuation of equilibria of system (26) with respect to the mobility restriction factor and the case detection ratio p, computed for the parameter values given in Table 1. Panels a and b depict the behavior of the basic reproduction number given by formula (8). Panels c and d present the behavior of (left vertical axis, in blue) and (right vertical axis, in red). Similarly, panels e and f plot the variation of (left vertical axis, in blue) and (right vertical axis, in red) with respect to the corresponding parameters. In these diagrams, solid and dashed lines stand for branches of stable and unstable equilibria, respectively. During the computations, a series of branching points are detected for (BP1), (BP2) (depicted in panels (a), (c) and (e)) and (BP3), (BP4) (depicted in panels (b), (d) and (f))

Fig. 6

a Two-parameter continuation of the branching point BP1 found in Fig. 5a with respect to p and . The resulting curve divides the parameter space into two regions: one for stable disease-free equilibria (yellow) and one corresponding to stable endemic equilibria (blue). b System responses obtained at the test points P1 (, ), P2 (, ), P3 (, ) and P4 (, ). The time plots present the behavior of (red), (black), (blue) and (green). All numerical simulations are calculated with the initial conditions specified in Table 1

Let us begin our study with the numerical continuation of the endemic equilibrium found above with respect to the mobility restriction factor . The result of this process can be observed in Fig. 5, panels (c) and (e). Specifically, panels (c) and (e) present the behavior of (left vertical axis, in blue), (right vertical axis, in red) and (left vertical axis, in blue), (right vertical axis, in red), respectively, as the parameter varies. Panel (a) shows the dependency on of the basic reproduction number , given by formula (27). In this diagram, it can be seen that for low values of , the basic reproduction number is smaller than one, due to which the system presents a stable disease-free equilibrium corresponding to the solid horizontal branches shown in Fig. 5c and (e). As increases, increases as well, and it crosses 1 from below at , where a branching point BP1 is detected. Here, the disease-free equilibrium loses stability and an endemic branch is born (via a forward bifurcation). Interestingly, at this point a COVID-19 outbreak occurs only for clusters Q3 and Q4, while clusters Q1 and Q2 remain disease-free. If increases further, however, the disease for clusters Q1 and Q2 develops for , where a branching point BP2 is found. From this point onward, the disease is present in all clusters, and the increment of the infected cases augments more rapidly as the mobility restriction factor grows.

A similar scenario is encountered when the case detection ratio p is considered as the bifurcation parameter, see Fig. 5b, d and f. A first branching point (from below) is found for (BP3), where a COVID-19 outbreak takes place, but only for clusters Q3 and Q4, as before. A full disease development is encountered at (BP4), where now clusters Q1 and Q2 show COVID-19 infection. This scenario is clearly depicted in Fig. 5d and f showing high infections for higher p (i.e., for inefficient testing campaigns). Cluster-oriented interpretation can be distinguished by locally targeted testing campaigns (higher p) and widespread random testing campaigns (lower p). Our model thus conjectures that it takes smaller reduction of from in order to clean up the active cases in Q1 and Q2 than in Q3 and Q4.

As can be seen from the numerical study discussed above, both the mobility restriction factor and the case detection ratio p play a crucial role in controlling the disease. For instance, Fig. 5c reveals that the branching point BP1 is responsible for a first COVID-19 outbreak, occurring in clusters Q3 and Q4. Therefore, our next concern will be to investigate how this critical point varies in the p- control space. For this purpose, we will carry out a two-parameter continuation of this critical point, see Fig. 6a. The black curve represents a locus of branching points on the p- plane. The resulting curve divides the control space into two regions: one for stable disease-free equilibria (yellow) and one corresponding to stable endemic equilibria (blue). In this way, for a specific disease control policy represented by the pair , we can determine a priori whether the policy will be effective or not in controlling a COVID-19 outbreak. This can be verified at the test points P1–P4 shown in Fig. 6a. For all these points, test trajectories are calculated using the data shown in Table 1, see Fig. 6b. As can be seen, the solutions computed at P1 and P3 (disease-free region, in yellow) decay to zero, while those computed for P2 and P4 (endemic region, in blue) settle down to an endemic equilibrium, where a long-term COVID-19 outbreak occurs.

Optimization of the disease control policies

In the previous section, we applied numerical continuation methods to study the effect of the mobility restriction factor and the case detection ratio p on the behavior of the modified COVID-19 model (26). In this way, we established critical values of the control parameters upon which a disease outbreak occurs. In this section, we will consider the effect of the control parameters on the average policy effect and the policy cost, as defined in Sec. 5.1. For this purpose, we will assume that the disease control policies represented by the pair are chosen from the yellow region in Fig. 6a.

To begin our study, we will carry out the numerical continuation of disease-free equilibria of model (26) with respect to and monitor the behavior of the average policy effect defined in (28). The result of this procedure can be seen in Fig. 7a. As can be expected, is a decreasing function of , since the average COVID-19 infections that can be avoided (as explained in Sect. 5.1) decrease if higher degrees of mobility between clusters are allowed. Moreover, panel (a) shows a series of points labeled P, which correspond to -values yielding . At these points, the resulting costs are shown in Fig. 7c, which depicts the behavior of the cost function (see (29)) with respect to . As can be seen in the diagram, this function grows as decreases, which is consistent with the fact that stricter contact restrictions lead to higher policy costs. This observation then raises the question if for a desired fixed value of , a more convenient control policy can be found in terms of cost reduction. To tackle this question, we will employ two-parameter continuation with respect to p and to find loci of control points yielding fixed values of , monitoring the corresponding cost function. The result can be seen in Fig. 7b. Here, a family of curves in the p- plane are shown for which is kept fixed. Panel (d) presents the behavior of the cost function along the curves obtained in panel (b). As can be seen, in all cases the cost function presents local minima, which can be interpreted as an optimal policy implementation for a desired fixed value of .

Fig. 7

a One-parameter continuation of equilibria as in Fig. 5c, showing the behavior of the average policy effect and the policy cost (panel c) with respect to . The points P correspond to -values yielding . These are found at (P), (P), (P), (P), (P) and (P). b Two-parameter continuation of equilibria of system (26) with respect to p and , keeping the average policy effect constant (at the values specified above). In this panel, the yellow and blue regions are the same as in Fig. 6a. Panel d shows the behavior of the policy cost computed along the curves obtained in panel (b), for fixed (grey curve), (blue curve), (red curve), (green curve), (purple curve) and (black curve). In this panel, the cost for fixed attains a minimum at (P), (P), (P), (P), (P) and (P)

Concluding remarks

Analysis in this study covers spatio-temporal aspects of COVID-19 transmission in Sri Lanka with the aid of basic reproduction number and path-following continuation pertained to the role of NPIs (contact restrictions and testing campaigns). The daily new cases have been widely used data; however, we processed normalized cases via the populations of RDHS divisions. It suppresses unbiased estimates as higher numbers of cases are reported in highly populated RDHS divisions. Subsequently, all RDHS divisions were categorized using Moran’s scatter into four clusters. Prioritization as well as route for interventions should be Q1 (high-high), Q4 (high-low), Q2 (low-high), and Q3 (low-low). One useful contribution is that the government can use such a route in vaccination programs started at the latter stage of the study period. Priority within a cluster may rely on logistics available within that cluster and temporary shift from the other clusters. Our result is also helpful in placing appropriate border controls for the sake of curtailing transmission waves. Even though Q1 and Q3 do not encounter different incidence levels in their spatial lags, Q2 is vulnerable for significant absorption while Q4 is responsible for significant diffusion. Therefore, border controls can be placed in every important intersecting point between two different clusters.

We extend the qualitative clustering analysis into a quantitative one by conducting an inverse problem using the cases data. A preliminary model of SIURSD type is proposed, carrying the metapopulation context with memory. Due to non-observable model variables, dimensional reduction leads us to an SI type. This final model may look parsimonious; however, it still explains essential mechanistic processes of COVID-19 transmission: cluster-wise contact matrix, viral shedding period, transmission scaling between detected and undetected cases, case detection ratio, contact restriction factor, and loss of immunity. Fitting to the data was done to reveal hidden dynamics including contract matrix, initial conditions for the active cases, case detection ratio, transmission scaling, and loss of immunity. Nonetheless, the SI type may provide beneficiary to big data analytics, especially when the observation period and network size are extended. Forward bifurcation for strongly connected network among clusters was found around the basic reproduction number 1. Numerical investigation was done for the case where the network is, according to rounding small -values to zero, not strongly connected. Time-varying effective local reproduction numbers for all clusters are also presented. Their appearance supersedes clueless cases data when it comes to localizing time at which the current transmission is high (reproduction number greater than 1), suggesting for immediate interventions.

An interesting result is evident from one-parameter continuation of equilibria. Recalling the analytical framework in Sec. 3.3, the initial direction of the continuum of endemic equilibria at is the Perron vector of the next generation matrix (see Eqs. (20) and (21)). As the network associated with the next generation matrix or the contact matrix is not strongly connected (see Fig. 8), Perron–Frobenius Theorem (cf. [80]) only guarantees the nonnegativity of the Perron vector. Particularly to our case, we obtainThis Perron vector indicates two findings: (1) the clusters Q1 and Q2 remain in the disease-free states when shortly exceeds 1; meanwhile (2) the long-term number of active cases in the cluster Q3 jumps to larger extent than that in the cluster Q4 as observed in Fig. 5. If we read the bifurcation diagrams backward in and p, then these findings mean that Q1 and Q2 achieve disease-free states quicker than Q3 and Q4 under the reduction of p and from and , respectively. The network in Fig. 8 explains that Q2 receives a relatively small “injection” from Q1 but returns with a large injection to Q1; meanwhile there is no essential self-injection in Q2. Equipped with a small self-injection, Q1 also injects Q3 and Q4 at comparable rates. Meanwhile, Q3 admits a relatively large self-injection but spares an injection to Q4. On the overall picture, it is arguable that Q1 and Q2 lose endemicity faster than Q3 and Q4 if the entire injection rates (the nonzero entries of the contact matrix ) are reduced simultaneously. At a certain stage, there comes, on the one hand, a scenario where the self-injection in Q1 and thus the injection to Q2 are negligible, making Q2 non-reproductive. On the other hand, the negligible injection from Q1 is compensated by the self-injection in Q3 that withstand both Q3 and Q4 in the endemic states. Notwithstanding this interesting finding, we also observe that disease-free equilibrium (DFE) can be found by reducing p and not so far away from and , respectively. Thus, we argue that the original interventions imposed by the government had been satisfactory during the observation period. From the point of view model transients, one should note that significant contact restrictions are both costly and not gainful in terms of average policy effect (APE). This is evident by concave behavior of APE and convex behavior of the cost against (see Fig. 7(a,c)). Therefore, reducing p and to arbitrarily small values does not make much sense. Scenarios for the optimal values of p and minimizing the cost under fixed magnitudes of APE were proposed. As expected, even optimal results come with a price, as the optimal -values walk toward the third quadrant by increasing APE values; see Fig. 7(b).

Fig. 8

Network based on the contact matrix . The arrow directed from the cluster Qi to the cluster Qj translates the statement “the susceptible humans in Qj contract infection through contact with infected humans in Qi” or shortly “Qi causes infection in Qj”

Finally, this study leaves us some gaps for further improvement. First, several attributes in the original model can be modified to capture more complexities. For example, the average viral shedding period for the undetected cases could have been different from that of the detected cases due to nonoccurrence of symptoms. Despite the averaging, taking the timely proportion of detected cases to be a constant can be too stringent owing to the unknown dark figures (undetected cases). Future improvement may include time-dependent noise for such parameters with given (under guidance of field experts) or computationally tested priors. That recovered and deceased cases preserve a constant ratio is also worth of improvement. Second, the control parameters and p actually represent adjustment of contact restrictions and testing campaigns on the national level, meaning that the scaled susceptible cases in all clusters are enforced the same way toward endemicity reduction, irrespective of their local resources. Meanwhile, the reduction of -values via the cluster-independent also serves as another limitation of the model. From the application point of view, this means that all actions entailed in the contact restrictions should simultaneously follow the adjustment of without proper consideration as to what actions are paramount among the others. For example, reduction of from to 0.75 means that those who go out for activities should reduce the intensities to , those who travel across clusters 4 days a week should change to 3 days a week, schools that are opened 4 days a week should be opened 3 days a week instead, working in the office 5 days a week should change to 3.75 equivalent working days, etc. Technically speaking, these changes may sound simple from the standpoint of the decision maker; however, different abilities and preferences among humans can make the implementation difficult to trace. Third, our SI model contains several parameters that multiply with others. A parameter identification analysis is worth considering if one were to reveal possible dependency among them and thus correct the model specification. Fourth, had regional data of confounding factors been there, we could have integrated these data e.g., in the -values from time to time to capture the different cluster peaks and fluctuations. This is due to the fact that the -values appear to be equivalent to the term of new cases. Incidence and meteorological data from other countries, for example, could be helpful toward this direction.