# COVID-19: Estimating spread in Spain solving an inverse problem with a
probabilistic model

Marcos Matabuena  CiTIUS (Centro Singular de Investigación en Tecnoloxías
Intelixentes), Universidade de Santiago of Compostela, Santiago de Compostela,
Spain<EMAIL_ADDRESS>Carlos Meijide-García  Universidade de Santiago
of Compostela, Santiago de Compostela, Spain  Pablo Rodríguez-Mier  Toxalim
(Research Centre in Food Toxicology), Université de Toulouse, INRAE, ENVT,
INP-Purpan, UPS, Toulouse, France  Víctor Leborán  CiTIUS (Centro Singular de
Investigación en Tecnoloxías Intelixentes), Universidade de Santiago of
Compostela, Santiago de Compostela, Spain

## Abstract

We introduce a new probabilistic model to estimate the real spread of the
novel SARS-CoV-2 virus along regions or countries. Our model simulates the
behavior of each individual in a population according to a probabilistic model
through an inverse problem; we estimate the real number of recovered and
infected people using mortality records. In addition, the model is dynamic in
the sense that it takes into account the policy measures introduced when we
solve the inverse problem. The results obtained in Spain have particular
practical relevance: the number of infected individuals can be $17$ times
higher than the data provided by the Spanish government on April $26$ $th$ in
the worst case scenario. Assuming that the number of fatalities reflected in
the statistics is correct, $9.8$ percent of the population may be contaminated
or have already been recovered from the virus in Madrid, one of the most
affected regions in Spain. However, if we assume that the number of fatalities
is twice as high as the official numbers, the number of infections could have
reached $19.5$%. In Galicia, one of the regions where the effect has been the
least, the number of infections does not reach $2.5$ %. Based on our findings,
we can: i) estimate the risk of a new outbreak before Autumn if we lift the
quarantine; ii) may know the degree of immunization of the population in each
region; and iii) forecast or simulate the effect of the policies to be
introduced in the future based on the number of infected or recovered
individuals in the population.

##  1 Introduction

The spread of the Coronavirus is generating an unprecedented crisis worldwide.
In Europe, for example, it is considered to be the most significant challenge
that the continent has faced since World War II. In the light of this
situation of emergency, the governments must work quickly to avoid the
collapse of the healthcare system, reduce the mortality associated with the
virus, and avoid the possible effects of an economic recession [1, 2].

Given the strong capacity of the virus to spread and the lack of preventive
measures, many countries have been systematically forced to lockdown the
population temporarily. Although these policies may be useful in controlling
the spread of the virus in the immediate future, they are unsustainable over
time from an economic point of view. 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 [3, 4]. However, to rigorously
predict these effects, it is necessary to assess the current spread of the
epidemic, which is often unknown.

At the beginning of the XX century, the first mathematical models to study the
dynamics of an epidemic were introduced. Probably the most well-known is the
susceptible-infected-recovered model (SIR). SIR models and its variations [5,
6] divide the population into compartments, and using differential
(deterministic) equations, the number of individuals in each of the
compartments over time is estimated. Since then, many new variations of these
models have been introduced in the literature (see for review [7, 8, 9, 10] or
more contemporary examples [11, 12]).

However, those approaches suffer from a critical shortcoming. In essence, they
focus on explaining the dynamics of the epidemic at the population level and
exclude the complex interactions that occur at the microscopic level between
individuals [13, 14, 15]. Furthermore, these models tend to be adjusted with
fixed parameters that are not updated over time, according to how the epidemic
evolves and the different policies introduced.

In the current era of precision medicine, where more efficient solutions are
being sought that optimize the health of the individual [16], it is surprising
not to find in the current state of the art models of epidemic spread that
follow more individual approaches [17, 18], as is rightly pointed out in [19].
Precisely in [19], the authors propose to adopt new models that exploit the
enormous individual information recorded by biosensors [20] and other devices.

A crucial aspect to consider in the success of mathematical modeling is the
need to adjust the parameters of the models with accurate data. However, this
is difficult to achieve in practice, because governments are often overloaded,
and information is only recorded on patients with symptoms. The observational
nature of this data [21], together with a delay or errors in the recording of
information [22], requires specific techniques to correct the biases related
to the data [23].

###  1.1 Main contributions and results

Motivated by the need to developed new models that overcome these drawbacks
and that estimate the expansion of the SARS-CoV-2 in a realistic way in Spain,
we introduce a new probabilistic model. The contributions and most significant
results of this paper are introduced below.

  * •

To the best of our knowledge, we propose one of the few probabilistic models
that estimate the evolution of an epidemic from each individual of the
population. This model combines the finite Markov model, along with continuous
probability distributions and a dynamic Poisson model.

  * •

The model is designed in such a way that it does not need to use statistics on
the number of infections that have a significant measurement error, and uses
mortality records instead.

  * •

The model estimates for each day the number of susceptible, contaminated,
recovered, and fatalities. For this, we solve an inverse problem, and we
introduce biological expert knowledge related to the SARS-CoV-2 disease into
the probabilistic model about the time each individual spends in the different
states and the transition of probabilities.

  * •

Using the data from several regions of Spain, we estimate the real number of
infections as of April $26$st. The results show a higher proportion of
infections than that provided by the Spanish government.

  * •

The effect of the Coronavirus has been uneven among regions of Spain. In
Madrid, more than $10$ percent of the population was infected, or it recovers
from the virus, while in Galicia, this percentage is less than $2.5$.

  * •

These results highlight differences between regions; i) the degree of immunity
of the population can be different; ii) the number of recovered people who can
return to normality is also different; iii) based on this analysis, a
reasonable strategy to be considered could be to lift the quarantine on a
case-by-case basis based on the degree of exposure of each region to the
virus.

  * •

Finally, we repeat the estimations assuming that the number of fatalities is
twice than reported by the official records. In this case, the number of
infections can increases between $40$ y $110$ percent by region. However, we
must interpret these results with caution.

###  1.2 Outline

The mathematical details of the model are relegated to the end of the article,
to facilitate the readability for a general target audience. The structure of
the article is as follows: First, we describe the evolution of Coronavirus in
Spain, together with some demographic and economic characteristics of the
Spanish population. Subsequently, we fit our model in different scenarios, and
we show the evolution of susceptible patients, infected and recovered ones
over time. Next, we discuss the results and limitations of the model. Finally,
we introduce the mathematical content of our model and the computational
details of the implementation and parameter optimization.

##  2 COVID-19 in Spain

Spain was one of the first countries in the world to experience the effects of
COVID-19 after China and Italy. However, despite the delay in the start of the
outbreak, with respect to these countries, the consequences are now more
dramatic. To give a better context to the evolution of Coronavirus in Spain
and to compare it with other countries, we introduce some historical
background:

  * •

January $31$st. The first positive result was confirmed on Spanish territory
in La Gomera. At that time, there were around $10,000$ confirmed cases
worldwide.

  * •

February $12$th. The Mobile World Congress, one of the most remarkable
technological congresses in the world to be held in Barcelona, was canceled.

  * •

March $8$th. Multitudinous marches were celebrated in Spain. Also, sports
competitions and other events were held as usual.

  * •

March $13$th. 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 $14$th. With the increase in the outbreak of infections, the government
declared a quarantine throughout the country.

  * •

March $21$st. Due to an overloaded health system, the first patients started
to arrive at new makeshift hospitals.

  * •

April $3$rd. Spain accounts for a total of $117,710$ confirmed cases,
surpassing Italy for the first time.

  * •

April $6$th. Spain becomes the country in the world with more deaths per
million inhabitants.

  * •

April $9$th. The FMI forecasts that $170$ countries are going to be into
recession this year in the worst crisis since the Great Depression.

  * •

April $18$th. The Spanish government changes protocols for the daily
statistics of COVID-19.

  * •

April $21$th. The president of the United States, Donald Trump, declared:
”It’s incredible what happened to Spain, it’s been shattered.”.

Figure 1, shows 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.

Figure 1: 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) [24].

| Galicia | País Vasco | Castilla y León | Cataluña | Madrid | España
---|---|---|---|---|---|---
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$

Table 1: Demographic and socioeconomics characteristics of the Spanish
population throughout some regions: Galicia, País Vasco, Castilla y León,
Cataluña, Madrid [25].

Following the statistics of the Population Reference Bureau, Spain is the 20th
country with the world’s oldest population [26]. The country demographic
structure, poverty rates, and epidemiological profiles is essential to compare
mortality between countries. In the Coronavirus disease, relative and absolute
case-fatality risk (CFR) [27] 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$ [28].

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 1 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 centers, and the País
Vasco is the region with the lowest poverty rate. 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 the spread of the
Coronavirus to other regions and countries.

Figure 2: Evolution of accumulated infected (left) and death patients (right)
in Galicia, País Vasco, Castilla y León, Madrid, Cataluña

Finally, in Figure 2, we show the evolution of infections and fatalities among
the regions under consideration. As we can see, Madrid is the most affected
region, while Galicia is the least affected, despite its older population.
However, it is important to note that the outbreak began later, and the
containment was carried out earlier than in Madrid.

##  3 Results

Next, we show the estimation performed with the model until April $26$st in
Galicia, País Vasco, Castilla y León, Madrid, Cataluña, assuming that these
two scenarios hold:

  1. 1.

We assume that the number of real deaths due to Coronavirus is the one
reflected by the official records.

  2. 2.

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.

We graphically represent along time the result of estimating how many
individuals belong to each of the following states:

  * •

$I_{1}(t)$: Number of infected individuals who are incubating the virus on day
$t$.

  * •

$I_{2}(t)$: Number of infected people who have passed the incubation period
and do not show symptoms on day $t$.

  * •

$I_{3}(t)$: Number of infected people who have passed the incubation period
and do show symptoms on day $t$.

  * •

$R_{1}(t)$: Number of recovered cases which are still able to infect on day
$t$.

  * •

$R_{2}(t)$. Number of recovered cases which are not able to infect anymore on
day $t$.

  * •

$M(t)$: Number of deaths on day $t$.

In addition, to understand the meaning of the results, we show: i) the number
of people who may be contaminated or who have already transmitted the virus as
a percentage of the population size; and ii) the rate of new infections each
day.

Finally, we introduce confidence bands of our estimations 111Our confidence
bands shall not be confused with classic frequentist statistics confidence
bands.

###  3.1 Scenario 1

#### Galicia

Figure 3: Results in Galicia: we assume that there are as many fatalities as
in the official records.

#### Castilla y León

Figure 4: Results in Castilla y León: we assume that there are as many
fatalities as in the official records.

#### País Vasco

Figure 5: Results in País Vasco: we assume that there are as many fatalities
as in the official records.

#### Madrid

Figure 6: Results in Madrid: we assume that there are as many fatalities as in
the official records.

#### Cataluña

Figure 7: Results in Cataluña: we assume that there are as many fatalities as
in the official records.

###  3.2 Scenario 2

#### Galicia

Figure 8: Results in Galicia: we assume there are twice as many deaths as the
official records.

#### Castilla y Leon

Figure 9: Results in Castilla y León: we assume there are twice as many deaths
as the official records.

#### País Vasco

Figure 10: Results in País Vasco: we assume there are twice as many deaths as
the official records.

#### Madrid

Figure 11: Results in Madrid: we assume there are twice as many deaths as the
official statistics.

#### Cataluña

Figure 12: Results in Cataluña: we assume there are twice as many deaths as
the official statistics.

###  3.3 Analysis of results

The most relevant results are described below:

  * •

Madrid is the most affected the region by COVID-19. $19.5$ % of the population
could have been infected or recovered from the virus if we assume that the
number of deaths is double that reported by the Government. At present, there
may be $120,000$ patients recovered.

  * •

Galicia is the region which suffered the mildest effects. The percentage of
infected people is less than $2.5$ %.

  * •

Castilla y León, Pais Vasco, Cataluña could have suffered the effects of
COVID-19 with a proportional magnitude. The percentage of infections could be
between $6-10$ % of the population.

  * •

The peak of new infections has probably occurred at the start of quarantine,
while the peak of people who can contaminate took place between March 17-24.

  * •

The number of new infections have been dramatically reduced after the
introduction of containment measures. Now, it seems that the situation is back
under control.

##  4 Discussion

The Coronavirus pandemic is causing an unprecedented crisis in health,
economic, and social terms throughout the world. To help designing better
politicies in the future, estimating the number of people who were infected by
the virus is a crucial priority. To address this issue, we have proposed a new
probabilistic model to estimate the spread of an epidemic by simulating the
behaviour of each individual in a population.

### COVID-19 in Spain

The predictions obtained over several regions of Spain have practical
relevance. The results reveal that the number of infections is much higher
than what is reported by the Spanish government. Two factors may explain these
discrepancies: i) The tests were carried out in a limited way among the
general population; and ii) no random sampling was done among the different
regions throughout the territory shading the real magnitude of the pandemic.
It is also worth noting that its impact was uneven across regions. Madrid was
the most affected region, while Galicia was the least affected. However,
according to our model, in Cataluña, Castilla y León, and the País Vasco, the
effects were surprisingly similar, even though the outbreak of infections
theoretically started at different time points. In addition, the geographical
dispersion of Castilla y León is more significant (see Table 1). Perhaps the
older population of Castilla y León (see Table 1) has a higher case fatality
rate than the other communities, and our model is overestimating the number of
infected people in that region. As previously mentioned, taking into account
the demographic structure of the regions, together with the epidemiological
profiles and other socioeconomic characteristics, is fundamental in the
analysis of results.

#### More deaths than those reported in the official records

Even though many patients die due to the virus, this is not the official cause
of death. This problem occurs due to a lack of diagnostic tests. Spanish and
International press echoed this problem in the current Coronavirus crisis, and
more than twice as many deaths than those officially announced have been
reported in certain regions. A recent study reported that the actual number of
deaths might be even three times higher in Italy [29].

Due to this fact, it is vital to take into account this uncertainty present in
the models to perform an accurate estimate of the number of actual infections.
However, this uncertainty is not probabilistic. To alleviate this problem, we
have assumed that there could be up to twice the number of official deaths.

If we consider this correction, the number of infected increases from $40$ to
$110$ percent by region.

#### Practical considerations of this results

These findings allow us have a better sense of the degree of exposure to the
virus of the population and to design a restoration to normality based on the
following factors.

  * •

The degree of immunity of each population.

  * •

The number of people who might contaminate today.

  * •

The number of recovered patients.

In this sense, it is vital to gradually de-escalate the confinement and always
proceed according to the individual characteristics of each region.

For the moment, there is no available vaccine for use with humans, and there
is uncertainty about treatments. Therefore it is advisable to take measures
such as social distancing [30] to prevent the infection rate from skyrocketing
and to avoid another outbreak.

#### The notorious peak

Wondering about when the peak of infection would have been reached was a hot
topic in mass media. However, such peak can be i) the maximum peak of new
infections. ii) the maximum peak of people who can contaminate. In the first
case, our models prove the existence of this peak between $9-16$ March, while
in the second case, it was estimated to be in March $17-24$.

From a perspective of virus containment, the most important is the second one
from as it indicates the number of potential transmitters.

### Methodological aspects of epidemic modeling

From a methodological standpoint, indeed, there are already simulators that
work at the individual level [31, 32, 33, 34, 35]. However, from our point of
view, the models analyzed have some of the following potential limitations: i)
some models use fixed parameters without prior optimization according to the
target population; ii) the parameters do not usually evolve dynamically
according to the policies introduced; and iii), many of these models were
designed without adding biological expert knowledge.

More generally we can find two philosophies when addressing epidemic
forecasting in the literature (see [18]): i) Mechanistic models: Using an
epidemiological approach, these models explain the evolution of an epidemic
from the causal mechanisms of transmissionc [13]; ii) Phenomenological models:
From entirely data-driven approaches such as regression models or time series,
the development of epidemics is predicted using historical data [18].

Phenomenological models can be an appropriate solution in the case of
epidemics such as influenza [36], for which there exists long-term data
collected on hospitalization and mortality rates. Additionally, the etiology,
together with the previous incidence rates, are better characterized. With all
this, extrapolations of the real impact of the epidemic in the whole
population can be done from data-driven approaches. In addition, current
evidence indicates that influenza incubation times are shorter than those of
COVID-19, which reduces an individual’s potential exposure to the virus. For
all these reasons, we believe that in the case of the new COVID-19, it is
necessary to use Mechanistic models as our proposal.

Often, defining a clear boundary between Mechanistic models and
Phenomenological models may not be simple, and even both methods can be seen
as complementary approaches. For example, SIR models and their variants from a
theoretical point of view try to explain the effects of an epidemic from the
mechanisms of transmission, however they are often used as time series [37].
At this point, we would like to point out that the underpinning assumptions of
these models are often unrealistic [5]. Moreover, from the perspective of
Epistemology and Epidemiology, the epidemic spread is based on the complex
interactions between individuals, which these models cannot capture due to
their population-based nature [18].

The use of individual models should be an aspiration in the development of new
methods. In the current era of big data, in developed countries, a substantial
part of the population is geolocalized through mobile devices. Moreover, it is
increasingly common to use biosensors that monitor patient health [20].
Knowledge about interactions between individuals and patient health can
provide a more realistic picture of the spread and evolution of the epidemic
[38]. However, in this context, it is necessary to systematically apply
techniques that correct the biases that appear in the observational nature of
these data. At the same time, it is also important to preserve the privacy of
the citizens.

The absence of reliable information is one of the main problems with the
COVID-19 pandemic. To alleviate this problem, when fitting model parameters,
it is important to understand the social, economic, and demographic
characteristics along with the epidemiological profiles of the region to
introduce expert knowledge into the model and reduce uncertainty in
statistical learning. Also, with this strategy, we reduce the computational
demands and we increase model interpretability.

### Collaborative attitude

Information is one of the most valuable weapons in combating a pandemic of
this magnitude. Understanding the mechanisms of virus transmission [39, 40],
the prevalence of asymptomatic patients in the population [41], risk and
prognostic factors [22] are topics that can reduce the number of deaths. Due
to the absence of prior knowledge of all these factors, a collaborative
approach between countries and institutions to share their knowledge should be
a top priority.

At an international level, we could highlight the initiative promoted by Eran
Segal from Israel, to which several countries have joined [42]. Segal and his
team designed a survey in which anonymous volunteers participated to estimate
the real evolution of the Coronavirus in Israel [43]. This questionnaire was
adapted and applied in other countries. There is now an international
consortium in action to share the results of these surveys, together with
other clinical information from patients. However, it is important to note
that the data obtained from the questionnaires do not represent a random
sample of the population, which limits the quality of the data collected.

### Expert knowledge in parameter fitting

From a general perspective of inverse problems and statistical learning, the
estimation of the model parameters is an ill-posed problem [44, 45, 46, 47,
48, 49]. In the case of epidemic spread models, these difficulties are even
more significant (see, for example, [50, 51]). In practice, the evolution of
several compartments is predicted simultaneously, and usually, there is only
information about the actual growth of one of the compartments.

In this regard, it is vital to follow Ockham’s Razor rule: the simplest
solution is most likely the right. Reducing the number of variables to be
fitted in the models and introducing expert knowledge into these models guided
by current biological evidence is for ill-posed estimation problems.

In the case of COVID-19, the current epidemiological evidence is inconsistent.
It is especially noteworthy about the case fatality rate (CFR) and the
proportion of asymptomatic patients [52]. In the first case, variations
between $0.4$ and $15$ percent have been reported in the literature.
Meanwhile, with asymptomatic patients, these range from $20$ percent to $80$
percent. The primary rationale for these discrepancies is the application of
non-observational data analysis methods for a sample composed largely of
elderly patients. In this sense, after carrying out our estimates taking into
account the demographic structure of Spain and using the published articles
that we believe reflect the reality in a more precise way, we estimate that
the rate of asymptomatic patients ranges between $75$ and $85$ percent.
Likewise, CFR may vary between $0.4$ and $0.8$ percent. A recent article has
estimated that asymptomatic patients to be $78$ percent [53].

### Combining a probabilistic model with a deterministic inverse problem

Our probabilistic model is non-stationary and has a complex dependency
structure. In addition, the time series of each day’s death records can be
seen as a realization of the stochastic process.

To fit the model, we estimate a deterministic inverse problem that makes the
observed deaths is the average of the random process. From there, we define a
window, and we try with a non-probabilistic approach to capture the
trajectories of fatalities close to the real one. In this way, we can obtain
trajectories of the random process that allow us to explain the real situation
of the infected in each region. In other words: From day one, with the chosen
configuration of parameters, we can have many different scenarios in our
probabilistic model; however, at the present moment, we already observed a
concrete trajectory. Then, using this trajectory as a reference, we look for a
set of close trajectories that measure other alternative scenarios that
approximate the real number of deaths.

### Model limitations

The main limitation of our proposal is that some of the parameters on which
the model depends on are set and may need to be tuned more carefully to the
study population. However, the current biological evidence does not allow
further refinement. Also, optimizing more parameters with a data-driven
approach is dangerous due to the quality of data about infections reported by
governments. Therefore, it is better to restrict oneself to current biological
expertise information.

The models are an approximation of reality, and as Cox would say: all models
are wrong, some are useful. We are predicting the unknown in this epidemic
without reliable data. However, with statistical simulation models [54, 55],
we can see the effect of a deviation of the input parameters, the extreme
cases, and the model predictions in other scenarios such as the end of
containment policies.

Another limitation of our model is that despite simulating data for each
individual in the population, we were not using information about the
demographic structure and epidemiological profile along regions. We also not
introduce geolocation information as in the interaction between each
individual in the population [38]. However, the inclusion of all these
aspects, may not lead to greater accuracy as no reliable model output data are
available, and, furthermore, we would be entering into privacy issues that are
beyond the scope of this study.

### Future work and model extensions

Our immediate work is to extend our estimations to other countries. Also, we
want to simulate the effect of an immediate end coronavirus lockdown [56] in
each of these location. We also want to carry out simulations to find out how
to minimize the effects of a possible outbreak in Autumn [4], for example,
with social distance [57] or weeks of partitioned work: combining remote
working with face-to-face activity. In a complementary way, it would be
interesting to predict the effects of these measures on the economy. Some
economists claim that the economic consequences of this crisis may be more
damaging to our health in the long term than the effects caused directly by
the virus [58].

From a methodological point of view, another exciting extension is to develop
a mixed model or an empirical Bayes framework (see an example of this
methodology [59]) to make a simultaneous estimation across many regions and to
be able to make a statistical inference about the parameters of the models. In
this case, it may be interesting to optimize the parameters with Bayesian
optimization [60], to alleviate the computational issues of this approach.
Another variation of the model may be to employ other probability
distributions such as the Conway-Maxwell-Poisson distribution [61] to handle
underdispersion situations in the number of newly infected. However, our
current election of Poisson distribution is well motivated given the
asymptotic properties of the Erdos-Rényi model [62].

Finally, another possible improvement of the model could be to introduce an
ensemble forecasting as they do in meteorology and manage better the possible
uncertainty in the initial conditions [63].

### Our future prediction in the world

Given the current emergency, and the need to provide predictions on the spread
of the epidemic worldwide, we have decided to set up a website
(https://github.com/covid19-modeling/forecasts) where we regularly update our
estimates across many regions of the world in a fully accessible manner. In
addition to this, we will be incorporating our simulations of what the effect
of current end coronavirus lockdown would be in each of these regions.

##  5 Mathematical model

###  5.1 Model elements

We suppose that $D=\\{0,1,\cdots,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$.

  * •

$I_{2}(t)$: Number of infected people who have passed the incubation period
and do not show symptoms on day $t$.

  * •

$I_{3}(t)$: Number of infected people who have passed the incubation period
and do show symptoms on day $t$.

  * •

$R_{1}(t)$: Number of recovered cases which are still able to infect on day
$t$.

  * •

$R_{2}(t)$. Number of recovered cases which are not able to infect anymore on
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\in D$ and $R(t)=R_{1}(t)+R_{2}(t)$ the number of
recovered people.

The above random processes describe the evolution of population in different
compartments. However, unlike the classics models [6, 10], we divided the
infected and recovered individuals in a more broader and realistic taxonomy
for the particular case of the COVID-19. There are two main reasons for this:
First, the patients tested by healthcare are usually those found in $I_{3}$.
In this case, there exists a corpus of prior knowledge about how they evolve
over time and in case of death, which is their survival time. Second, there is
evidence that there are recovered patients who can still infect others.

###  5.2 Model definition

$I_{1}$$I_{3}$$R_{2}$$I_{2}$$R_{1}$$\alpha$$1-\beta$$\beta$$1-\alpha$ Figure
13: Diagram of state changes in our model.

The causal mechanism of new infected individuals (see Figure 13) is introduced
below. For each day $t\in D$, we assume that the new infected population
$I_{1}(t)-I_{1}(t-1)$, is generated by the individual interaction of the
susceptible people with infected patients and the recovered cases who can
still contaminate. Formally, we assume that if a person can contaminate, it
does so according to a random variable
$X\sim\mathbb{\text{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 important
to remark, that it is natural to assume that the function $R_{i}(t)$ follows a
decreasing trend over time in our setting, basically for 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. This all makes it more
complicated to interact with a non-infected people.

Once a new infected person arrives to the model (see Figure 13), 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 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 $\alpha,\beta\in[0,1]$, we
have: $\mathbb{P}(I_{2}|I_{1})=\alpha$, $\mathbb{P}(I_{3}|I_{1})=1-\alpha$,
$\mathbb{P}(M|I_{3})=\beta$, $\mathbb{P}(R_{1}|I_{3})=1-\beta$,
$\mathbb{P}(I_{3}|I_{2})=1$, $\mathbb{P}(R_{2}|R_{1})=1$; all other
transitions take a value equal to zero in probability. More schematically, the
transition matrix between events is shown in the Equation 1.

$P=\begin{bmatrix}&$I1$&$I2$&$I3$&$R1$&$M$&$R2$\\\
$I1$&0&1-\alpha&\alpha&0&0&0\\\ $I2$&0&0&0&1&0&0\\\
$I3$&0&0&0&1-\beta&\beta&0\\\ $R1$&0&0&0&0&0&1\\\ $M$&0&0&0&0&1&0&\\\
$R2$&0&0&0&0&0&1\\\ \end{bmatrix}.$ (1)

Additionally, the Table 2 shows the random variables that model the time
between transitions together with the references that have been used. In
particular cases, as survival time, we have made our estimates based on data
from some of these researchers and incorporating some expert knowledge of
others.

Transition | Random variable | Used references
---|---|---
$I_{1}\to I_{2}$ | $Gamma(5.807,0.948)$ | [64, 65]
$I_{1}\to I_{3}$ | $Gamma(5.807,0.948)$ | [64, 65]
$I_{2}\to R_{1}$ | $Uniform(5,10)$ |
$I_{3}\to R_{1}$ | $Uniform(9,14)$ | [65]
$I_{3}\to M$ | $Gamma(6.67,2.55)$ | [66, 67, 65, 68]
$R_{1}\to R_{2}$ | $Uniform(7,14)$ | [69, 70]
Table 2: Random variables of the time of each transition

Finally, in a similar way, in Table 3, we show the values used to model the
transition probabilities.

Coefficient | Value | Used references
---|---|---
$\alpha$ | $0.8$ | [53, 71, 72, 73]
$\beta$ | $0.06$ | [74, 52, 75, 76, 77, 78]
Table 3: Probability of each transition

###  5.3 Model implementation

The proposed model is not known to have an closed-form solution. In a real-
world setting, it is necessary to use statistical simulation methods to
approximate the mean trajectory or quantile functions over time. 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}$,$\cdots$, $M_{s}$ of each day in the timeframe
$O=\\{1,\cdots,s\\}$.

Next, we suppose that our model ($M$) is dependent on a vector of parameters
$\theta=(\theta_{1},\theta_{2})\in\mathbb{R}^{p_{1}}\times\mathbb{R}^{p_{2}}$
(with $p_{1}+p_{2}=p$), where $\theta_{1}$ is a vector of dimension $p_{1}$,
defined in beforehand, and $\theta_{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))\in\mathbb{N}^{7}$
and
$T=(T_{1}(0),T_{2}(0),T_{3}(0),T_{4}(0),T_{5}(0),T_{6}(0))\in\mathbb{N}^{m}\times\cdots\times\mathbb{N}^{m}$.
$S$ has the 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 integer number
that represents the maximun numbers of registred days.

To simplify the notation, from now on we denote the average dead trajectory by
the function $Mean(\theta_{1},\theta_{2},S,T)(t)$ for each day $t\in D$

The next step is to estimate $\hat{\theta}_{2}$. To do this, we solve the
following optimization problem:

$\hat{\theta}_{2}=\arg\min_{\theta_{2}\in
S\subset\mathbb{R}^{p_{2}}}\sum_{i=1}^{s}\omega^{i}(M_{i}-Mean(\theta_{1},\theta_{2},S,T)(i))^{2},$
(2)

where $\omega=(\omega^{1},\cdots,\omega^{s})$ is a weighted vector that can
help to improve model estimation. Examples of this weights may be:

$\omega^{i}=M_{i}/\sum_{i=1}^{s}M_{i}$ or
$\omega^{i}=(1/M_{i})/(\sum_{i=1}^{s}1/M_{i})$ $(i=1,\dots,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
space of parameters. It is important because we have prior knowledge of what
the range of the parameters is. We may know constraints involving them as
well.

In Equation (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$
differents simulations and we denote by
$\overline{M}(\theta_{1},\theta_{2}^{0},S,T)=\frac{1}{B}\sum_{i=1}^{B}M^{i}(\theta_{1},\theta_{2},S,T)$,
the estimated mean trajectory. $M^{i}(\theta,\theta_{2},S,T)$ $(i=1,\cdots,B)$
denote the result of each simulation.

So the optimization problem to be solved is:

$\hat{\theta}_{2}=\arg\min_{\theta_{2}\in
S\subset\mathbb{R}^{p_{2}}}\sum_{i=1}^{s}\omega^{i}(M_{i}-\overline{M}(\theta_{1},\theta_{2},S,T)(i))^{2}.$
(3)

Schematically, the overall optimization process is described below.

  1. 1.

Define an initial $\theta_{2}^{0}$ and run $B$ times
$M(\theta_{1},\theta_{2},S,T)$. We denote by
$M^{1}(\theta_{1},\theta_{2}^{0},S,T)$, $\cdots$,
$M^{B}(\theta_{1},\theta_{2}^{0},S,T)$ to each of the obtained results.

  2. 2.

Estimate the mean trajectory
$\overline{M}(\theta_{1},\theta_{2}^{0},S,T)=\frac{1}{B}\sum_{i=1}^{B}M^{i}(\theta_{1},\theta_{2}^{0},S,T)$

  3. 3.

Estimate the mean square error
$\hat{RSS}^{0}=\sum_{i=1}^{s}\omega^{i}(\overline{M}(\theta_{1},\theta_{2}^{0},S,T)(i)-M_{i})^{2}$.

  4. 4.

To construct a succession of vectors $\\{\theta_{2}^{j}\\}^{M+1}_{j=1}$ so
that $\hat{RSS}^{0}$ $>$ $\hat{RSS}^{1}$ $>$ $\hat{RSS}^{2}$ $>$ $\cdots$ $>$
$\hat{RSS}^{M+1}$. For example with a stochastic optimization algorithm.

  5. 5.

Stop after $M+1$ iterations and return $\theta_{2}^{M+1}$ as the optimal
parameter of the problem.

In our particular setting, $\theta_{2}$ contain the parameter of a $R_{i}(t)$
function defines in Section 5.2. From now on, we asume that
$R_{i}(t)=\min\\{C,ae^{-(bt+ct^{2}+dt^{3}+et^{4}+ft^{5})}\\}$ where
$a\in[0,3]$, $b\in[-1,1]$, $c\in[0,1]$, $d\in[0,1]$, $e\in[0,1]$, $f\in[0,1]$
with
$\theta_{2}=(a,b,c,d,e,f)\in[0,3]\times[-1,1]\times[0,1]\times[0,1]\times[0,1]\times[0,1]$
and $C$ is a positive constant fixed $0.005$.

###  5.4 Model optimization

In our setting, we have found optimal parameters taking into account
randomness in the approximation of the mean. To do this, we must resort to
stochastic optimization algorithms.

Many algorithms of this kind can be found in literature, but based on the good
results obtained in a preeliminary analysis we have decided to use a state-of-
the-art evolutionary algorithm: the CMA-ES [79]. CMA-ES is an evolutionary-
based derivative-free optimization technique that can optimize a wide variety
of functions, including noisy functions as the one we use in our method. One
survey of Black-Box optimizations found it outranked 31 other optimization
algorithms, performing especially keen on “difficult functions” or larger
dimensional search spaces [80]. From a theoretical point of view, CMA-ES can
be seen as a particular case of Expectation-Maximization algorithm (EM) [81].

###  5.5 Model Inference

Our model has been fitted independently for each region. Therefore, we cannot
make statistical inference in the usual sense since we only observe a
trajectory of a stochastic process with a complex dependence structure. Hence,
to identify the model, we assume that the observed trajectory of accumulated
deaths is the mean of our random process.

Given a day $t^{\prime}\in O$, $M_{t^{\prime}}$ accumulated number of observed
deaths in $t^{\prime}$ and the multivariate random process
$\\{X(t)=(I_{1}(t),I_{2}(t),I_{3}(t),R_{1}(t),R_{2}(t),M(t))\\}_{t\in D}$, the
death trajectories of the probabilistic model close to the observed path on
$t^{\prime}$, are for us a set of possible explanations of the spread of the
coronavirus in a study region. We now formalize this idea. Let
$I_{\alpha}=[M_{t^{\prime}}-\alpha M_{t^{\prime}},M_{t^{\prime}}+\alpha
M_{t^{\prime}}]$, where $\alpha\in(0,1)$. Considerer,
$\\{\omega\in\Omega:M(t^{\prime},\omega)\in I_{\alpha}\\}$, being $\Omega$
sample space of process $\\{X(t)\\}_{t\in D}$. We build our confidence bands
of level $\alpha$ as paths that belong
$\\{\omega\in\Omega:M(t^{\prime},\omega)\in I_{\alpha}\\}$.

Our confidence bands contain the trajectories of infected and recovered
patients that generate many deaths similar to that which occurs in reality.
Therefore, in practice, $\alpha=0.1$ is usually a fair value, and $t^{\prime}$
is selected as one of the latest datums in the historical records.

###  5.6 Software details and resources

Our proposal has been implemented in several programming languages: C++,
Python, and R, although the results shown in this article have been obtained
with Python. We optimized the parameters using library ` pycma` [82], and
`Numpy` has been used for mathematical operations.

In the different statistical analyses performed, 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 can be downloaded from at
[83] and [84].

In the immediate future we are going to release the code used in this paper
for the benefit of the scientific community at
(https://github.com/covid19-modeling).

## Acknowledgements

The authors thank Oscar Hernan Madrid Padilla for their helpful suggestions.
We are also grateful to Felipe Muñoz López for his help in the computer
implementations at the early stages of this research.

This work has received financial support from the Consellería de Cultura,
Educación e Ordenación Universitaria (accreditation 2019-2022 ED431G-2019/04)
and the European Regional Development Fund (ERDF), which acknowledges the
CiTIUS-Research Center in Intelligent Technologies of the University of
Santiago de Compostela as a Research Center of the Galician University System.

## Competing Interests

The authors declare no competing interests.

## References

  * [1] Ilona Kickbusch and Gabriel Leung.  Response to the emerging novel coronavirus outbreak.  BMJ, 368, 2020.
  * [2] Scott P. Layne, James M. Hyman, David M. Morens, and Jeffery K. Taubenberger.  New coronavirus outbreak: Framing questions for pandemic prevention.  Science Translational Medicine, 12(534), 2020.
  * [3] Martin Enserink and Kai Kupferschmidt.  Mathematics of life and death: How disease models shape national shutdowns and other pandemic policies.  Science Magazine, 2020.
  * [4] Stephen M Kissler, Christine Tedijanto, Edward Goldstein, Yonatan H Grad, and Marc Lipsitch.  Projecting the transmission dynamics of sars-cov-2 through the postpandemic period.  Science, 2020.
  * [5] A Huppert and G. Katriel.  Mathematical modelling and prediction in infectious disease epidemiology.  Clinical Microbiology and Infection, 19(11):999–1005, Nov 2013\.
  * [6] Matt J Keeling and Pejman Rohani.  Modeling infectious diseases in humans and animals.  Princeton University Press, 2011.
  * [7] Sandip Mandal, Ram Rup Sarkar, and Somdatta Sinha.  Mathematical models of malaria-a review.  Malaria journal, 10(1):202, 2011.
  * [8] Linda JS Allen.  Some discrete-time si, sir, and sis epidemic models.  Mathematical biosciences, 124(1):83–105, 1994.
  * [9] Junkichi Satsuma, R Willox, A Ramani, B Grammaticos, and AS Carstea.  Extending the sir epidemic model.  Physica A: Statistical Mechanics and its Applications, 336(3-4):369–375, 2004.
  * [10] Frank Ball, Catherine Larédo, David Sirl, and Viet Chi Tran.  Stochastic Epidemic Models with Inference, volume 2255.  Springer Nature, 2019.
  * [11] Boseung Choi and Grzegorz A Rempala.  Inference for discretely observed stochastic kinetic networks with applications to epidemic modeling.  Biostatistics, 13(1):153–165, 2012.
  * [12] Dave Osthus, Kyle S Hickmann, Petruţa C Caragea, Dave Higdon, and Sara Y Del Valle.  Forecasting seasonal influenza with a state-space sir model.  The annals of applied statistics, 11(1):202, 2017.
  * [13] J. S. Koopman and J. W. Lynch.  Individual causal models and population system models in epidemiology.  American journal of public health, 89(8):1170–1174, Aug 1999.  10432901[pmid].
  * [14] J. Gomez-Gardees, L. Lotero, S. N. Taraskin, and F. J. Perez-Reche.  Explosive contagion in networks.  Scientific Reports, 6(1):19767, 2016.
  * [15] Matt J. Keeling and Ken T. D. Eames.  Networks and epidemic models.  Journal of the Royal Society, Interface, 2(4):295–307, Sep 2005\.  16849187[pmid].
  * [16] Michael R Kosorok and Eric B Laber.  Precision medicine.  Annual review of statistics and its application, 6:263–286, 2019\.
  * [17] Nicholas G. Reich, Logan C. Brooks, Spencer J. Fox, Sasikiran Kandula, Craig J. McGowan, Evan Moore, Dave Osthus, Evan L. Ray, Abhinav Tushar, Teresa K. Yamana, Matthew Biggerstaff, Michael A. Johansson, Roni Rosenfeld, and Jeffrey Shaman.  A collaborative multiyear, multimodel assessment of seasonal influenza forecasting in the united states.  Proceedings of the National Academy of Sciences, 116(8):3146–3154, 2019.
  * [18] Logan C. Brooks, David C. Farrow, Sangwon Hyun, Ryan J. Tibshirani, and Roni Rosenfeld.  Nonmechanistic forecasts of seasonal influenza with iterative one-week-ahead distributions.  PLOS Computational Biology, 14(6):1–29, 06 2018.
  * [19] Cécile Viboud and Alessandro Vespignani.  The future of influenza forecasts.  Proceedings of the National Academy of Sciences, 116(8):2802–2804, 2019.
  * [20] Xiao Li, Jessilyn Dunn, Denis Salins, Gao Zhou, Wenyu Zhou, Sophia Miryam Schussler-Fiorenza Rose, Dalia Perelman, Elizabeth Colbert, Ryan Runge, Shannon Rego, Ria Sonecha, Somalee Datta, Tracey McLaughlin, and Michael P. Snyder.  Digital health: Tracking physiomes and activity using wearable biosensors reveals useful health-related information.  PLoS biology, 15(1):e2001402–e2001402, Jan 2017.  28081144[pmid].
  * [21] Sander Greenland.  Multiple-bias modelling for analysis of observational data.  Journal of the Royal Statistical Society: Series A (Statistics in Society), 168(2):267–306, 2005.
  * [22] Marc Lipsitch, Christl Donnelly, Christophe Fraser, Isobel Blake, Anne Cori, Ilaria Dorigatti, Neil Ferguson, Tini Garske, Harriet Mills, Steven Riley, Maria Van Kerkhove, and Miguel Hernan.  Potential biases in estimating absolute and relative case-fatality risks during outbreaks.  PLoS neglected tropical diseases, 9:e0003846, 07 2015.
  * [23] Heejung Bang and James Robins.  Doubly robust estimation in missing data and causal inference models.  Biometrics, 61:962–73, 01 2006.
  * [24] John Hopkins Coronavirus Resource Center.  https://coronavirus.jhu.edu/map.html.
  * [25] Instituto Nacional de Estadística (INE).  https://www.ine.es/.
  * [26] Population Reference Bureau.  https://www.prb.org/countries-with-the-oldest-populations/.
  * [27] AC Ghani, CA Donnelly, DR Cox, JT Griffin, C Fraser, TH Lam, LM Ho, WS Chan, RM Anderson, AJ Hedley, et al.  Methods for estimating the case fatality ratio for a novel, emerging infectious disease.  American journal of epidemiology, 162(5):479–486, 2005.
  * [28] Xianxian Zhao, Bili Zhang, Pan Li, Chaoqun Ma, Jiawei Gu, Pan Hou, Zhifu Guo, Hong Wu, and Yuan Bai.  Incidence, clinical characteristics and prognostic factor of patients with covid-19: a systematic review and meta-analysis.  medRxiv, 2020.
  * [29] Paolo Buonanno, Sergio Galletta, and Marcello Puca.  News from the covid-19 epicenter.  Available at SSRN 3567093, 2020.
  * [30] Spencer Woody, Mauricio Garcia Tec, Maytal Dahan, Kelly Gaither, Michael Lachmann, Spencer Fox, Lauren Ancel Meyers, and James G Scott.  Projections for first-wave covid-19 deaths across the us using social-distancing measures derived from mobile phones.  medRxiv, 2020.
  * [31] Keith R. Bisset, Ashwin M. Aji, Madhav V. Marathe, and Wu-chun Feng.  High-performance biocomputing for simulating the spread of contagion over large contact networks.  BMC genomics, 13 Suppl 2(Suppl 2):S3–S3, Apr 2012.  22537298[pmid].
  * [32] Dennis L. Chao, M. Elizabeth Halloran, Valerie J. Obenchain, and Ira M Longini, Jr.  Flute, a publicly available stochastic influenza epidemic simulation model.  PLOS Computational Biology, 6(1):1–8, 01 2010.
  * [33] Keith Bisset, Jiangzhuo Chen, Xizhou Feng, Sritesh Kumar, and Madhav Marathe.  Epifast: A fast algorithm for large scale realistic epidemic simulations on distributed memory systems.  pages 430–439, 01 2009.
  * [34] Y Ma, K. R. Bisset, J. Chen, S. Deodhar, and M. V. Marathe.  Formal specification and experimental analysis of an interactive epidemic simulation framework.  In 2011 IEEE International Conference on High Performance Computing and Communications, pages 790–795, Sep. 2011.
  * [35] Liang Mao and Ling Bian.  Spatial-temporal transmission of influenza and its health risks in an urbanized area.  Computers, Environment and Urban Systems, 34(3):204–215, May 2010\.  PMC7124201[pmcid].
  * [36] Centers for Disease Control and Prevention.  https://www.cdc.gov/flu/about/keyfacts.htm.
  * [37] Ottar N Bjørnstad, Bärbel F Finkenstädt, and Bryan T Grenfell.  Dynamics of measles epidemics: estimating scaling of transmission rates using a time series sir model.  Ecological monographs, 72(2):169–184, 2002.
  * [38] Lars Lorch, William Trouleau, Stratis Tsirtsis, Aron Szanto, Bernhard Schölkopf, and Manuel Gomez-Rodriguez.  A spatiotemporal epidemic model to quantify the effects of contact tracing, testing, and containment.  arXiv preprint arXiv:2004.07641, 2020.
  * [39] Juanjuan Zhang, Maria Litvinova, Wei Wang, Yan Wang, Xiaowei Deng, Xinghui Chen, Mei Li, Wen Zheng, Lan Yi, Xinhua Chen, et al.  Evolving epidemiology and transmission dynamics of coronavirus disease 2019 outside hubei province, china: a descriptive and modelling study.  The Lancet Infectious Diseases, 2020.
  * [40] Lei Luo, Dan Liu, Xin-long Liao, Xian-bo Wu, Qin-long Jing, Jia-zhen Zheng, Fang-hua Liu, Shi-gui Yang, Bi Bi, Zhi-hao Li, Jian-ping Liu, Wei-qi Song, Wei Zhu, Zheng-he Wang, Xi-ru Zhang, Pei-liang Chen, Hua-min Liu, Xin Cheng, Miao-chun Cai, Qing-mei Huang, Pei Yang, Xin-fen Yang, Zhi-gang Huang, Jin-ling Tang, Yu Ma, and Chen Mao.  Modes of contact and risk of transmission in covid-19 among close contacts.  medRxiv, 2020.
  * [41] Yan Bai, Lingsheng Yao, Tao Wei, Fei Tian, Dong-Yan Jin, Lijuan Chen, and Meiyun Wang.  Presumed Asymptomatic Carrier Transmission of COVID-19.  JAMA, 02 2020.
  * [42] Eran Segal, Feng Zhang, Xihong Lin, Gary King, Ophir Shalem, Smadar Shilo, William E. Allen, Yonatan H. Grad, Casey S. Greene, Faisal Alquaddoomi, Simon Anders, Ran Balicer, Tal Bauman, Ximena Bonilla, Gisel Booman, Andrew T. Chan, Ori Ori Cohen, Silvano Coletti, Natalie Davidson, Yuval Dor, David A. Drew, Olivier Elemento, Georgina Evans, Phil Ewels, Joshua Gale, Amir Gavrieli, Benjamin Geiger, Iman Hajirasouliha, Roman Jerala, Andre Kahles, Olli Kallioniemi, Ayya Keshet, Gregory Landua, Tomer Meir, Aline Muller, Long H. Nguyen, Matej Oresic, Svetlana Ovchinnikova, Hedi Peterson, Jay Rajagopal, Gunnar Ratsch, Hagai Rossman, Johan Rung, Andrea Sboner, Alexandros Sigaras, Tim Spector, Ron Steinherz, Irene Stevens, Jaak Vilo, Paul Wilmes, and .  Building an international consortium for tracking coronavirus health status.  medRxiv, 2020.
  * [43] Hagai Rossman, Ayya Keshet, Smadar Shilo, Amir Gavrieli, Tal Bauman, Ori Cohen, Esti Shelly, Ran Balicer, Benjamin Geiger, Yuval Dor, and Eran Segal.  A framework for identifying regional outbreak and spread of covid-19 from one-minute population-wide surveys.  Nature Medicine, 2020.
  * [44] Finbarr O’Sullivan.  A statistical perspective on ill-posed inverse problems.  Statistical science, pages 502–518, 1986.
  * [45] Martin Hanke and Otmar Scherzer.  Inverse problems light: Numerical differentiation.  The American Mathematical Monthly, 108(6):512–521, 2001.
  * [46] Peter J Bickel, Bo Li, Alexandre B Tsybakov, Sara A van de Geer, Bin Yu, Teófilo Valdés, Carlos Rivero, Jianqing Fan, and Aad van der Vaart.  Regularization in statistics.  Test, 15(2):271–344, 2006.
  * [47] Vladimir Vapnik.  The nature of statistical learning theory.  Springer science & business media, 2013.
  * [48] Ernesto De Vito, Lorenzo Rosasco, Andrea Caponnetto, Umberto De Giovannini, and Francesca Odone.  Learning from examples as an inverse problem.  Journal of Machine Learning Research, 6(May):883–904, 2005.
  * [49] Vladimir Vapnik and Rauf Izmailov.  V-matrix method of solving statistical inference problems.  Journal of Machine Learning Research, 16(51):1683–1730, 2015.
  * [50] Huili Xiang and Bin Liu.  Solving the inverse problem of an sis epidemic reaction-diffusion model by optimal control methods.  Computers & Mathematics with Applications, 70(5):805–819, 2015\.
  * [51] Alexandra Smirnova and Gerardo Chowell.  A primer on stable parameter estimation and forecasting in epidemiology by a problem-oriented regularized least squares algorithm.  Infectious Disease Modelling, 2(2):268–275, May 2017.  29928741[pmid].
  * [52] Anthony S. Fauci, H. Clifford Lane, and Robert R. Redfield.  Covid-19 - navigating the uncharted.  New England Journal of Medicine, 382(13):1268–1269, 2020.
  * [53] Michael Day.  Covid-19: four fifths of cases are asymptomatic, china figures indicate.  BMJ, 369, 2020.
  * [54] Satyajith Amaran, Nikolaos V. Sahinidis, Bikram Sharda, and Scott J. Bury.  Simulation optimization: a review of algorithms and applications.  4OR, 12(4):301–333, 2014.
  * [55] Satyajith Amaran, Nikolaos V. Sahinidis, Bikram Sharda, and Scott J. Bury.  Simulation optimization: a review of algorithms and applications.  Annals of Operations Research, 240(1):351–380, 2016.
  * [56] Corey M Peak, Rebecca Kahn, Yonatan H Grad, Lauren M Childs, Ruoran Li, Marc Lipsitch, and Caroline O Buckee.  Modeling the comparative impact of individual quarantine vs. active monitoring of contacts for the mitigation of covid-19.  medRxiv, 2020.
  * [57] Stephen M Kissler, Christine Tedijanto, Marc Lipsitch, and Yonatan Grad.  Social distancing strategies for curbing the covid-19 epidemic.  medRxiv, 2020.
  * [58] Martin McKee and David Stuckler.  If the world fails to protect the economy, covid-19 will damage health not just now but also in the future.  Nature Medicine, 2020.
  * [59] Ziyue Liu and Wensheng Guo.  Government responses matter: Predicting covid-19 cases in us under an empirical bayesian time series framework.  medRxiv, 2020.
  * [60] Bobak Shahriari, Kevin Swersky, Ziyu Wang, Ryan P Adams, and Nando De Freitas.  Taking the human out of the loop: A review of bayesian optimization.  Proceedings of the IEEE, 104(1):148–175, 2015.
  * [61] Galit Shmueli, Thomas P. Minka, Joseph B. Kadane, Sharad Borle, and Peter Boatwright.  A useful distribution for fitting discrete data: revival of the conway-maxwell-poisson distribution.  Journal of the Royal Statistical Society: Series C (Applied Statistics), 54(1):127–142, 2005.
  * [62] Thomas House.  Modelling epidemics on networks.  Contemporary Physics, 53(3):213–225, 2012.
  * [63] ZHAN Zhang and TN Krishnamurti.  Ensemble forecasting of hurricane tracks.  Bulletin of the American Meteorological Society, 78(12):2785–2796, 1997.
  * [64] Stephen A Lauer, Kyra H Grantz, Qifang Bi, Forrest K Jones, Qulu Zheng, Hannah R Meredith, Andrew S Azman, Nicholas G Reich, and Justin Lessler.  The incubation period of coronavirus disease 2019 (covid-19) from publicly reported confirmed cases: estimation and application.  Annals of internal medicine, 2020.
  * [65] AG Abdel-Salam and M Mollazehi.  Modeling survival time to recovery from covid-19: A case study on singapore.  2020\.
  * [66] Coronavirus Pneumonia Emergency Response Epidemiology Novel et al.  The epidemiological characteristics of an outbreak of 2019 novel coronavirus diseases (covid-19) in china.  Zhonghua liu xing bing xue za zhi= Zhonghua liuxingbingxue zazhi, 41(2):145, 2020.
  * [67] Robert Verity, Lucy C. Okell, Ilaria Dorigatti, Peter Winskill, Charles Whittaker, Natsuko Imai, Gina Cuomo-Dannenburg, Hayley Thompson, Patrick G. T. Walker, Han Fu, Amy Dighe, Jamie T. Griffin, Marc Baguelin, Sangeeta Bhatia, Adhiratha Boonyasiri, Anne Cori, Zulma Cucunubá, Rich FitzJohn, Katy Gaythorpe, Will Green, Arran Hamlet, Wes Hinsley, Daniel Laydon, Gemma Nedjati-Gilani, Steven Riley, Sabine van Elsland, Erik Volz, Haowei Wang, Yuanrong Wang, Xiaoyue Xi, Christl A. Donnelly, Azra C. Ghani, and Neil M. Ferguson.  Estimates of the severity of coronavirus disease 2019: a model-based analysis.  The Lancet Infectious Diseases, 2020/04/19 XXXX.
  * [68] Henrik Salje, Cécile Tran Kiem, Noémie Lefrancq, Noémie Courtejoie, Paolo Bosetti, Juliette Paireau, Alessio Andronico, Nathanaël Hoze, Jehanne Richet, Claire-Lise Dubost, et al.  Estimating the burden of sars-cov-2 in france.  2020\.
  * [69] Qifang Bi, Yongsheng Wu, Shujiang Mei, Chenfei Ye, Xuan Zou, Zhen Zhang, Xiaojian Liu, Lan Wei, Shaun A Truelove, Tong Zhang, et al.  Epidemiology and transmission of covid-19 in shenzhen china: Analysis of 391 cases and 1,286 of their close contacts.  MedRxiv, 2020.
  * [70] Katrin Zwirglmaier Ehmann, Christian Drosten, Clemens Wendtner, MD Zange, Patrick Vollmar, DVM Rosina Ehmann, Katrin Zwirglmaier, MD Guggemos, Michael Seilmaier, Daniela Niemeyer, et al.  Virological assessment of hospitalized cases of coronavirus disease 2019\.
  * [71] Hiroshi Nishiura, Tetsuro Kobayashi, Takeshi Miyama, Ayako Suzuki, Sungmok Jung, Katsuma Hayashi, Ryo Kinoshita, Yichi Yang, Baoyin Yuan, Andrei R. Akhmetzhanov, and Natalie M Linton.  Estimation of the asymptomatic ratio of novel coronavirus infections (covid-19).  medRxiv, 2020.
  * [72] Sakiko Tabata, Kazuo Imai, Shuichi Kawano, Mayu Ikeda, Tatsuya Kodama, Kazuyasu Miyoshi, Hirofumi Obinata, Satoshi Mimura, Tsutomu Kodera, Manabu Kitagaki, Michiya Sato, Satoshi Suzuki, Toshimitsu Ito, Yasuhide Uwabe, and Kaku Tamura.  Non-severe vs severe symptomatic covid-19: 104 cases from the outbreak on the cruise ship “diamond princess” in japan.  medRxiv, 2020.
  * [73] Kenji Mizumoto, Katsushi Kagaya, Alexander Zarebski, and Gerardo Chowell.  Estimating the asymptomatic proportion of coronavirus disease 2019 (covid-19) cases on board the diamond princess cruise ship, yokohama, japan, 2020\.  Eurosurveillance, 25(10), 2020.
  * [74] Dimple D Rajgor, Meng Har Lee, Sophia Archuleta, Natasha Bagdasarian, and Swee Chye Quek.  The many estimates of the covid-19 case fatality rate.  The Lancet Infectious Diseases, 2020.
  * [75] Robert Verity, Lucy C. Okell, Ilaria Dorigatti, Peter Winskill, Charles Whittaker, Natsuko Imai, Gina Cuomo-Dannenburg, Hayley Thompson, Patrick G. T. Walker, Han Fu, Amy Dighe, Jamie T. Griffin, Marc Baguelin, Sangeeta Bhatia, Adhiratha Boonyasiri, Anne Cori, Zulma Cucunubá, Rich FitzJohn, Katy Gaythorpe, Will Green, Arran Hamlet, Wes Hinsley, Daniel Laydon, Gemma Nedjati-Gilani, Steven Riley, Sabine van Elsland, Erik Volz, Haowei Wang, Yuanrong Wang, Xiaoyue Xi, Christl A. Donnelly, Azra C. Ghani, and Neil M. Ferguson.  Estimates of the severity of coronavirus disease 2019: a model-based analysis.  The Lancet Infectious Diseases, 2020/04/06 XXXX.
  * [76] Joseph T. Wu, Kathy Leung, Mary Bushman, Nishant Kishore, Rene Niehus, Pablo M. de Salazar, Benjamin J. Cowling, Marc Lipsitch, and Gabriel M. Leung.  Estimating clinical severity of covid-19 from the transmission dynamics in wuhan, china.  Nature Medicine, 2020.
  * [77] Elisabeth Mahase.  Covid-19: death rate is 0.66% and increases with age, study estimates.  BMJ, 369:m1327, Apr 2020.
  * [78] Christian Dudel, Tim Riffe, Enrique Acosta, Alyson A van Raalte, and Mikko Myrskyla.  Monitoring trends and differences in covid-19 case fatality rates using decomposition methods: Contributions of age structure and age-specific fatality.  medRxiv, 2020.
  * [79] Nikolaus Hansen.  The cma evolution strategy: A tutorial.  arXiv preprint arXiv:1604.00772, 2016.
  * [80] Nikolaus Hansen, Anne Auger, Raymond Ros, Steffen Finck, and Petr Povsik.  Comparing results of 31 algorithms from the black-box optimization benchmarking bbob-2009.  In Proceedings of the 12th Annual Conference Companion on Genetic and Evolutionary Computation, GECCO 10, pages 1689–1696, New York, NY, USA, 2010. Association for Computing Machinery.
  * [81] David H Brookes, Akosua Busia, Clara Fannjiang, Kevin Murphy, and Jennifer Listgarten.  A view of estimation of distribution algorithms through the lens of expectation-maximization.  arXiv preprint arXiv:1905.10474, 2019.
  * [82] niko, Youhei Akimoto, yoshihikoueno, Dimo Brockhoff, Matthew Chan, and ARF1.  Cma-es/pycma: r3.0.3, April 2020.
  * [83] COVID-19 Rubén Fernández Casal.  https://github.com/rubenfcasal/COVID-19.
  * [84] Datadista.  https://github.com/datadista/datasets/tree/master/COVID%2019.