Mathematical modelling of the second wave of COVID-19 infections using deterministic and stochastic SIDR models.

Recently, various countries from across the globe have been facing the second wave of COVID-19 infections. In order to understand the dynamics of the spread of the disease, much effort has been made in terms of mathematical modeling. In this scenario, compartmental models are widely used to simulate epidemics under various conditions. In general, there are uncertainties associated with the reported data, which must be considered when estimating the parameters of the model. In this work, we propose an effective methodology for estimating parameters of compartmental models in multiple wave scenarios by means of a dynamic data segmentation approach. This robust technique allows the description of the dynamics of the disease without arbitrary choices for the end of the first wave and the start of the second. Furthermore, we adopt a time-dependent function to describe the probability of transmission by contact for each wave. We also assess the uncertainties of the parameters and their influence on the simulations using a stochastic strategy. In order to obtain realistic results in terms of the basic reproduction number, a constraint is incorporated into the problem. We adopt data from Germany and Italy, two of the first countries to experience the second wave of infections. Using the proposed methodology, the end of the first wave (and also the start of the second wave) occurred on 166 and 187 days from the beginning of the epidemic, for Germany and Italy, respectively. The estimated effective reproduction number for the first wave is close to that obtained by other approaches, for both countries. The results demonstrate that the proposed methodology is able to find good estimates for all parameters. In relation to uncertainties, we show that slight variations in the design variables can give rise to significant changes in the value of the effective reproduction number. The results provide information on the characteristics of the epidemic for each country, as well as elements for decision-making in the economic and governmental spheres.

Introduction

Since the end of 2019, the world has been experiencing the consequences of the increase in cases of COVID-19. In political, economic, and social context, various measures have been adopted distinctively by the countries, in an attempt to mitigate the effects of the pandemic [45, 54]. After a significant increase in the number of infections and deaths, the proposed measures had an effect to a certain extent, causing the number of new infections to be reduced in some locations. However, the relaxation of such measures, in view of the supposed control of the spread of the disease, caused the number of cases to rise again [20, 29, 34]. Other factors, such as the identification of reinfections and the seasonal immunity, can also be associated with such an increase [15, 25, 33, 53]. In the sense of understanding the burden caused by the spread of diseases, numerous studies in different fields of science have been conducted, among which we may mention especially the development and improvement of mathematical models [42], environmental studies [3, 19, 30, 48, 49] and geopolitical aspects [46].

Some relevant works that analyze the early stages of the spread of the disease have helped to understand important aspects of the pandemic. In this context, Perc et al. [42] proposed a flexible model to predict the spread of COVID-19 based on minimal and maximal growth rates (for the number of infected individuals) in a certain period for the USA, Slovenia, Iran, and Germany. Hâncean et al. [21] studied the initial spread of COVID-19 in Romania, presenting an analysis of imported cases (mainly Romanian citizens returning from Italy) using human-to-human transmission networks. Ribeiro et al. [46] analyzed the relationship between the size of cities and the spread of COVID-19 in Brazil, using power-law models relating the size of the population and the number of cases/deaths. The results indicate that small cities are more affected in the early stages of the pandemic; on the other hand, large cities are more favorable environments for the spread of SARS-CoV-2 in the long term.

The analysis of the effects of environmental factors and the use of spatial data infrastructure, such as the Geographic Information System (GIS), have also been explored in research on the spread of COVID-19. Roher et al. [48] considered the influence of weather scenarios—for instance, dust of particulate matter—in the spread of COVID-19 in some locations in Europe, such as Paris, London, and the Canary Islands. The influence of aerosols on the infection fatality rate of COVID-19 in Germany, Italy, and Spain was studied by Li et al. [30], focusing on the aerosol optical depth. In the same vein, Gupta et al. [19] analyzed the effect of weather on the spread of the disease in India with a detailed statistical approach. Using a support vector machine approach, Gupta et al. [19] established an association between daily transmission and weather parameters. Al-Kindi et al. [3] employed GIS techniques to investigate the spatiotemporal pattern of the COVID-19 pandemic in the territory of Oman. Stoffel et al. [49] analyzed the effect of the pandemic on carbon dioxide emissions during 2020, presenting an interesting comparison between solar radiation management and an antipyretic medicine (referring to global warming as a fever). The authors warn that, until the end of the pandemic, governments must do even more to reduce greenhouse gas emissions.

The so-called second wave of a pandemic is characterized by an expressive increase in the number of cases (and consequently also in the number of deaths caused by the disease), after a significant drop in the number of new infections during the first wave. Research show that the virus associated to COVID-19 has been mutating [44]. However, the occurrence of new cases, as well as a possible increase in the severity of the disease or in the transmissibility, is not yet associated with such biological alteration [14]. According to Ghanbari [18], the conditions for the occurrence of a second wave of infections are due to the increase in the number of cases within specific groups, such as people in the same age-group, in specific places, or even the persistence of the spread of the first wave.

Vaccines have been widely studied since the onset of the pandemic [4, 7, 17, 32]. Several countries have started vaccination, but so far, only a few have achieved high vaccination coverage, as the UK [28] and Israel [35]. Libotte et al. [31] proposed a strategy to determine an optimal control strategy for vaccine administration in COVID-19 pandemic considering a hypothetical scenario in relation to the Chinese population. Mukandavire et al. [39] estimated the basic reproduction number and critical vaccination coverage to control the disease for different hypothetical vaccine efficacy scenarios in South Africa. However, some countries face a gradual increase in the number of cases of the disease, while vaccination is still a long way off. Priesemann et al. [45] claimed for a pan-European effort to mitigate the effects of further waves of COVID-19, aiming at (i) reducing the number of cases, (ii) keeping the number of cases at a low level, and (iii) develop a longer-term common vision to guide the management of the pandemic.

In view of the threat of a new wave of infections, several efforts have been concentrated on trying to predict their behavior and thus provide mechanisms for decisions aimed at mitigating the disease. In this context, Cacciapaglia et al. [8] proposed a robust approach based on epidemic renormalization group framework to simulate the second wave for different countries. The authors conducted statistical analyses on different levels of human interaction to obtain an estimate of the occurrence of the second wave in Europe, which would most likely take place between July 2020 and January 2021. Nori Junior et al. [40] proposed a two-wave statistical model based on the superposition of normal distributions, where future scenarios can be predicted based on parameters estimated for the first wave. Castro [9] proposed a SEIR-based model with time-dependent transmission and mortality rates to represent the second wave considering data from Spain and Germany. The proposed model was employed to simulate the post-confinement epidemic under several scenarios. Ghanbari [18] proposed a SID-based model with two sub-populations (which are classified by immunological criteria) of infected individuals, in order to simulate the second wave in Iran. Pedro et al. [41] proposed a SEIR model that considers a time-varying parameter to evaluate the influence of schools and workplaces closure on the transmission of SARS-CoV-2.

Usually, such models and related studies do not consider the uncertainties associated with simulations and data. In this context, stochastic models have been widely used as a way to incorporate such uncertainties into computational simulations. Chanu and Singh [10] studied a stochastic SEQIR model to analyze the epidemic in India. In turn, He et al. [22] proposed a discrete-time stochastic SEIHR model considering the governmental measures to control the spread of the disease, as well as the clinical condition of individuals. Rihan et al. [47] incorporated a time delay in the stochastic SIRC model. Faranda and Alberti [16] proposed a stochastic SEIR model to analyze the impact of the relaxation of isolation measures in France and Italy. Tang et al. [51] proposed a stochastic SEIHRQ model to investigate the imported cases on local epidemic in the Shaanxi-China province. This model takes into account the effects of intermittent inflow with a Poisson distribution on the likelihood of multiple outbreaks.

Typically, compartmental models have constant parameters, which may not be sufficiently descriptive, affecting the goodness-of-fit statistics and, consequently, the long-term simulations. In addition, such models may not be able to describe the behavior of the second wave of an epidemic, without taking advantage of arbitrary choices. The major contributions of this work are to provide a generalized methodology that is capable of identifying the interface between two epidemiological curves, so that compartmental models with time-dependent parameters can fit to more complex epidemiological data accurately. Our goal is to determine the parameters of the SIDR model that best describe the behavior of the second wave in Germany and Italy. We adopt a piecewise time-dependent transmission rate for each wave. The determination of the parameters of the model considers a constraint related to the effective reproduction number, in order to provide more realistic results. In addition, we propose to analyze the uncertainties of the simulations using a stochastic SIDR model, taking into account the incorporation of additional differential equations to represent the variation of some parameters. This work is organized as follows. The description of the SIDR model is presented in Sect. 2, considering the deterministic approach and a stochastic version of the model. The formulation of the proposed inverse problem is presented in Sect. 3. Section 4 describes the data employed, the numerical methods used to solve the deterministic and stochastic problems and the definition of some error metrics—used to evaluate the performance of the proposed approach. The results obtained by solving the inverse problem are presented and discussed in Sect. 5. Finally, conclusions are drawn in Sect. 6.

Epidemiological models

In this section we describe the deterministic and stochastic versions of the SIDR compartmental model used in the mathematical modelling, as well as the definition of the initial conditions and some considerations regarding the basic reproduction number and the effective reproduction number. For the stochastic SIDR model, we present the stochastic differential equations that describe the uncertainties in the control parameters.

Deterministic SIDR model

In order to simulate the dynamic behavior of COVID-19 epidemic, we adopt the SIDR model. This model considers the interaction among individuals in the population, which are classified as susceptible, infected, dead and recovered individuals. We employ a normalized version of the model, in the sense that the quantities S, I, D and R are divided by the population size N. The variations of S, I, D and R over time in this model are represented byParameters , and represent, respectively, the transmission rate, death rate, and recovery rate. The initial conditions are , , , and , since the constraint given by holds for all t. The term represents the first element in the vector of cases reported (i.e., the number of infected individuals at ), as we will discuss in Section 3. In turn, is a quantity related to the population size and the fraction of cases reported. We use a new variable to represent I(0),since this parameter will be an argument of the optimization procedure in the inverse problem.

In order to predict the potential spreading of an infectious disease in a population, the basic reproduction number () may be considered as an epidemiological metric. This parameter represents the average number of secondary cases resulting from a single case, taking into account a population totally susceptible to the disease [5]. When , the disease has a spreading potential, with a tendency to increase the number of new cases.

The basic reproduction number for the SIDR model is given by [5]It can be a large value when is small compared to . Thus, in the early stage of COVID-19 spread in different countries, we can observe high values for considering different numerical approaches [8, 24, 26, 27, 36, 43, 57]. As mentioned by Marimuthu et al. [36], this parameter encompasses the rate of contact for the individuals in the population, the probability of transmission by contact (not every contact is capable to transmit the disease) and the estimated period in which the population is affected by the disease. On the other hand, a population may not be entirely susceptible to the disease, as consequence of immunity, for instance. In practice, the effective (or time-dependent) reproduction number () can provide more realistic results. represents the expected number of secondary cases arising from a primary case of infection at time t [52].

To evaluate this parameter, various methodologies can be used, such as exponential growth and time-dependent methods [36]. In this work, we adopt defined as [13]where S(t) is the susceptible population at a given time t.

Stochastic SIDR model

In the context of modeling the dynamics of epidemics, there are several factors that incorporate uncertainty in the parameter estimation process. Considering such uncertainties means obtaining results that are consistent with the real scenario. In order to evaluate the influence of uncertainties in such epidemiological models, the control parameters can be stochastically perturbed [16] and, for this purpose, the original deterministic model is rewritten as a stochastic model, where the parameters that characterize the dynamics of the epidemic are defined by specific probability distributions. In this case, for each control parameter (in our case, , and ), a new differential equation with stochastic contribution is formulated. Thus the following stochastic differential equations are employed [16]:where is the stochastic control parameter, which in our case corresponds to . In turn, is the initial condition related to the corresponding parameter, is the standard deviation for the stochastic control parameter, and is a random number represented by a specific distribution type. In this case, the stochastic SIDR model is given bywith appropriated initial conditions (discussed in Sect. 4.2.2).

Formulation of the inverse problem

In the SIDR model described by Eq. (1), the parameters , , and must be calculated in order to allow the numerical solution of the system of ordinary differential equations. It is important to mention that since is a quantity related to the effective number of infected individuals (see Eq. (2)), it is employed to define the initial conditions for susceptible, dead, and recovered individuals, i.e., , and , respectively. For this purpose, it is necessary to formulate and solve an inverse problem through the minimization of the difference between calculated values and reported data.

Commonly, parameters , and are taken to be constant during the integration of the differential system, for the sake of simplicity. However, according to Castro [9], both the transmission rate and death rate can be chosen as time-dependent parameters to better capture the effects of policies defined to mitigate the pandemic. Following Cheynet [12], we define as:where (), () and (day) are parameters that define the transmission rate of a particular epidemic.

The other parameters of the model are considered to be constant. Thus, for each wave, the parameters , , , and must be calculated. In addition, for the first wave, the quantity related to the initial value for the infected individuals () must be estimated, as well as the time instant where each set of parameters are effective, i.e., the interface between the end of the first wave and the start of the second wave. The first wave is defined within , and the second wave holds for , where is the final time. It is important to mention that the model continuity is guaranteed at .

We will define a vector of parameters, (to be estimated), aswhere the superscripts and refer, respectively, to the first and the second waves.

The numerical solution of the system given by Eq. (1) for a discrete time set and a particular approximation of the vector produces vectors of responses for each compartment of the model, i.e., , , and . We adopt reported data related to infected and dead individuals, which are respectively denoted by and , in the parameter estimation procedure. Then, we are interested in an approximation of the vector that minimizes the error between the model responses and , and the reported data, and , respectively. In other words, we aim to minimize the quantitywhere and are the highest values in the vectors and (just for scaling purposes).

Furthermore, in practice, the basic reproduction number and the effective reproduction number cannot assume arbitrarily large values. Considering this situation, we introduce an inequality constraint to ensure that the solution of the inverse problem will represent the adherence to the observed values for in the pandemic, as followswhich holds . S(t) is the simulated normalized susceptible population at time t and is the maximum value that the effective reproduction number can take (which is further discussed). We also consider that , according to Eq. (7), indicating the time dependency for . Obviously, the use of or even as a parameter to analyze the spread of the disease must be carefully conducted [2]. The main idea here is to avoid unrealistic values for .

Mathematically, the inverse problem is formulated as [38]:subject to Eq. (9), where is the optimal value of .

In this case, the normalized SIDR model must be simulated considering the parameters calculated by the optimizer, in order to obtain the number of infected and dead individuals estimated by the model and, consequently, the value of the objective function (). It is important to mention that the constrained optimization problem is transformed into an equivalent unconstrained problem by using the static penalty function method [55].

Data, numerical algorithms and error metrics

Data

The data employed for parameter estimation are the daily number of infected individuals and the cumulative number of deaths caused by COVID-19 in Germany and Italy. We use 297 records for each time series, both ranging from February 02, 2020 to November 25, 2020 [56].

Numerical algorithms

Specific algorithms for deterministic and stochastic models are presented below.

Deterministic approach

Each inverse problem is solved by using the differential evolution algorithm [50] (DE) in the optimization step, considering the following parameters: population size , crossover probability , amplification parameter , maximum number of generations , and strategy rand/1/bin (see Storn and Price [50] for a detailed description of the algorithm). The search space is given by , , , , , , and . The evolutionary process is interrupted if a given maximum number of generations is reached.

To evaluate the inequality constraint, is considered to be equal to 10, based on recent studies regarding in Germany [27]. The penalization parameter (required to use the static penalty function method) is equal to . The deterministic SIDR model is numerically solved using the fourth-order Runge–Kutta method [6] considering 1000 equally spaced points. In addition, for each wave, the model is simulated for equal to 297 days.

Stochastic approach

To solve a stochastic differential equation, several numerical methods with different characteristics can be found in the literature [23]. In this work, we use the Milstein method [37]. We do not estimate the parameters of the stochastic model, but only use it for simulation purposes, and therefore, the optimization procedure is not employed in this step. In order to assess the influence of uncertainties, we consider the results regarding the first wave, which are previously obtained using the deterministic approach. In this case, we investigate the sensitivity of the calculated parameters only considering the second wave (when ), whereas the first wave is simulated considering the parameters estimated by means of the deterministic approach. In practice, this means that, for the second wave, the uncertainties are associated only to , , and .

In this context, to simulate the stochastic SIDR model given by Eq. (6) for the second wave, the initial conditions for S, I, D and R are taken as the simulated values at considering the estimated parameters for the first wave (this strategy ensures the continuity of the waves). In addition, the initial conditions for , and are defined as the average values for the first wave denoted as , , and , respectively.

For each parameter, it is required to define a type of distribution to represent the stochastic contribution. As suggested by Faranda and Alberti [16], and are normally distributed, and follows a lognormal distribution, defined aswhere is a proportionality constant, and , , and are the corresponding standard deviations, which are essentially defined as a proportion of the average value of the parameter obtained in the first wave, that is, , and . We also consider in the Milstein method to simulate the process, with final time equal to 600 days and 100 independent runs to obtain the average values.

Error metrics

In order to compare the quality of solution, two metrics are considered in this work. The first is the coefficient of determination (), used to determine the goodness-of-fit of the model. Mathematically, this metric is given by [11]where is the mean value of and or D, or , or . The vector represents a vector of ones with same length of . In practice, the quantity gives the percentage variation in Z compared to the set of variables considered.

The second metric represents the root-mean-squared error (), which is calculated as [11]Similarly to Eq. (12), or . This metric quantifies the differences between the reported data and the responses of the model.

Results

In this section, we present the results for the proposed inverse problem considering the deterministic approach, as well as the analyzes regarding the stochastic approach, for both Germany and Italy.

Deterministic Inverse Problem

The first problem addressed is the solution of the deterministic inverse problem for the two waves. The results are shown below.

Germany

Initially, we are interested in verifying the model’s ability to capture information about both waves using the proposed approach. To that end, we run DE for 20 times, in such a way that each new result obtained by solving the deterministic inverse problem is independent of the others, that is, the initial population is different in each run, and the randomness of the genetic operators guarantees the diversity of the results. This strategy ensures that favorable results are not obtained by chance. Table 1 presents the best result and the corresponding standard deviation for each estimated parameter. In view of the value of the standard deviation corresponding to each parameter, it is clear that the proposed methodology is robust in the sense of estimating the model parameters for the analyzed data set. The visual inspection of Fig. 1a, which shows the simulation of the model compared to the reported data, demonstrates that the set of parameters is in good agreement with the dynamics of the disease. The value obtained for the objective function, evaluated at the best estimator for (the best results in Table 1), is equal to . This corroborates the similarity between the profiles resulting from the integration of and , with the analyzed data. The error metrics indicate that, considering the reported data on the number of infected individuals, . For data on the number of dead individuals, . In relation to the root-mean-squared error, is equal to 546.89075 and 334.28045, respectively for infected and dead individuals.

Table 1

Results obtained for the proposed inverse problem using reported data from Germany. The values shown correspond to the best result and the standard deviation calculated from the computed values

First wave		Second wave
β0fday^-1	0.31349±0.00244	β0sday^-1	0.09517±0.00425
β1fday^-1	0.02190±0.00122	β1sday^-1	0.00013±0.00075
β2fday	0.10140±0.00989	β2sday	0.00033±0.00003
γ^fday^-1	0.00290±0.00011	γ^sday^-1	0.00058±0.00001
α^fday^-1	0.04556±0.00352	α^sday^-1	0.00854±0.00043
δ	1.91×10^-5±2.97×10^-8	τday	166.18015±5.96786

Fig. 1

a Simulation of the SIDR model with parameters obtained from the solution of the inverse problem given by Eq. (10), using epidemiological data from Germany in the period between February 2 and November 25, 2020. The readings for compartments S, I and R are on the left and for compartment D are on the right; b Transmission rate of the disease in the same period, where the vertical blue dashed line represents the transition between the first and second waves

Figure 1a shows the drop over time in the population of susceptible individuals, for each wave. In addition, between the two waves, it is observed a phase where this group of individuals remains nearly constant, i.e., the number of new infections between the end of the first wave and the start of the second wave remains very low. On the other hand, the number of recovered individuals increases during the epidemic period, as expected. Such results demonstrate that the system presents a behavior that is consistent with the dynamics of an epidemic.

Regarding , whose behavior is shown in Fig. 1b, initially it assumes the maximum value within the analyzed range, since at the beginning of the epidemic, the number of immune individuals is very small and, therefore, new infections are more likely to occur. As the time goes by and the gradual mitigation of the disease, this value decreases, tending to zero, showing the proximity of the end of the first wave. During the second wave, this parameter is nearly constant since and, therefore, (see Eq. 7). A possible criticism to this approach is related to the discontinuity in the profile of at . It must be borne in mind that this behavior is a consequence of the choice of a “double-wave” model, where two profiles of are expected. Moreover, with this assumption we avoid discontinuities in the profiles for the individuals in the compartments.

The effective reproduction number , with being calculated using the values shown in Table 1, is presented in Fig. 2. In general, for each wave, decreases significantly compared to the initial value. The high value at the initial point represents the start of the epidemic period. Thus, as different measures are taken to mitigate the epidemic, the natural tendency of this parameter is to reduce until the epidemic is contained . In addition, Fig. 2 shows the behavior of obtained by Petrova et al. [43] using the SID model, and by Khailaie et al. [27] using the SECIR model, both for the first wave. The results are quite similar, mainly in terms of the behavior of the curves. On average, we obtained approximately equal to 1.5 for the first wave and 6.50 for the second wave. According to the Robert Koch Institute [1], preliminary information indicates that, for the second wave, this value is close to 3. Such difference may be due to type of approach used to determine each parameter. The value calculated here is obtained considering a particular set of reported data.

Fig. 2

Estimate of the effective reproduction number of the disease in Germany, in the period between February 2 and November 25, 2020. The values are calculated using Eq. (4) and considering the transmission rate shown in Fig. 1b. The results are compared with the estimates obtained with SECIR [27] and SID [43]

Italy

For the data from Italy, we adopt the same methodology described for the German case: we run the optimizer for 20 times and obtain a set of independent results, of which the best result and the standard deviation corresponding to each parameter are shown in Table 2. Again, it is clear that the simulated profiles, considering the best parameters obtained, have good agreement with the analyzed data, as illustrated in Fig. 3. Qualitatively, the results are similar to those obtained for the data from Germany (see Fig. 1), but differ in terms of the magnitude of the data and results. The error metrics point out that the fit between the responses of the model and the data is also satisfactory, with and for data related to the number of infected individuals, and and for dead individuals. Even with a greater deviation than in the results for the German case, the fitting still reflects the behavior observed in the data set fairly well.

Table 2

Results obtained for the proposed inverse problem using reported data from Italy. The values shown correspond to the best result and the standard deviation calculated from the computed values

First wave		Second wave
β0fday^-1	0.47815±0.02454	β0sday^-1	0.78693±0.02434
β1fday^-1	0.02830±0.00112	β1sday^-1	0.00558±0.00009
β2fday	0.38865±0.06656	β2sday	0.17129±0.02332
γ^fday^-1	0.00502±0.00012	γ^sday^-1	0.00071±0.00001
α^fday^-1	0.05472±0.00155	α^sday^-1	0.02349±0.00001
δ	1.09×10^-6±1.64×10^-9	τday	187.08035±11.74343

Fig. 3

a Simulation of the SIDR model with parameters obtained from the solution of the inverse problem given by Eq. (10), using epidemiological data from Italy in the period between February 2 and November 25, 2020. The readings for compartments S, I and R are on the left and for compartment D are on the right; b Transmission rate of the disease in the same period, where the vertical blue dashed line represents the transition between the first and second waves

In turn, Fig. 4 presents the behavior of for each wave, considering the parameters shown in Table 2. Again, our results are in good agreement with those obtained by Petrova et al. [43] and Khailaie et al. [27] using the SID and SECIR models, respectively. This shows that the proposed methodology is capable of obtaining good estimates for the model parameters and, consequently, for , as illustrated in Fig. 4. In this case, the average values regarding the effective reproduction number are, approximately, 1.31 and 4.7 for the first and second waves, respectively. This result is in line with that obtained by Faranda and Alberti [16] (the basic reproduction number obtained by these authors is, approximately, equal to 4, using a stochastic SEIR model). Due to the stochastic nature of this model, the value of this parameter can fluctuate up to 6. In turn, Cacciapaglia et al. [8] reported that the maximum value that can reach is close to 10 for the Italian case.

Fig. 4

Estimate of the effective reproduction number of the disease in Italy, in the period between February 2 and November 25, 2020. The values are calculated using Eq. (4) and considering the transmission rate shown in Fig. 3b. The results are compared with the estimates obtained with SECIR [27] and SID [43]

Stochastic analysis

In order to analyze the influence of uncertainties, this section considers the simulation of the stochastic SIDR model. For this purpose, at first we consider the results obtained for the first wave, in both cases (Germany and Italy). Our focus is to analyze the dynamics of the population and the behavior of the parameters of the model for the second wave, when . The first wave is simulated with the same parameters obtained in the previous sections (Tables 1 and 2 for Germany and Italy, respectively). In practice, for the second wave, uncertainties are incorporated into the parameters , , and . The initial conditions for the stochastic SIDR model are those obtained for the first wave using the deterministic approach. In this sense, for Germany, we have , , and ; for Italy, , , and . In addition, since is a time-dependent function, its average value is taken in this analysis.

Fig. 5

a Stochastic simulations of the SIDR model for the COVID-19 epidemic in Germany and Italy. The variation of the mean of the distribution corresponding to the transmission rate in the second wave is considered, according to Eq. (11), by assigning different values for ; b Variation of the effective reproduction number over time for each particular value of represented by the area in light blue. The curve shown in dark blue represents the average values for each time instant

Figure 5 presents the average profiles and the effective reproduction number (only for the second wave) considering the influence of uncertainties, with different values for , where , as defined in Eq. (11). The fluctuation of each curve in relation to the deterministic solution (where is the reference profile) is shown in Fig. 5a. For both Germany and Italy, the increase in the value of implies in more distancing with respect to the reference profile, i.e., the peak localization referring to the second wave gets further from the reference value. This variation also affects the maximum value observed for each peak. We may also note that a high value for at a given time represents an increase in the number of susceptible individuals (for both countries), and thus requires a proportional decrease in the populations of infected, recovered and dead individuals, according to the hypotheses of the SIDR model.

In addition, for both cases analyzed, decreases over time (see Fig. 5b). Such profiles represent only the second wave and, therefore, it is expected a high value in the start of the wave (which will produce the second peak, as discussed previously). Furthermore, the higher is the value of , the lower is the initial effective reproduction number for the second wave, which is consonant to the displacement of the second peak to the right in the number of infected individuals, i.e., with a more quantity of susceptible individuals—as a consequence of a higher value of —we obtain a higher value for (for the same time instant t).

Figure 6 provides another perspective on the variation of the average profiles and the values of in Germany and Italy considering uncertainties, as shown in Fig. 5. Initially, note that for , the mean value of is strictly decreasing, given the particular choice of . Such monotonicity provides an indication of the behavior of the average profiles at a given time : the number of infected and dead individuals tends to decrease, whereas the number of susceptible and recovered individuals increases, preserving the total size of the population in both countries, as increases. Clearly, the shift of the average profiles is reflected in the value of the effective reproduction number, as can be seen by comparing the corresponding frames in Fig. 5.

Fig. 6

Variation of the mean of the distribution in relation to the variation of . The gray squares show the four values of used in the stochastic simulations (see Fig. 5) and, for each of these values, 2000 samples of are drawn, where the dashed horizontal lines represent the corresponding mean of the distribution for the particular value of

The behavior of the profiles for different values of in Fig. 5a is in agreement with the mean of the samples in Fig. 6. Even with the values of being equally spaced, the profiles in Fig. 5a do not shift linearly as increases. Note that in Eq. (11) and, therefore, the mean and standard deviation of the samples are coupled and cannot be analyzed individually. In addition, follows a lognormal distribution, which is not symmetrical around the mean value. Therefore, the behavior of the profiles is a consequence of the adopted distributions. Also note that values greater than can produce negative transmission rates (see the samples in frame 4 of Fig. 6).

Conclusions

In this work we proposed a framework to simulate the second wave of COVID-19 pandemic. We analyze data from Germany and Italy, which have well-defined behavior regarding the second wave. In addition, a sensitivity analysis considering a stochastic SIDR model is conducted. For this purpose, the parameters of the SIDR compartmental model, piece-wise defined (for each wave), are obtained using DE. The inverse problem considers a piece-wise time-dependent transmission rate for each wave and an inequality constraint is adopted, in order to avoid unrealistic values for the effective reproduction number. In general, the parameters obtained for the first wave, especially regarding the effective reproduction number, are in agreement with those reported by other authors.

In relation to uncertainty analysis, we can conclude that small perturbations in the design variables vector can result in significant variations in the value of the effective reproduction number. In this case, this fluctuation can help to understand different scenarios of the pandemic, as well as to define policies to mitigate the disease. Finally, the formulation of the inverse problem with an inequality constraint, associated with a piece-wise time-dependent transmission rate for each wave and the sensitivity analysis represent the main contributions of this work.