COVID-19: Estimation of the transmission dynamics in Spain using a stochastic simulator and black-box optimization techniques

Computer Methods and Programs in Biomedicine 211 (2021) 106399 Contents lists available at ScienceDirect Computer Methods and Programs in Biomedicine journal homepage: www.elsevier.com/locate/cmpb COVID-19: Estimation of the transmission dynamics in Spain using a stochastic simulator and black-box optimization techniques Marcos Matabuena a , ∗, Pablo Rodríguez-Mier b , Carlos García-Meixide c , Victor Leborán a a CiTIUS (Centro Singular de Investigación en Tecnoloxías Intelixentes), Universidade de Santiago of Compostela, Santiago de Compostela, Spain b Toxalim (Research Centre in Food Toxicology), Université de Toulouse, INRAE, ENVT, INP-Purpan, UPS, Toulouse 31300, France c Universidade de Santiago de Compostela, Santiago de Compostela, Spain a r t i c l e i n f o Article history: Received 19 April 2021 Accepted 31 August 2021 Keywords: Epidemic models COVID-19 Computing science Stochastic processes Evolutionary computations a b s t r a c t Background and objectives: Epidemiological models of epidemic spread are an essential tool for optimizing decision-making. The current literature is very extensive and covers a wide variety of deterministic and stochastic models. However, with the increase in computing resources, new, more general, and flexible procedures based on simulation models can assess the effectiveness of measures and quantify the current state of the epidemic. This paper illustrates the potential of this approach to build a new dynamic probabilistic model to estimate the prevalence of SARS-CoV-2 infections in different compartments. Methods: We propose a new probabilistic model in which, for the first time in the epidemic literature, parameter learning is carried out using gradient-free stochastic black-box optimization techniques simulating multiple trajectories of the infection dynamics in a general way, solving an inverse problem that is defined employing the daily information from mortality records. Results : After the application of the new proposal in Spain in the first and successive waves, the result of the model confirms the accuracy to estimate the seroprevalence and allows us to know the real dynamics of the pandemic a posteriori to assess the impact of epidemiological measures by the Spanish government and to plan more efficiently the subsequent decisions with the prior knowledge obtained. Conclusions: The model results allow us to estimate the daily patterns of COVID-19 infections in Spain retrospectively and examine the population’s exposure to the virus dynamically in contrast to seroprevalence surveys. Furthermore, given the flexibility of our simulation framework, we can model situations —even using non-parametric distributions between the different compartments in the model— that other models in the existing literature cannot. Our general optimization strategy remains valid in these cases, and we can easily create other non-standard simulation epidemic models that incorporate more complex and dynamic structures. ©2021 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY license ( http://creativecommons.org/licenses/by/4.0/ ) 1. Introduction The spread of SARS-CoV-2 is generating unprecedented health and socio-economical crisis worldwide, being one of the most significant challenges in Europe since World War II. In the light of this emergency, the governments ought to organize an appropriate schedule and optimize political decisions based on scientific evidence to avoid the collapse of the healthcare system, reduce virus-related mortality and minimize the potential effects of an economic recession [39,44,50,51] . ∗Corresponding author. E-mail address: [email protected] (M. Matabuena). Given the vital capacity of the virus to spread and the lack of effectiveness of preventive measures, many countries have been systematically forced to lock down the population temporarily. Although these policies may help control the spread of the virus, they are economically unsustainable over time. In this regard, forecasting the evolution and consequences of the pandemic based on the exposure of the population becomes a critical factor in decision-making [23,40] . However, it is first necessary to assess the current spread of the epidemic to rigorously predict these effects, which is often unknown due to the limited tracking of new infections and active cases. At the beginning of the 20 th century, the first mathematical models to study the dynamics of an epidemic were introduced. Probably the best-known method is the susceptible-infectedhttps://doi.org/10.1016/j.cmpb.2021.106399 0169-2607/© 2021 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY license ( http://creativecommons.org/licenses/by/4.0/ ) M. Matabuena, P. Rodríguez-Mier, C. García-Meixide et al. Computer Methods and Programs in Biomedicine 211 (2021) 106399 recovered model (SIR). SIR model and its variations [33,38] divide the population into compartments, and using differential (deterministic) equations, the number of individuals in each of the compartments over time are estimated. Since then, many new variations of these models that also involved stochastic versions have been introduced in the literature (see for review [3,5,49,70] or other contemporary examples [14,59] ). Despite the enormous progress with these models in recent decades, their direct applications can be limited in several settings. First, most models explain the dynamics of the epidemic at the population level [28,37,41] excluding relevant individual interactions. Second, model-specific assumptions can be restrictive and abstracted from practice. For example, practitioners use Poisson’s homogeneous process to handle the mechanism of new infections or parametric distributions that determine time transitions [29,42] . Third, introducing model reformulations in practice can be challenging and time-consuming with the current optimization strategies of the literature-based primary on designed specific procedures with likelihood equations [10] . We believe this is a critical factor limiting the performance of initial experiments and the use of novel and non-standard formulation of epidemic models in a routine and straightforward manner. As in statistical learning theory, we can say that there is no universal model for all scenarios. Instead, we probably have to design specific models following the existing epidemiological evidence for each situation and introduce the prior knowledge obtained into models. Simulation techniques are a prominent alternative method to build complex and more realistic epidemiological models at a high computational cost. However, their use is not new, and several agent models have appeared in the literature [75,79] , which allow modeling the possible impact of different interventions on the evolution of a pandemic. For instance, we can study the impact of vaccination, social distance, or lockdown policies in the reduction of infections or mortality [32] . More specifically, some of the specific advantages are summarized below: • We can introduce a wide variety of distributions in the components of the model that can be specified with intractable complex likelihood equations [17] or even non-parametric assumptions. • Simulation models allow the introduction of personal information of individuals, such as age and other covariates relevant to disease manifestation, without introducing challenges in the model implementation, unlike classical epidemic models. • Adding some constraints into the model, such as the social interactions between individuals, is not complicated from a model design perspective and only increases computation demands. A cornerstone in expanding this area of research is the ability to obtain reliable solutions to the underlying optimization problem without resorting to problem-specific optimization strategies. Advances in computational power and the field of Black-Box optimization [66] can be an essential milestone in achieving such a goal and being able to examine different models without consuming much time using general procedures. However, sometimes, this strategy requires high-computing environments. Only by evaluating an objective function can these algorithms learn reasonable solutions performing multiple simulations. In this paper, we explore this idea. Using a flexible yet straightforward dynamic probabilistic model that we designed based on the biological evidence of the onset of the pandemic, we estimate the seroprevalence in different regions of Spain along different waves. We also reconstruct the dynamics of infections and recoveries in different compartments to answer specific epidemiological questions, such as when the famous infection peaks happened. To do this, we solve an inverse problem with the mortality records to estimate some specific model parameters, such as the daily rate of infections. In this task, we use, for the first time in this area, the CMAES algorithm [30] , one of the state-of-the-art Black-Box stochastic optimization methods that have been in our previous tests more competitive than other existing algorithms. We must note that our primary purpose in mathematical modeling is not to make forecasts about the dynamic evolution of the pandemic. Instead, the aim of our proposal is to perform backcasting: to retrospectively reconstruct the dynamics of infections while estimating seroprevalence in the different compartments of the model. By estimating this information, we can better characterize the concrete mechanisms of virus transmission in the territories analyzed. Thus, for example, we can guide political decisions in a more refined sense by establishing more advanced and personalized epidemiological thresholds to determine lock-down policies, according to each territory’s specific socio-economic and healthcare factors and the dynamic evolution of the number of infections drawn by our model in the different compartments. 1.1. Outline The article structure is as follows: First, we introduce our new mathematical model to estimate the spread of COVID-19 in regions and countries together with the model optimization strategy used. Then, we introduce some historical background on the evolution of the COVID-19 pandemic in Spain. Also, some demographic and economic characteristics of the Spanish population are presented. Next, we evaluate the behavior of the model and we illustrate its usefulness, performing different analyses across several Spanish regions, reporting the day-to-day evolution of susceptible, infected, and recovered patients. Finally, we discuss the results, the model limitations, and the power and value of the new methodology presented in the existing literature. 1.2. Aims of the analysis In order to show the usefulness and broad potential of our proposal for practitioners, we perform different analyses that allow answering the following epidemiological questions: 1. What was the spread of the virus in the first wave in different regions of Spain like?. For example, when did the peak of infections occur?. How many infected people were there in Spain at the end of the lockdown policies? 2. Using a longer time frame, until March 1, 2021, how were the overall dynamics of SARS-CoV-2 in the Spanish population as a whole?. For example, the healthcare situation was critical in October of 2020, and there were discussions about applying a national lockdown; What could be the real epidemiological situation at that time? 3. Given that, from a theoretical point of view, we can reconstruct the dynamics of infections with our model, how was the actual day-to-day capacity to detect new cases in Spain? 2. Mathematical model and optimization strategy 2.1. Model elements Suppose that D = { 0 , 1 , . . . , n } is the set of days under study. Consider the following random processes whose domain is defined on D. • S(t) : Number of people susceptible to become infected on day t. • I 1 (t) : Number of infected individuals who are incubating the virus on day t. 2 M. Matabuena, P. Rodríguez-Mier, C. García-Meixide et al. Computer Methods and Programs in Biomedicine 211 (2021) 106399 Fig. 1. Diagram of state changes in our model. • I 2 (t) : Number of infected people who have passed the theoretical incubation period and who: (i) don’t show symptoms or (ii) symptoms are mild on day t. • I 3 (t) : Number of infected people who have passed the incubation period and do show moderate or severe symptoms on the day t. • R 1 (t) : Number of recovered cases which are still able to infect on the day t. • R 2 (t) . Number of recovered cases that are not able to infect anymore on the day t. • M (t) : Number of deaths on day t. Henceforth, we will denote by I(t) = I 1 (t) + I 2 (t) + I 3 (t) the number of infected people at time t ∈ Dand R (t) = R 1 (t) + R 2 (t) the number of recovered people. The above random processes describe the dynamic of population individuals in separate compartments. We divided the infected and recovered individuals in a broader and specific taxonomy for the particular case of the COVID-19 than the classical epidemiological models [5,38] . There are two main reasons for this. First, the patients tested by healthcare are usually those found in I 3 . In this case, there is an essential corpus of prior knowledge about how they evolve, and in case of death, their survival time. Second, there is evidence that there are recovered patients who can still infect others. 2.2. Basic model definition The causal mechanism of newly infected individuals is introduced below. For each day t ∈ D, we assume that the new infections I new 1 (t) are generated by the individual interaction of the susceptible people with infected patients and the recovered cases that they can still contaminate (patients that belong of states I 1 , I 2 , I 3 , R 1 ). Formally, we assume that if an individual can contaminate, it does so according to a random variable X ∼Poisson (R i (t)) , being R i (t) the average number of new infections that can cause each person in the day t. It is natural to assume that the function R i (t) follows a decreasing trend in the first months of the epidemic, basically due to two reasons: (i) quarantine policies have been systematically introduced along with different countries and regions. (ii) the number of susceptible people decreases over time, while the number of infected people can increase. These facts indicate that in our particular setting, it is more complicated to interact with non-infected people. Once a new infected person arrives to the model (see Fig. 1 ), we assume that the transitions between the different graph states are modeled by a probability law that verifies the following conditions: (i) the transition probabilities are independent of the absolute instant when such transition takes place Table 1 Random variables of the time of each transition. Transition Random variable Used references I 1 → I 2 Gamma (5 . 807 , 0 . 948) Abdel-Salam and Mollazehi [1] , Lauer et al. [43] I 1 → I 3 Gamma (5 . 807 , 0 . 948) Abdel-Salam and Mollazehi [1] , Lauer et al. [43] I 2 → R 1 Uniform (5 , 10) I 3 → R 1 Uniform (9 , 14) Abdel-Salam and Mollazehi [1] I 3 → M Gamma (6 . 67 , 2 . 55) Abdel-Salam and Mollazehi [1] , Novel et al. [57] , Salje et al. [68] , Verity et al. [81] R 1 → R 2 Uniform (7 , 14) Bi et al. [7] , Ehmann et al. [22] Table 2 Probability of each transition. Coefficient Value Used references α0.8 Day [19] , Mizumoto et al. [53] , Nishiura et al. [56] , Tabata et al. [77] β0.06 Dudel et al. [20] , Fauci et al. [24] , Mahase [48] , Rajgor et al. [65] , Verity et al. [80] , Wu et al. [83] (ii) the probabilities depend only on the current state of the patient regardless of the previous path in the graph. In particular, given the States = {I 1 , I 2 , I 3 , R 1 , R 2 , M} , and α, β∈ [0 , 1] , we have: P (I 2 |I 1 ) = α, P (I 3 |I 1 ) = 1 −α, P (M|I 3 ) = β, P (R 1 | I 3 ) = 1 −β, P (I 3 |I 2 ) = 1 , P (R 2 |I 1 ) = 1 ; all other transitions take a value equal to zero in probability. More schematically, the P , probability transition matrix, between events is shown in the Eq. (1) . P = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ I1 I2 I3 R 1 M R 2 I1 0 α1 −α0 0 0 I2 0 0 0 1 0 0 I3 0 0 0 1 −ββ 0 R 1 0 0 0 0 0 1 M 0 0 0 0 1 0 R 2 0 0 0 0 0 1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ . (1) Additionally, Table 1 shows the random variables that model the time between transitions together with the references used in our elections for real examples of COVID-19 in Spain. Table 2 shows the values used to model the transition probabilities. The Supplementary Material provides specific details about how the mentioned parameters and functions were selected. Finally, as we defined the model above, we have a continuoustime probabilistic model. However, the surrogate variables to fit and compare the model results are recorded daily in real-world situations. Consequently, in our implementation we perform the simulation between transitions and new infections on a daily basis and we truncate the corresponding continuous time in days. 2.3. Stochastic model implementation Our model does not have a closed-form solution. Therefore, in a real-world setting, it is necessary to use statistical simulation methods to approximate specific population characteristics of the stochastic process as quantile functions. Also, we must fit some parameters of the model to characterize the behavior of the study population. For this purpose, we use a sample of the deceased patients {M 1 , M 2 , . . . , M s } along the set of days O = { 1 , . . . , s } . Next, we suppose that our model ( M) is dependent on a vector of parameters θ= (θ1 , θ2 ) ∈ R p 1 ×R p 2 (with p 1 + p 2 = p), where θ1 is a vector of dimension p 1 , defined in beforehand, and θ2 must be estimated from the sample. Furthermore, let us assume that the initial state of the system is characterized by S = (S(0) , I 1 (0) , I 2 (0) , I 3 (0) , R 1 (0) , R 2 (0) , M (0)) ∈ N 7 and T = (T 1 (0) , T 2 (0) , T 3 (0) , T 4 (0) , T 5 (0) , T 6 (0)) ∈ N m ×. . . ×N m . Shas the 3 M. Matabuena, P. Rodríguez-Mier, C. García-Meixide et al. Computer Methods and Programs in Biomedicine 211 (2021) 106399 number of elements for each compartment of the model on day 0. T also contains the amount of remaining days to complete the transition they are in for each individual in the initial state, being m a natural number that represents the maximum number of registered days. To simplify the notation, for each day t ∈ D, we denote the average dead trajectory by the function Mean (θ1 , θ2 , S, T )(t) . The next step is to estimate ˆ θ2 . To do this, we propose to solve the following optimization problem: ˆ θ2 = arg min θ2 ∈ S⊂R p 2 s  i =1 ω i (M i −Mean (θ1 , θ2 , S, T )(i )) 2 , (2) where ω = (ω 1 , . . . , ω s ) is a weighted vector that can help to improve model estimation. Examples of these weights may be: ω i = M i /  s i =1 M i or ω i = (1 / M i ) / (  s i =1 1 / M i ) (i = 1 , . . . , s ) . At this point, it is relevant to note that the above optimization problem (2) , as formulated, includes the possibility of introducing constraints in the parameters’ space. The previous fact is essential because we often have prior knowledge of the range of parameters, or we can even invoke the biological interpretation of parameters to answer this question. Thus, by introducing this knowledge, we can spread up the speed of optimization Black-Block techniques significantly. In Eq. (2) , we have used the real mean trajectory. However, in practice, this is unknown, and we must approximate it using simulation. Next, we run B different simulations, and we denote by M (θ1 , θ2 , S, T ) = 1 B  B i =1 M i (θ1 , θ2 , S, T ) , the estimated mean trajectory. M i (θ, θ2 , S, T ) (i = 1 , . . . , B ) denotes the result of simulation number i . So the optimization problem to be solved is: ˆ θ2 = arg min θ2 ∈S⊂R p 2 s  i =1 ω i (M i −M (θ1 , θ2 , S, T )(i )) 2 . (3) Schematically, the overall optimization process is described below. 1. Define an initial θ0 2 and run B times M(θ1 , θ2 , S, T ) . We denote by M 1 (θ1 , θ0 2 , S, T ) , . . . , M B (θ1 , θ0 2 , S, T ) , different results are obtained. 2. Estimate the mean trajectory M (θ1 , θ0 2 , S, T ) = 1 B  B i =1 M i (θ1 , θ0 2 , S, T ) 3. Estimate the mean square error ˆ RSS 0 =  s i =1 ω i ( M (θ1 , θ0 2 , S, T )(i ) −M i ) 2 . 4. To construct a succession of vectors { θj 2 } R +1 j=1 so that ˆ RSS 0 > ˆ RSS 1 > ˆ RSS 2 > . . . > ˆ RSS R +1 . For example with a stochastic optimization solver. 5. Stop after R + 1 iterations and return θR +1 2 as the optimal parameter of the problem. In our particular setting, θ2 contains the parameter of the R i (t) function that are defined in Section 2.2 . In our preliminary experiments, we assumed that their functional form is equal to R i (t) = min { C, ae −(bt + ct 2 + dt 3 + et 4 + ft 5 ) } where a ∈ [0 , 3] , b ∈ [ −1 , 1] , c ∈ [0 , 1] , d ∈ [0 , 1] , e ∈ [0 , 1] , f ∈ [0 , 1] with θ2 = (a, b, c, d, e, f) ∈ [0 , 3] × [ −1 , 1] ×[0 , 1] ×[0 , 1] ×[0 , 1] ×[0 , 1] and Cis a positive constant fixed 0.005. However, our final election after several experiments and sensitivity analysis with different family of functions that include the previously exponential or inverse sigmoid/logit functions is, R i (t) = d + a 1+ b −(t−c) , where a ∈ [0 , 5] , b ∈ [0 , 1] , c ∈ [ −30 , 30] , d ∈ [0 , 0 . 1] . 2.4. Structural limitations of the model in COVID-19 pandemic The behavior of our model is primarily determined by the parameters α, β, and the function R i (t) , while small variations in the distributions function of the transitions times should not have a significant impact on the seroprevalence estimations. R i (t) is estimated from observed mortality records. However, the parameters α, βare fixed, with statics values over time, and perhaps, in practice, their value should vary in successive waves. In addition, α, β, determine the infection fatality rate ( IF R ), the quotient between fatalities, and the number of infections. In particular, in our model IF R , is given, IF R % = ((1 −α) ∗β) ∗100 . Some studies have shown that in the current pandemic, IF R is the gold standard epidemiological indicator for monitoring the severity with which the virus has affected different countries [58,60] . However, dynamic estimation of IF R is challenging to perform since a precise approximation of the real number of infected people is generally only possible in a single time point, thanks to seroprevalence studies. Several studies have investigated which factors influence the value of the IF R . The primary sources of variations are age and sex [76] . If the distribution of new infections is uniform along time regarding these variables, we can assume constants and statics values for αand βin different time-spans. Testing that assumption is necessary to determine if we must vary the coefficients α, and β, over different waves. However, since we do not know about the true infections in successive waves, it is not trivial to validate this hypothesis empirically, and some statistical estimations are needed. In Spain, heath institutions performed an ambitious and unique longitudinal epidemiological study to know the patterns of virus expansion in several time points that allow us to draw estimation about IF R . In particular, we have information on infections at the beginning of June and the end of November 2020. Using this information, we can estimate the differences in IF R between the first and successive waves. We must note that the Spanish situation in the first wave was critical, with many problems in the elderly population, particularly in nursing homes, and and so changes in the IF R along time are expected. 2.5. Model implementation to handle multi-waves Analogous to an intervention analysis in the context of time series, we need to update the daily infection function to be able to model well the reality in at least the following two situations: at the end of each lockdown, and a new return to normal, or shortly before an explosive growth in the number of new cases or deaths occurs, and no lockdown policies are applied. In these situations, there are abrupt changes in daily infection trends, and therefore the functional form of the function R i (t) needs to be modified. Consider t 0 = 0 < t 1 < t 2 < . . . < t m , m temporal points, defined with expert knowledge, and, in which, we hypothesize that the trend of the daily infection function is modifiable between different { t s } m s =0 , for example between two waves. Then, we propose to define the function R i (t) , as a piecewise function, dependent on the local functions R t 1 i (t) , ... , R t m i (t) , that is, R i (t) = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ R t 1 i (t) t ∈ [0 , t 1 ) R t 2 i (t) t ∈ [ t 1 , t 2 ) . . . R t m i (t) t ∈ [ t m −1 , t m ) . In our fits in successive waves, we assume that the functional form of each R t s i (t) (s = 1 , . . . , m ) is identical, and equal to prior function R t s i (t) = d s + a s 1+ b −(t−c s ) s (for all t ∈ [ t s −1 , t s ) , and any s ∈ { 1 , . . . , m } ), where, the sub-index s , denotes the dependence of parameters, to time-period [ t s −1 , t s ) . Then, the number of initial freemodel parameters is multiplied by the number of periods, m , considered. 4 M. Matabuena, P. Rodríguez-Mier, C. García-Meixide et al. Computer Methods and Programs in Biomedicine 211 (2021) 106399 Regarding parameters αand β, our implementation allows varying these parameters between { t s } m s =0 . Let be αs and βs , the mentioned parameter for the interval [ t s −1 , t s ) . In this paper α= α1 = . . . = αm . However, in the global Spanish analysis along different waves, βwill be vary dynamically. 2.6. Conformal simulation bands to quantify model uncertainty Quantify model uncertainty is a critical point in order to interpretate the results obtained together with their natural limits. In epidemic modeling, according to [21] , we can decompose the model uncertainty in three different sources of error: • Data uncertainty : Uncertainty in specified model parameters, estimated externally from data, or in data to which models are fitted. • Stochastic uncertainty : Uncertainty derived from the method of simulation. • Structural uncertainty : Uncertainty in the optimal model structure, or derived from the use of more than one model structure for a given question. In our settings, we use the dynamic evolution of mortality as a source of information. Then, directly reporting the “Stochastic Uncertainty” drawn for the stochastic simulator as a measure of uncertainty is unrealistic, since we can restrict the number of possible scenarios that happen in practice, according to mortality residuals. More formally, our source of information is a correlated time series that determines the possible reality that occurred, and for a given configuration of the vector parameter θ, we should be removed a significant fraction of unrealistic simulation scenarios. “Data uncertainty” is not trivial to incorporate in this type of model, and the non-parametric bootstrap solution proposed in D’Agostino McGowan et al. [21] have important limitations, since it does not consider the dynamic and correlation structure of daily reports. More general, bootstrap strategies that can handle the specific dependence structure of daily reports can be adapted as our setting, such as block or wild bootstrap, but it is something out of the scope of this work. Testing “Structural uncertainty” is also challenging. There may be no universal model in an epidemic modeling context, and perhaps the best model for each situation is a combination of broad dictionary of simpler models that change dynamically as the pandemic evolves. In this paper, we propose to use a solution similar to the one proposed in Shen et al. [73] , which consists of selecting a small fraction of simulations, according to error criteria, e.g., mean squared error, that shape the possible real infection trajectories. With the remaining trajectories, we apply specific conformal inference techniques that, to the best of our knowledge, have not been applied previously to the context of stochastic simulation models. Conformal inference methods are a general methodology to quantify model uncertainty [72] , with well-established theoretical foundations [45] , and were used in an extensive list of machine learning and statistics problems (see for example a contemporany application [46] ). Below, we introduce the specific mathematical detail of the used conformal simulation bands that can handle heterocedastic noise. First, suppose that θ= (θ1 , ˆ θ2 ) is the optimal parameter configuration, where ˆ θ2 was estimated according to methodology proposed in the Section 2.3 . Then, 1. Perform B = 10 , 0 0 0 simulation of the model M(θ1 , ˆ θ2 , S, T ) and evaluate the mean square error metric (RSS). ˆ RSS s and M s (θ1 , ˆ θ2 , S, T ) ( s = 1 , . . . , B ), denote the results for iteration s of mean square error and the estimation of mortaly records in the simulation model respectively. 2. Let Sel = { i ∈ { 1 , . . . , B } : ˆ RSS i ≤ˆ RSS (10 0 0) } , the set of index of simulation with the lesser or equal 10 0 0 value of RSS estimations. ˆ RSS (10 0 0) denote the element 10 0 0, considering the order sample of { ˆ RSS s } B s =1 . 3. Using the subsample of death simulation trajectories { M i (θ1 , ˆ θ2 , S, T ) } i ∈ Sel , estimate pointwise the standart deviation ˆ σ(t) , ∀ t ∈ O. 4. Define the conformal score Score i = max t∈O | M i (θ1 , ˆ θ2 , S, T )(t) −M t | ˆ σ(t) if ˆ σ(t) > 0 , ∀ i ∈ Sel. Otherwise Score i is equal to 0. 5. Calculate the quantile q α= arg min t∈ R + {  10 0 0 s =1 1 { Score s ≤t} 10 0 0 ≥α} , with α= 0 . 95 , to guarantee distributional intervals that cover a confidence level of 90% . 6. Define for each t ∈ O, the confidence interval prediction as [ M (θ1 , ˆ θ2 , S, T )(t) −ˆ σ(t) q α, M (θ1 , ˆ θ2 , S, T )(t) + ˆ σ(t) q α] , where M (θ1 , ˆ θ2 , S, T ) denote the simulation mean that correspond with our given mortality estimations. 7. Finally, to build confidence bands for the rest of the stochastic process that makes up our epidemic model, we must select simulation trajectories that lead to the mortality outcome falling within the mortality band calculated in step 6). A key element of conformal inference is the selection of conformal scores. In this paper, the conformal score has been estimed using the geometry of the supreme norm || ·|| ∞ . || ·|| ∞ is often used in the analysis of stochastic processes in the field of functional data analysis to estimate confidence bands (see for example [27] ). 2.7. Our probabilistic model proposal in the literature The literature on epidemic modeling is very broad and includes both mechanistic models built from causal epidemiological knowledge and more statistical approaches that exploit information from historical data with purely predictive models [9,41] . Nowadays, it is not easy to establish a boundary between both approaches since, on many occasions, both methodologies are used from the same point of view or even jointly. Our model definition is not very complicated from a mathematical point of view. However, it introduces new challenges from the computational and modeling point of view: the simulation of the trajectory of each individual along the population, the introduction of probability distributions that go beyond the exponential law as traditional models do [5] , and the use of latent infections models through a non-homogeneous Poisson process, extending in this sense the Markovian property [5] . In [25] , the authors define a model similar to ours and propose a resolution framework with Bayesian estimation methods. However, we assume that specific parameters are known from the epidemiological scientific evidence. Our approach using mortality records [61] allow us to obtain reliable seroprevalence estimates, but the difficulty of the estimation increases. We also do not introduce a Bayesian approach but a frequentist approach with its advantages and disadvantages as the need in the Bayesian paradigm of selecting prior functions. Finally, with the philosophy of our simulation model, we can consider complex extensions without making too many changes in the implementation. 2.8. Model optimization with a CMA-ES black-box solver In our setting, we must find optimal parameters taking into account randomness in approximating the mean. To do this, we should resort to stochastic optimization algorithms. Many stochastic algorithms are available in the literature. Still, based on the excellent preliminary results, we have decided to use 5 M. Matabuena, P. Rodríguez-Mier, C. García-Meixide et al. Computer Methods and Programs in Biomedicine 211 (2021) 106399 a state-of-the-art evolutionary algorithm: the CMA-ES [30] . CMAES is an evolutionary-based derivative-free optimization technique that can optimize a wide variety of functions, including noisy functions, like the one we use in our method. One survey of Black-Box optimization strategies found that CMA-ES outranked 31 other optimization algorithms, performing exceptionally well on “difficult” functions or larger dimensional search spaces [31] . From a theoretical point of view, CMA-ES can be seen as a particular case of the Expectation-Maximization algorithm (EM) [8] . We provide specific algorithm steps in Supplementary Material. 2.9. Inverse problem behavior The model fit with the mortality records invokes new challenges in model identification so that the inverse problem is welldefined. We fixed some model parameters according to existing scientific evidence to address this issue to make the problem more regularized. A potential alternative to fitting more parameters with the data is to transform the objective function into a multiobjective optimization problem, considering the daily cases or other ICU indicators as additional sources of information. However, in the early stages of the pandemic, and even nowadays, there were essential doubts about the real capacity of detection of new infections. 2.10. Tuning parameter We performed multiple experiments with CMA-ES to check how the optimization solver behaves. At the same time, through statistical simulation, we estimate the variance of the empirical mean by varying fatalities in different settings. After those initial experiments, we decided to estimate the mean at 300 repetitions, that is, B = 300 . In addition, we have allowed CMA-ES to run 30 0 0 iterations, starting the optimization algorithm from different random points. To obtain the results of this paper, CMA-ES, was able to find the optimal solution with the function R i (t) selected in less than two hours. Finally, the loss function used is Mean Square Error (MSE) with w i = 1 (i = 1 , . . . , s ) (see Section 2.3 ). 2.11. Software details and resources Our proposal has been implemented in several programming languagesC++, Python and Ralthough the results that are shown in this article have been obtained with Python. We optimized the parameters using library pycma [2] , and numpy has been used for mathematical operations. In the different performed statistical analyses, we have used R. Plots have been made both in R with ggplot2 library and in Python with matplotlib . Finally, the training data used to fit the models in the first wave can be downloaded at [67] , and [18] , that represent the daily Spanish statistics of COVID-19 fatalities. In the most extensive analysis of the overall Spanish population, we use the excess of mortality as a source of information. The raw data to estimated excess of mortality can be obtained in the public web interface related to MoModaily Spanish mortality surveillance system ( https://momo.isciii.es/ public/momo/dashboard/momo _ dashboard.html ), coordinated by National Institute of Health Carlos III (ISCIII). We release the code used in this paper for the benefit of the scientific community at ( https://github.com/covid19-modeling ). 3. COVID-19 in Spain Spain was one of the first countries worldwide to experience the effects of COVID-19, after China and Italy. However, the consequences were more dramatic despite the delayed outbreak start with respect to these countries. To give a better context to the evolution of Coronavirus in Spain, in the first wave, and compare it with other countries, we introduce some historical background: • January 31st. The first positive result was confirmed on Spanish territory in La Gomera. At that time, there were around 10,0 0 0 confirmed cases worldwide. • February 12th. The Mobile World Congress, one of the most remarkable technological congresses in the world, to be held in Barcelona, was cancelled. • March 8th. Multitudinous marches were celebrated in Spain. Also, sports competitions and other events were held as usual. • March 13th. Madrid reported 500 new cases of Cov-19 in one day (64 deaths total). Wuhan had gone into lockdown with 400 new cases per day (17 deaths total). • March 14th. With the increase in the outbreak of infections, the government declared a quarantine throughout the country. • March 21st. Due to an overloaded health system, the first patients started to arrive at new makeshift hospitals. • April 3rd. Spain accounts for a total of 117,710 confirmed cases, surpassing Italy for the first time. • April 6th. Spain becomes the country in the world with more deaths per million inhabitants. • April 9th. The FMI forecasts that 170 countries are going to fall into recession this year in the worst crisis since the Great Depression. • April 18th. The Spanish government changes protocols for the daily statistics of COVID-19. Fig. 2 , shows the accumulated number of cases and deceases respectively in the previous periods in Spain, Italy, China, United Kingdom, and the United States according to the data supplied by the different governments. Following the statistics of the Population Reference Bureau, Spain is the 20th country with the world’s oldest population [12] . The country demographic structure, poverty rates, and epidemiological profiles are essential to compare mortality between countries. In the Coronavirus disease, relative and absolute case-fatality risk (CFR) [26] increases dramatically with age and with comorbidity, as evidenced by the current literature. Relative risk can increase by more than 900% in patients over 60 [85] . Subsequently, we perform a descriptive analysis in the regions of Spain that we analyze in this paper: Galicia, País Vasco, Castilla y León, Madrid, Cataluña. Table 3 contains the essential demographic and socioeconomic characteristics of these regions. We can see that Castilla y León is the region with the highest proportion of elderly people. At the same time, Castilla y León has the most delocalized population centres, and the País Vasco is the region with the lowest poverty rate. The national poverty rate is higher than the other analyzed regions because we do not include the most poverty regions. Spain is a multicultural country where there are significant economic, geographical, social, and demographic differences throughout the regions. All these peculiarities make Spain an interesting country to extrapolate the effects of COVID-19 spread to other regions and countries. Finally, in Fig. 3 , we show the evolution of infections and fatalities among the regions under consideration in the first wave. As we can see, Madrid is the most affected region, while Galicia is the least affected, despite its older population. However, it is essential to note that the outbreak began later, and the containment was carried out earlier than in Madrid. 6 M. Matabuena, P. Rodríguez-Mier, C. García-Meixide et al. Computer Methods and Programs in Biomedicine 211 (2021) 106399 Fig. 2. Spread and number of deaths of Coronavirus in Spain, Italy, China, and the United States. Number of accumulated infected patients (left) and the number of accumulated deaths (right) [13] . Table 3 Demographic and socioeconomics characteristics of the Spanish population throughout some regions: Galicia, País Vasco, Castilla y León, Catalu ´ na, Madrid [34] . Galicia País Vasco Castilla León Catalu ´ na Madrid Spain Population 2,698,763 2,181,916 2,553,301 7,609,497 6,685,470 47,100,396 At-risk-of-poverty rate 18.8 8.6 16.1 13.6 16.1 21.5 Population density 91.28 305.19 25.47 239.01 830.02 93.08 Percentage of population by age group 0–9 7.50 8.93 7.51 9.78 9.92 9.28 9–18 7.56 8 . 78 7.79 9.74 9.44 9.37 18–30 10.39 10.66 10.67 12.72 12.80 12.42 30–45 21.18 20.28 19.47 22.24 23.24 22.07 45–60 22.77 23.28 23.60 21.77 22.23 22.61 60–80 22.60 21.41 22.30 18.24 17.33 18.70 from 81 on 8.01 6.67 8.65 5.50 5.03 5.56 Fig. 3. Evolution of accumulated infected (left) and death patients (right) in Galicia, País Vasco, Castilla y León, Madrid, Cataluña. 4. Results 4.1. First wave analysis In order to explore the limits of the model in a more challenging scenario, we start the analysis with the first wave. In this period, most Spanish seroprevalence surveys were performed, and therefore we have a reliable estimation of the number of infections accross different Spanish regions in a single time point, allowing us to evaluate our model performance. In addition, the information is of poor quality, and epidemiological evidence is scarce; thus, making estimations in this scenario is more complicated. More specifically, we restrict the model analysis to April 26st in Galicia, País Vasco, Castilla y León, Madrid and Cataluña. As there is considerable uncertainty about mortality records, we assuming that these two scenarios hold: 1. We assume that the number of real deaths due to Coronavirus is reflected in official records. 7 M. Matabuena, P. Rodríguez-Mier, C. García-Meixide et al. Computer Methods and Programs in Biomedicine 211 (2021) 106399 Fig. 4. Results in Madrid. 2. We suppose a more pessimistic scenario. We assume that many people have died of Coronavirus, but they have not been included in the records because a diagnostic test was not performed. In particular, we shall suppose that the number of deaths is twice as high as those indicated in the official records each day. To display results in an easy-to-view format, we graphically represent the evolution of some states defined at the beginning of Section 2.1 in the two cases considered. In addition, to gain further insights into the results, we show (i) the number of people who may be contaminated or have already transmitted the virus as a percentage of the population size; and (ii) the rate of new infections each day (denoted in the Figures as λt ). Finally, we introduce confidence bands of our estimations using methodology described in Section 2.6 . Here, we only show the Figures that contain results in Galicia and Madrid. The rest of the Figures are available in Supplementary Material. 4.2. Multi-wave analysis To show more recent and informative results on Coronavirus dynamics in Spain, we adjusted the model for the total Spanish population until 1 March 2021. To avoid choosing between the two scenarios above, we use excess mortality as a source of information to feed our model. αhas been selected with the same criteria as the first wave. However, βis fixed with a value equal to 0.085-in the first period, while the rest with 0.0425. These values were established to guarantee an IF R of 1 . 7% in the first wave and 0 . 85% in successive periods. Specific details about IF R estimations are relegated to the Supplementary Material. Daily infection function R i (t) was fitted as a piecewise function (see Section 2.5 for details). In particular, the cut-off points selected for the jumps are as specified below ( Figs. 4–6 ). 1. 1 March to 30 May. 2. 1 June to 5 July. 3. 6 July to 15 August. 8 M. Matabuena, P. Rodríguez-Mier, C. García-Meixide et al. Computer Methods and Programs in Biomedicine 211 (2021) 106399 Fig. 5. Results in Galicia. 4. 16 August to 30 September. 5. 1 October to 24 December 6. 25 December to 1 March. It is important to note that these periods correspond to critical events in pandemic evolution, such as a change in lockdown politics, holidays, or other events that led to abrupt changes in the dynamics of infections. 4.3. Analysis of results The most relevant results in the first wave ( 1 st March 2020 to 26 th April 2020) are outlined below: • Madrid was the most affected region by COVID-19. If we consider an extreme setting (e.g., the number of deaths is double that reported by the Government), 22.5% of the population could have been infected or recovered from the virus. On April 26st, there may have been almost 1,2 million of patients recovered. • Galicia was the region that suffered the mildest effects. The percentage of infected people was less than 2.9%. • Castilla y León, Pais Vasco and Cataluña could have suffered the effects of COVID-19 with a proportional magnitude. In those regions, the percentage of infections could have been between 6 and 12% of the population. • The peak of new infections probably occurred at the start of quarantine, while the peak of people who can contaminate took place between March 17 and 24. • The number of new infections have been dramatically reduced after the introduction of containment measures. • The most accurate scenario is the pessimistic scenario. The analysis of the excess mortality reported in the Supplementary 9