Outbreaks of infectious diseases have been a common occurrence throughout history, often linked or followed by disruptions in societies and human activities 1 . There are several ways to measure the impact of outbreaks but death tolls are the most relevant ones whenever the disease can threaten lives. For instance, 32 million persons have died between 1981 and 2019 in the ongoing HIV epidemics, 700 000 in 2018 alone 2 . Aside from this large scale epidemic, the world has experienced several other recent outbreaks with varying degrees of severity and scale: Zika fever, whose symptoms are mild but can produce long-lasting effects in newborns (microcephaly) 3, 4 ;
Ebola virus disease, with a high mortality rate estimated between 20 and 75% 5, 6 ; Swine flu/H1N1, which became a pandemic in 2009-2010 although with a lower mortality rate than regular flu 7 .
In 2019-2020, the severe acute respiratory syndrome COVID-19 has emerged as the most recent pandemic, caused by the virus named SARS-CoV-2 8, 9 . Due to its novelty and lack of previous exposition, humans have no immunity against this threat, leading to increased number of infections.
At the time of this writing, the specifics of the pathogen transmission are still being investigated, as well as the complete infection process once the virus enters the host. However, it has been shown that the main human-to-human transmission mode occurs by the spreading of contaminated droplets, similar to other flu-like diseases 8 . In sharp contrast with H1N1, however, the mortality rate of COVID-19 is estimated in the range 1-4%, with higher ill-outcomes among persons of old age 10, 11 . This situation creates a unique scenario where healthcare facilities and workers can be overwhelmed in a short period of time, ultimately leading to untreated patients of COVID-19 as well as other diseases 12 .
To make matters even worse, asymptomatic patients can spread the pathogen for an extended period of time, showing none or mild symptoms during the course of the infection. As a result, laboratory tests to detect the viral load are necessary to identify the correct number of cases outside hospital environments. Similar to the H1N1 pandemic, the required number of tests far exceeds the current amount available in most countries. Without timely tracking of new cases, contact tracing becomes a challenging task, hindering estimates of new cases per infection summarized by the basic reproduction number, R 0 . The significance of this parameter lies in the fact that it provides a way to glimpse the values of the transmission rates. Those can then be used in compartmental models -mathematical models that describe the evolution of epidemics assuming nearly homogeneous populations 13, 14 . Earlier estimates for R 0 using epidemiological data from Wuhan, China, set R 0 between 1.5 to 5.7 without additional measures to restrict the spreading [15] [16] [17] . With measures in place -such as lockdown, self-isolation, and social distancing -R 0 was estimated to be around 1.05 15 . More importantly, the insufficient number of tests, in addition to the long waiting time for lab results, affects the quality of epidemic models.
In the absence of mass testings, the death toll can be used as an alternative metric to probe the extension of the epidemic. The medical staff can assess the cause of death from clinical reports, which may contain test results or not, using the best of their knowledge. Additional tests may also be appended to reports to further specify the cause of death.
Both numbers of cases and deaths are publicly available as part of a global effort to tackle the pandemic. Here, we study the evolution of COVID-19 deaths in order to reduce the issues caused by the limited number of laboratory confirmed tests. We show that the accumulated deaths can be effectively described by simple functions, namely, sigmoids whose parameters are explained in terms of the SIR epidemic model. The SIR model is a compartmental model with the following health states: susceptible, infective, and removed 18 . The removed state represents those who have passed away or have developed immunity, either by recovering from the disease or any other method such as vaccination. The SIR model was chosen in detriment of epidemic models with additional health states or reinfections because it is the simplest model that addresses immunity.
Among our results, we show that crude mortality rates can be computed from the parameters of sigmoids and that the rates can change up to one order of magnitude, depending on the severity of the outbreak in a given region. The paper is organized as follows. Sec. 2 contains the description 2/20 of the data and variables used along the text. Sec. 3 explains how the SIR model is reduced from a system of differential equations to a single non-linear differential equation, with emphasis on the expansion around equilibrium. Data is modeled in Sec. 4 via sigmoidal functions, whose parameters are explained in terms of the epidemiological parameters of the SIR model. Time windows are addressed in Sec. 5, with special emphasis on R 0 , crude mortality rate, and quantitative effects of the confinement. Final comments and conclusions are listed in Sec 6.
The European Center for Disease Prevention and Control (ECDC) provides COVID-19 data updated in a daily schedule 19 . The daily reports portray the distribution of new confirmed cases and new deaths presented as time series. The dataset also displays the population size N according to 2018 World Bank census for each geographical region (see Table 1 ).
To better grasp the nature of the data, consider the number of accumulated deaths in France, as shown in Fig. 1 . Similar to other European countries, France was heavily affected by the pandemic nearing the mark of 30 000 deaths, with a sharp increase in deaths of infected patients around April-May. These deaths are linked to infections which took place between 3 to 6 weeks prior. Before diving into modeling, we seek for general features that can be used to model outbreaks.
These features concern the spreading regime of the disease among the population, excluding spatial effects and temporal delays. As noted above, the death toll experienced a rapid growth between March and April (see Fig. 1 ). This regime is the hallmark of epidemics and denotes the exponential phase. In general, the growth rate in compartmental epidemic models is summarized by R 0 that is After the exponential phase, the system relaxes toward an equilibrium state with no infective people, with a characteristic time scale, namely, the relaxation time τ . We shall investigate the relationship between τ with epidemiological parameters in the next section. Considering these three aspects, and using Fig. 1 as reference, it becomes clear that the number of accumulated deaths is a monotonic function, with an early exponential growth only to be replaced by a smooth relaxation towards equilibrium. Therefore, sigmoids appears as ideal candidates to model the data.
For the sake of simplicity, let us assume the spreading dynamics can be approximated by the standard SIR model in a homogeneous population of size N. The model comprises a population whose subjects can be classified in three distinct health states, namely, susceptible, infective, and removed. The removed state includes individuals that have either died or recovered from 4/20 the disease. The latter are assumed to not be infective anymore, nor susceptible to become sick again because of some kind of immunity. The fraction of individuals in each compartment is, respectively, S(t), I(t), and R(t), at time instant t. The dynamic goes as follows. Infective subjects in the population transmit the pathogen to susceptible ones, under adequate conditions.
The transmission occurs with rate α, and we assume the homogeneous mixture of the population, that is, each person in the population is statistically equivalent to each other 14 . Once infected, the person remains infective for an average period 1/λ, where λ is the removal rate. The dynamical equations that describe the model are: can be shown 24 that the equation can be further reduced to
whose general solution can be obtained by quadrature
The stationary condition in (2) gives the value R ∞ as a solution of the transcendental equation
The solution R(t) can be obtained by inverting (3), for which a general method remains unknown. To circumvent the issue, one may expand (2) The curve that represents the infective fraction of the population in the SIR model is symmetrical, with peak at the center, as long as parameters remain constant.
around R = 0 or R = R ∞ . Each expansion has advantages and issues. The linear term in the expansion around R = 0 dictates the exponential growth of R(t), with an effective rate αS 0 − λ.
We use the adjective effective rather loosely here because at some point the curve should change its curvature and converge to an equilibrium value. Unfortunately, the competition between the early exponential growth and relaxation towards equilibrium is often difficult to assess near the onset of the outbreak, requiring higher order contributions in the expansion.
Alternatively, we can get a better picture of the problem by expanding R(t) around the equilibrium. By doing so, the expansion only requires contributions up to second order, as it already carries the information regarding R ∞ . Furthermore, the expansion near the equilibrium also ensures that τ describes the dominant relaxation time, rather than combinations of several decay modes, each one with its own timescale. Define
Keeping terms up to o(δR 2 ), one converts (4) into a Bernoulli equation with relaxation time
Solving for δR and transforming back to R, we find the approximate solution (see Fig. 2 )
The sigmoidal expression in (5) satisfies the requirements listed in Sec. 2 and it is an excellent candidate to model the data of accumulated deaths. However, the removed compartment corresponding to R(t) holds both recovered and deceased fractions of the population, i.e., all the infected who are unable to spread the disease. We can simplify this issue by assuming the existence of a simple relation between R(t) and the accumulated number of deaths divided by the population size, g(t),
at the time instant t. To keep the model as simple as possible, we neglect temporal delays and impose
where the crude mortality rate f is the ratio between deceased and infected. As such, it also can be used as an estimator for the likelihood to die after contracting the disease.
The equilibrium value g ∞ varies from country to country and can be used to characterize the outbreak, more specifically, to assess the impact of the outbreak in the afflicted population. The monotonic nature of cumulative quantities together with the upper bound g ∞ restrict the possible functional forms for g(t). Sigmoids are natural candidates to describe g(t) since they are bounded and monotonic. Here we consider the following expression in tandem with (5):
The timescale τ dictates the relaxation of g(t) towards g ∞ , with inflection at t c = τ ln b. Equation (7) has 3 or 4 parameters depending on whether g eff (t) must pass by the initial entry g eff (0) = g(0) = g 0 or not. The former implies the constraint a = g ∞ − g 0 (1 + b), whereas the later g eff (0) = (g ∞ − a)/(1 + b) with a g ∞ . Since the data are noisy, we do not require the fitting curve to pass by g 0 .
The relationship between g(t) and R(t) is consistent if the following equations hold: the beginning of the outbreak is unknown, the assumption S 0 ≈ 1 is no longer valid. Alternatively, one may consider an estimate for the removal rate λ, or its probability distribution. In this case, the input parameter to solve (8a-8c) depends solely on the characteristics of the disease and local demographics. It does not require details concerning the disease spreading nor the moment in which the outbreak begins. Therefore, evidence-based values for λ from patient data are far more suitable for the purposes of this study, and shall be used hereafter. is held fixed during the optimization (not shown), then the optimal set of parameters become highly sensitive to the initial guess of the optimization algorithm.
The observation above inquires whether the SIR model describes the data or not. However, this issue can be understood by inspecting the number of daily deaths (see Fig. 1a ). An asymmetry is observed around the peak of daily deaths (see Fig. 1b ), which is in sharp contrast with the corresponding curve in the theoretical SIR model (see Fig. 2b ). The asymmetry can be explained by a variation in the transmission rate. Unlike other recent epidemics, several countries have adopted control measures such as lockdowns of non-essential workers, flight and other travel restrictions.
The efforts effectively reduce the transmission rate of COVID-19 once they are in place. Thus, the data must be divided in non-overlapping time windows, each with its own set of epidemiological parameters.
As a first example of data where the necessity to consider several time windows is obvious, we consider the global deaths between February 29 and May 31 (see Fig. 4 ). This case is interesting while increasing the number of data points. r 2 forms a plateau for sizes between 29 to 41 days followed by a sudden increase for 42 days, indicating that 41 days is the optimal size of the time window (February 29 -April 9).
In the following, we decide to divide the observation time in two consecutive time windows.
Let us explain how we determine the optimal separation time (optimal duration of the first time We can now move to more complicated cases, in which the data themselves contains artifacts and the fitting procedure can be tricky. That is the case of France, for instance. The lockdown was where the notation R 0,k with k = 1 and 2 refers to the first and second time windows, respectively.
In addition, the crude mortality rate in the second time window reads f 2 = (4.67 ± 1.16) × 10 −3 , one order of magnitude higher than Denmark or Germany (see Table. 2).
The analysis in the first time window requires more care. The fitting procedure requires bounds, otherwise the fit may incorrectly converge to an equilibrium value that is much lower than the one obtained for the second time window. In addition, multiple solutions can be found for (8a- 
The tragic developments in the COVID-19 pandemic have exposed flawed aspects in protocols used to assess large scale epidemics. Despite the various improvements in the global capacity to produce laboratory tests, usually based on reverse-transcription polymerase chain reaction (rRT-PCR), the majority of afflicted countries were unable to enforce mass testing policies. This shortcoming has also been experienced in 2002 with the SARS epidemic but in a smaller scale. The lack of mass 15/20 testing contributed in keeping the number of cases unknown, affecting the accuracy of disease spreading models. In this paper, we resort to the death toll from April 4 to May 31 instead of reported cases as our primary data source as death certificates are mandatory and may contain other medical assessments linking the cause of death with COVID-19 outbreaks. We emphasize this methodology is not immune itself from sub-notification, nor unforeseen delays on death certificates.
We model the death toll via the sigmoid curve in (7), which requires the parameters g ∞ (capacity) and τ (time scale). The capacity limits the growth of the outbreak, which is expected to vary with geographic region, healthcare quality, and efforts to control the spreading of the virus. Together with the crude mortality rate f , which also depends on healthcare and population demographics, they give access to a far more credible estimate for the infected fraction of the population R ∞ = g ∞ /f . To keep the model as simple as possible, we also neglect temporal delays between g(t) and R(t). The inclusion of said delay may introduce additional effects, but those are not expected to be dominant in this case since spatial effects are not being investigated.
The advantage of our approach relies on the curve fitting of a monotonic curve (death toll) to a sigmoid in sharp contrast to complex optimization schemes for models with multiple health states 20, 21 . Epidemiological parameters are extracted by solving the algebraic system of equations (8a-8c) for transmission rate, crude mortality rate and initial condition, respectively, α, f, and Concerning deviations in parameters estimated through (8a-8c), they can be tracked to the uncertainties in the removal rate. In general, the inverse of the removal rate describes the average time required for an infective person to change to the removed state. For COVID-19, virus shedding occurs more prominently between 3 and 7 days, with a peak at day 5 since the onset of symptoms 25 .
The exact probability distribution for λ remains an open issue so our analysis assume a uniform distribution.
Finally, the SIR model rests on the random mixing hypothesis but deviations are expected with stronger effects in population with varying demographics or populations at risk. In particular, communicable respiratory diseases become major issues in correctional facilities, given the lack of adequate environmental and sanitary conditions. The combination of higher transmission rate of pathogens and reduced removal rate with reduced healthcare can increase the crude mortality rate for incarcerated individuals. By a similar argument, outbreaks in nursing homes may affect estimates of epidemiological parameters because the disease becomes disproportionately more lethal for older patients. In this study, both effects are neglected in hopes to understand the disease spreading of the average population with the minimal number of parameters possible.