key: cord-0914119-v3b1jejy
authors: Tiwari, Sankalp; Vyasarayani, C. P.; Chatterjee, Anindya
title: Data suggest COVID-19 affected numbers greatly exceeded detected numbers, in four European countries, as per a delayed SEIQR model
date: 2021-04-14
journal: Sci Rep
DOI: 10.1038/s41598-021-87630-z
sha: bfea42e16b8ba874a12578b9dffef016b56834aa
doc_id: 914119
cord_uid: v3b1jejy

People in many countries are now infected with COVID-19. By now, it is clear that the number of people infected is much greater than the number of reported cases. To estimate the infected but undetected/unreported cases using a mathematical model, we can use a parameter called the probability of quarantining an infected individual. This parameter exists in the time-delayed SEIQR model (Scientific Reports, article number: 3505). Here, two limiting cases of a network of such models are used to estimate the undetected population. The first limit corresponds to the network collapsing onto a single node and is referred to as the mean-[Formula: see text] model. In the second case, the number of nodes in the network is infinite and results in a continuum model wherein the infectivity is statistically distributed. We use a generalized Pareto distribution to model the infectivity. This distribution has a fat tail and models the presence of super-spreaders that contribute to the disease progression. While both models capture the detected numbers well, the predictions of affected numbers from the continuum model are more realistic. Our results suggest that affected people outnumber detected people by one to two orders of magnitude in Spain, the UK, Italy, and Germany. Our results are consistent with corresponding trends obtained from published serological studies in Spain, the UK and Italy. The match with limited studies in Germany is poor, possibly because Germany’s partial lockdown approach requires different modeling.

www.nature.com/scientificreports/ These models are obtained by considering two limiting cases of a time-delayed network SEIQR model motivated by the model of Young et al. [13] [14] [15] . In the network model, the whole population is divided into N subpopulations based on their net infectivity ( β ) values, and each node represents a sub-population or group. In the first limiting case that we adopt, which is the same as a mean-β model 15 , the entire network 14 is collapsed into a single node ( N = 1 ). This model was originally proposed by Young et al. 13 , and for a fast pandemic some simplifications and approximations are possible 15 . In the second limiting case 14 , which is a continuum model, we take N → ∞ . Here, the infectivity is treated as a continuously distributed parameter in the population. Note that these models are not network models themselves, but are the extreme limits of an underlying network model. Please see the books by Barrat et al. 16 and Barabási 17 for more details on network models.

For Italy, Germany, the UK, and Spain, we fit these two models to the data reported under the heading 'total cases' on the Worldometer website 18 . We consider the data from February 15-June 18 for fitting (125 days). Beyond mid-June, all the countries seemed to be experiencing a second wave of COVID-19 after relaxing social distancing norms, or perhaps due to increased testing rates. Therefore, the constancy of model parameters beyond mid-June may not be a reasonable assumption.

Both these models include a parameter called the probability of detecting an infected individual. This parameter, upon fitting from detected population data, allows us to indirectly estimate the affected but undetected population. We will find that the continuum model fits the data better than the mean-β model. The continuum model predicts that the affected people outnumber the detected people by 8, 22, 48 , and 130 times in Spain, the UK, Italy, and Germany, respectively. However, it is emphasized that only the officially detected numbers are used for data fitting. The continuum model also indirectly suggests the presence of 'super-spreaders' (or 'super-spreading events') in all the countries, in the form of a fat tail in the distribution of the infectivity β in the population.

We point out the differences between our network model and complex network models. In our model, each node k represents a portion of the susceptible population with infectivity β k . Further, each node k is connected to every other node r = 1, 2, . . . , N with coupling coefficient β kr = √ β k √ β r . In complex network models, these coupling coefficients are referred to as weights 19 . In the continuum limit of our model, when N → ∞ , the initial infectivity corresponding to each node (population group) is assumed to follow a Pareto distribution φ(β).

In complex scale-free networks, the number of nodes are large but finite. Usually, their degree distribution follows a power-law 20 . Sometimes, the distribution of weights of the network may also have a fat tail 19 . Super spreading events are known to occur in complex networks 20 . For example, in the work of Saumell-Mendiola et al. 21 , it was showed that by connecting two complex scale-free networks by a small number of connections (super spreaders), it is possible to achieve endemic equilibrium. Without this small number of connections, the disease did not progress in either of their networks. Finally, we emphasize that our network model's continuum limit allowed us to make certain mathematical simplifications 14 , which led to a collapse in the dimensionality of the system. The rest of this paper is organized as follows. In "SEIQR models" section, we briefly discuss the two models (mean-β and continuum) used in this work. "The case of Italy" section, we present and discuss in detail the results of the optimization calculations (i.e., parameter fitting) for Italy. In "The cases of Germany, UK, and Spain" section, we present the results for the remaining three countries: Germany, the UK and Spain. In "Conclusions" section, we present our conclusions. In the supplementary material for this paper, we present a detailed study investigating other infectivity distributions in the continuum model, and some sensitivity analysis results.

A detailed description of the mean-β model 15 ( N = 1 ) and the continuum model 14 ( N → ∞ ) can be found in the literature. In this section, we describe them briefly for clarity and completeness. Mean-β model. The mean-β model 15 can be derived from the five state SEIQR model of Young et al. 13 by assuming no loss of immunity after recovery. This assumption is valid for a fast pandemic like COVID-19. In the mean-β model, exposed ( E m ), quarantined ( Q m ), and recovered ( R m ) states becomes slave variables of susceptible ( S m ) and infected ( I m ) states whose dynamics are governed by the following DDEs:

The parameters p m , γ m , τ m , β m , and σ m are described in Table 1 . The subscript m in all the quantities serves to distinguish them from those used in the continuum model ( N → ∞ ). By defining and integrating Eq. The complete dynamics of the pandemic in the mean-β can be captured by the first-order nonlinear DDE given by Eq. Continuum model. The other limit of the network model 14 is for the case of N → ∞ , which implies that the infectivity ( β ) is now distributed continuously over the population. The governing differential equations for the states S and I in this case are as follows:

where p, γ , τ , and σ are described in Table 2 . It has been shown 14 that S(β, t) admits a solution of the form:

Therefore, if the initial distribution of infectivity in the population, φ(β) , is specified, the subsequent variation of S is simply through f(t). Using algebraic manipulations, it has been shown that f(t) satisfies the following non-linear DDE 14 : It is useful to introduce a new random variable u = √ β with probability density function ψ(u) , such that ψ(u) = 2uφ(u 2 ) . Then, and

The quantities corresponding to Eqs. (8) and (9) now are and In the present work, we assume ψ(u) to be of the form which is a generalized Pareto distribution 25 . This means that Note that

In the supplementary material, we have presented a detailed study investigating other distributions for ψ(u) . That discussion supports our choice of the generalized Pareto distribution.

A summary of the parameters used in the continuum model is presented in Table 2 . Note that m > 2 ensures that u has a finite mean. If m > 3 , then u has finite variance as well.

The fitting error. In the mean-β model, for a given set of parameter values, we compute h m (t) and fit it with the data for the total number of detected cases as reported on the Worldometer website 18 . This is done by minimizing the fitting error We see from Table 1 that there are four parameters to be identified in the mean-β model. The fitting error for the continuum model is defined as We see from Table 2 that there are five parameters to be identified in the continuum model.

To test the sensitivity of both the models to both fixed and fitted parameters, the parameters are varied by ±2% around the values corresponding to the optimum fit. We generate 10,000 parameter sets using the Latin Hypercube sampling command lhsdesign in MATLAB 26 . Including the two externally specified parameters σ and γ , the hypercube is six dimensional in the mean-β model, and seven dimensional in the continuum model. The fits corresponding to these samples are plotted using bands of lighter shades around the optimum fits plotted in darker shades. Results will be presented in "The case of Italy" and "Results for the mean-β model" sections below. We also study the sensitivity of the fitted parameters to the externally specified parameters σ and γ by varying them by ±2% around σ = 3 and γ = 0.07 . The results of this latter analysis are presented in the supplementary material.

The data for the detected cases match very well for Italy. In this section, we present detailed results for Italy obtained from the two models. The results for other countries will be presented in the next section.

Results for the mean-β model. We minimize the fitting error E 0m (see Eq. (20)) using the optimization routine fminsearch in MATLAB. Since there are four free parameters in the mean-β model, the input variable for the optimization code is a four-by-one column vector, suitably transformed so that the constraints in Table 1 are automatically satisfied. We have performed several hundred optimization calculations with random initial conditions and have found many converging solutions. Several of these solutions correspond to nearly identical and low values of E 0m . Several other local minima yielded significantly higher E 0m values, and were discarded.

The parameter set that yields the lowest E 0m in all the random trials is reported in the first row of Table 3 . The fit generated using these parameters, along with the reported data, is shown in the top-left panel of Fig. 1 . The reported data, which records the percentage of detected cases in Italy from February 15 for the following 125 days, is plotted in green circles. For easier visibility, only the data of alternate days is plotted. To account for the initial uncertainty in the reporting, we neglect initial data where the number of cases is less than 1% of the number reported on the 125th day. The fitted h m (t) is plotted for a longer duration using a dashed line to depict the saturation value clearly. In the figure, the percentage of detected cases saturates at 0.3966% of Italy's population.

We also plot the percentage of infected population ( w m (t) ) in the top-right panel of Fig. 1 with a dashed curve. From the mean-β model, the percentage of infected (affected) people early during the progression of the pandemic is around 9% , while the saturation value is 90.4331% . Both these numbers seem too high, as will be discussed further below. For this same model, the ratio of the population affected to population detected ( A/D = w m /h m ) saturates at 228. The basic reproduction number 13 R 0 , for the mean-β model, is found from fitted parameters to be

The mean-β model does offer some further useful insights into data fits, as follows. Upon inspection of the local minima obtained from the fminsearch runs, we noted that all the minima corresponding to low values of E 0m have nearly identical β m and V 0 (equal to the values reported in the first row of Table 3 ), but different values for p m and τ m . The fitted values of p m were consistently low, however. To investigate further, we fix the values of β m and V 0 , and plot E 0m in the p m − τ m plane, in the mid-left panel of Fig. 1 , for low values of p m . We see that the lowest values for E 0m are obtained on a thin band cutting across the p m − τ m plane, which spans the entire assumed range of τ m and a relatively much smaller range of p m . Finally, upon plotting E 0m in the p m − τ m plane in the bottom-left panel of Fig. 1 , we observe that the thin band corresponds to almost fixed value of p m ≈ 0.0048 , indicating that p m can be robustly identified, along with β m and V 0 . In contrast, τ m is essentially indeterminate. We now consider the continuum model.

Results for the continuum model. In this case, there are five free parameters to be estimated. Therefore, the input variable to the optimization code is a 5 × 1 vector, suitably transformed so that the constraints in Table 2 are automatically satisfied. The parameter set that results in the lowest value of E 0 for all the optimization trials is reported in the first row of Table 4 . The corresponding fit h(t), again plotted for a longer duration to show saturation, is plotted using a solid line in the top-left panel of Fig. 1 . The figure indicates that the fit to the detected data from the continuum model is slightly better than that from the mean-β model (numerically, E 0 = 0.6800 for the continuum model, while E m0 = 1.8770 for the mean-β model). The optimizing value of m is found to be 2.6688 and corresponds to a fat tail in φ(β) 's distribution as shown in Fig. 2 . The percentage of detected cases saturate at 0.4051% , which is only slightly more than that predicted by the mean-β model. We also note that a variation of ±2% in the parameters (both externally specified and fitted) changes this number by about 20% , as opposed to about 7% in the mean-β model. This means that the continuum model is more sensitive to small changes in parameters, and hence, the fitted parameters are somewhat more robustly determined, as compared to the mean-β model. When comparing the continuum model with the mean beta model, a large difference is seen in the estimated number of affected people. We plot w(t) in the top-right panel of Fig. 1 using a solid line. The continuum model predicts that the percentage of population infected, or affected, in Italy saturates at 19.4350% , with an affected-todetected ratio (A/D) of 48 at saturation. Also the percentage of affected population at the start (on the 18th day) is very small. These numbers are intuitively more satisfactory, and we now compare them with a serological study.

We could not find any nationwide serological study for Italy. A big serological study in Italy was carried out in the northeastern region of the country, which is one of the most affected regions in Italy, between May 5 and May 15 27 . It included approximately 6000 participants. This study reported a seroprevalence of 23.1%(95% CI 22.0 − 24.1) , while our continuum model predicts that 16.79% and 17.74% of the Italian population was infected on May 5 and May 15, respectively. Note that we have used the nation-wide data to fit the continuum model, while the serological study was carried out in a region with high incidence. It is not surprising then that our continuum model gives a lower estimate.

A small-f expansion of Eq. (13) with 2 < m < 3 , yields terms of the form f 3−m . This fractional power makes linear stability analysis impossible, but significantly improves the fit in the early stages of the pandemic, as shown in the supplementary material. Thus our fitted value of m = 2.6688 for Italy provides indirect evidence for a fat tail in the distribution of u, which corresponds to an active role of superspreaders in the progression of the pandemic.

For distributions with faster decaying tails, or with no tails at all (i.e., with finite support), the dominant term on the right hand side of Eq. (13) turns out to be O(f ) , for small f. Linearization is then possible but as explained in the supplementary material, the fit in the initial stages is visibly poorer. One consequence of m < 3 is that the idea of the R 0 is not applicable to this model. Note that as soon as the pandemic progresses even a little bit, and f takes a strictly positive value, G of Eqs. (14) and (15) has an exponentially decaying envelope, and the subsequent dynamics is better behaved. This is why in the latter stages of the pandemic, all the distributions studied in the supplementary material give equally good fits. A detailed theoretical investigation of the consequence of m < 3 is left for future work. Here we restrict ourselves to noting from numerical simulations that good fits require m < 3 , i.e., a fat tail in the distribution of u.

Upon inspecting the local minima obtained from the fminsearch runs, we found that all the minima corresponding to low values of E have nearly identical values of a, m and f 0 (reported in the first row of Table 4 ) but different values of p and τ . These observations are similar to those from the mean-β model. Moreover, the values of p are low, while the values of τ vary over its entire assumed range. For more insight, we fix a, m and f 0 , and plot E 0 in the p − τ plane (for low values of p) as shown in the mid-right panel of Fig. 1 . We see that the lowest values of E 0 are obtained on a thin band cutting across the p − τ plane, spanning the entire assumed range of τ and a relatively much smaller range of p. Similar to the estimation results of the mean-β model, τ and p remain indeterminate even for the continuum model. Upon plotting E 0 in the p − τ plane in the bottom-right panel of Fig. 1 , we observe that the thin band of minimum values corresponds to an almost fixed value of p ≈ 0.0207. Figure 2 . Plots of the fat-tailed distribution ( φ(β) ) as used in the continuum model for Italy, Germany, the UK, and Spain.

Scientific Reports | (2021) 11:8106 | https://doi.org/10.1038/s41598-021-87630-z www.nature.com/scientificreports/ In the next section, we report results for Germany, the UK, and Spain. We will see that the main features of the results reported for the case of Italy hold for these countries as well.

The best fits obtained using the two models for these three countries are presented in Fig. 3 . For Germany (see the top-left panel of Fig. 3 ), the continuum model under-predicts and for Spain (see the bottom-left panel of Fig. 3 ), the continuum model over-predicts the actual data near the end. However, for the UK (see the mid-left panel of Fig. 3) , the fit is excellent. The variation of E 0m obtained from the mean-β model in the p m − τ m plane (with β m and V 0 fixed at values reported in Table 3 ) is plotted on the left panels of Fig. 4 . The variation of E 0 obtained from the continuum model in the p − τ plane (with a, m and f 0 fixed at values reported in Table 4 ) is plotted on the right panels of Fig. 4 . The affected-to-detected (A/D) ratio corresponding to the best fits for the mean-β model and the continuum model are reported in Tables 3 and 4 , respectively.

We see that the fitting results are qualitatively similar to those obtained for Italy. However, there are a few observations that stand out in the results for the continuum model as highlighted below:

• We see in Table 4 that the optimum value of m for Germany, the UK, and Spain is around 2.5. This indicates that the infectivity distribution φ(β) for each of these countries has a fat tail as can be in Fig. 2 . We see from the figure that most of the population has small infectivity (low value of β ). However, these curves decay to zero very slowly, since the first moment of φ(β) is infinite, i.e.

Such a distribution is consistent with the presence of 'super-spreaders' (people who have, or events that lead to, a high value of β ). Our results are in line with findings of Wong et al. 28 that the distribution of infections caused by an index case is fat-tailed. The role played by different kinds of super-spreading events, and empirical evidence for their role in the COVID-19 pandemic, are discussed in detail by Althouse et al. 29 (see also references therein). See the works by Adam et al. 30 and Britton et al. 31 for mathematical models discussing the role of super-spreading in COVID-19. Super-spreading events have been widely reported for COVID-19 in the medical literature as well 32,33 . • We see from Table 4 that the affected-to-detected ratio (A/D) is high for all the countries and varies between 8 for Spain and 130 for Germany. • The ratio of symptomatic cases to detected cases can be approximated from the value of p. Symptomatic cases outnumber the reported cases by about 5 times in Spain, 12 times in the UK, 25 times in Italy, and 60 times in Germany. 34 . Our continuum model predicts this number to be 8.63% . We could not find a comparable nationwide study for Germany. Several studies are underway, but their results are not as yet published 35 . However, the match with limited studies 36, 37 in Germany is poor, possibly because Germany's partial lockdown approach requires different modeling, e.g., a model with spatial structure.

In this work, we fit the data for the total number of infected people in four western European countries. We use two limiting cases of the time-delayed network SEIQR model: the mean-β model and the continuum model. In earlier works, it was shown that for fast pandemics, each of these two models reduces to one non-linear delay differential equation.

After fixing the values of the biological parameters σ and γ , we need to identify four parameters in the meanβ model and five parameters in the continuum model. In both the cases, we see that there are many parameter sets that minimize the fitting error, yielding almost identical values of the objective function. All these sets have almost identical values of all parameters other than p and τ . Other subsidiary quantities such as the total number of infected people, the affected-to-detected ratio, and the basic reproduction number are also close to each other in value. By plotting the fitting error in the p − τ plane (with other parameters fixed at their identified values), we see that a narrow band yielding minimum error cuts across this plane, spanning the entire assumed range of τ and a small range of p comprising low values. In this light, the τ values reported in Tables 3 and 4 may not be reliable.

We see from the results that the continuum model yields superior fits in comparison to the mean-β model. Even the worst fit obtained from the continuum model, for Spain, has 2-norm fitting error of only 1.28% . Moreover, it gives reasonable and physically realizable values for all the epidemiological quantities. We conclude from the results that the predicted number of infected people depends on the assumed distribution of the infectivity rate. A single homogeneous infectivity rate overestimates the seroprevalence in the countries examined.

Both the models are relatively insensitive to small changes in the input parameters. The continuum model is more sensitive than the mean-β model, which indicates indirectly that its parameter estimates are somewhat more robust.

The most important prediction from the models is that the total number of affected people far outnumber the people detected with COVID-19 in all the four countries. The continuum model predicts that the affected-to-detected ratio, in increasing order, is 8, 22, 48, and 130 for Spain, the UK, Italy, and Germany, www.nature.com/scientificreports/ respectively. The first three of these numbers are consistent with the serological surveys conducted in the corresponding countries. In particular, our estimated number of infected people in Spain in early May was about 5% as per our model, in agreement with a nation-wide seroprevalence study 8 . We emphasize that the detailed work done in this seroprevalence study retains primary importance. Our work provides a mathematical supporting view of consistency, and is not intended to replace such detailed seroprevalence studies. Our numbers for affected people in Italy and the UK match reasonably well with seroprevalence data from these countries 27, 34 . In contrast however, our predicted numbers for Germany are too high. This may be because the partial lockdown approach of Germany requires spatial structure within the model. Table 3 . The right side shows the variation of E 0 in the p − τ plane (for low values of p ) obtained from the continuum model. The parameters a, m, and f 0 are fixed at the values reported in Table 4 .

Substantial undocumented infection facilitates the rapid dissemination of novel coronavirus (SARS-CoV-2)

Estimating the effects of non-pharmaceutical interventions on COVID-19 in Europe

State-level tracking of COVID-19 in the United States

INDSCI-SIM A State-level Epidemiological Model for India. Ongoing study at https:// indsc icov

Our World in Datahttps:// ourwo rldin data. org/ coron avirus

A demographic scaling model for estimating the total number of COVID-19 infections

How the pandemic might play out in 2021 and beyond

Prevalence of SARS-CoV-2 in Spain (ENE-COVID): a nationwide, population-based seroepidemiological study

Seroprevalence of antibodies to SARS-CoV-2 in 10 sites in the United States

Full Lockdown Policies in Western Europe Countries Have No Evident Impacts on the COVID-19 Epidemic

Consequences of delays and imperfect implementation of isolation in epidemic control

Complete dimensional collapse in the continuum limit of a delayed SEIQR network model with separable distributed infectivity

New approximations, and policy implications, from a delayed dynamic model of a fast pandemic

Dynamical Processes on Complex Networks

The architecture of complex weighted networks

Epidemic processes in complex networks

Epidemic spreading on interconnected networks

Outbreak dynamics of COVID-19 in Europe and the effect of travel restrictions

Reconciling early-outbreak estimates of the basic reproductive number and its uncertainty: framework and applications to the novel coronavirus (SARS-CoV-2) outbreak

Modelling the COVID-19 epidemic and implementation of population-wide interventions in Italy

Wiley StatsRef: Statistics Reference. Online 1-10

A network model of Italy shows that intermittent regional strategies can alleviate the COVID-19 epidemic

Prevalence of SARS-CoV-2 IgG antibodies in an area of northeastern Italy with a high incidence of COVID-19 cases: a population-based study

Evidence that coronavirus superspreading is fat-tailed

Superspreading events in the transmission dynamics of SARS-CoV-2: opportunities for interventions and control

Clustering and superspreading potential of SARS-CoV-2 infections in Hong Kong

A mathematical model reveals the influence of population heterogeneity on herd immunity to SARS-CoV-2

Unlocking UK COVID-19 policy

Serology-and PCR-based cumulative incidence of SARS-CoV-2 infection in adults in a successfully contained early hotspot (CoMoLo study)

Infection fatality rate of SARS-CoV2 in a super-spreading event in Germany

A.C. conceptualized the problem and supervised the overall research. S.T. and C.P.V generated all the figures and wrote the preliminary manuscript. A.C. reviewed and edited the manuscript.

The authors declare no competing interests.