Time-adjusted Analysis Shows Weak Associations Between BCG Vaccination Policy and COVID-19 Disease Progression

Abstract

In this study, we ascertain the associations between BCG vaccination policies and the disease progression in the initial phases of COVID-19 pandemics through analysis of various time-adjusted indicators either directly extracted from the incidence and death reports, or estimated as parameters of disease progression models. We observe weak correlation between BCG vaccination status and indicators related to disease reproduction characteristics. We did not find any associations with case fatality rates (CFR), but the differences in CFR estimates were likely dominated by differences in testing and case reporting between countries. Supplementary material is available through GitHub at https://github.com/fmfi-compbio/ bcg-supplement.

Introduction

The reports on a possible use of the well-established and widely used Bacillus Calmette–Guérin (BCG) vaccine as a protection against COVID-19 (de Vrieze, 2020) raised a lot of interest and media coverage. Several clinical trials have been designed to evaluate the potential of BCG for protection against the SARS-CoV-2 infection in healthcare workers (Bonten, 2020; Khattab, 2020; Curtis, 2020; Cirillo and DiNardo, 2020) . These studies are driven by the so called non-specific effects of BCG vaccine on viral infections, observed in animal models, as well as in humans, although the molecular basis of this phenomenon is not completely understood (Moorlag et al., 2019) .

The associations between BCG vaccination policy and COVID-19 disease progression have also been a subject to contraversy in data analysis, with some studies claiming significant effects on the number of cases and case fatality rates (Miller et al., 2020; Berg et al., 2020) , while others criticizing weaknesses of those studies and claiming no statistically significant differences (Szigeti et al., 2020; Hensel et al., 2020; Fukui et al., 2020; Singh, 2020) .

While correcting for many covariate factors (such as population size, population age distribution, etc.), most of these studies, however, failed to correct for the differences in time progression of the epidemics in each country. COVID-19 epidemic in each country started from relatively few imported cases and in its initial phases spread quickly through exponential growth with high reproduction numbers. At unchecked growth rates, a significant percentage of the country population would be infected before the disease would subside. However, this growth rate only continues until effective measures, such as lockdowns or social distancing policies, are introduced, changing the dynamics of the epidemics substantially, with infection rates rarely reaching a significant percentage of the whole population in the first wave (Flaxman et al., 2020) . In this study, we have estimated a variety of indicators adjusting for time since the beginning of the epidemics in each country, and found that several key indicators show weak, but statistically significant, associations with the BCG vaccination status.

Results

To compare the COVID-19 disease progression between countries with recent universal BCG vaccination policy and those without, several parameters derived from the case and death reports in each country were selected. The parameters reflect early-stage disease spread characteristics (when they are likely not yet affected by social distancing policies), early-stage case fatality rates (before potential effects from overwhelmed health care system), and progression of the disease after the changes characteristic for social distancing policies take effect.

In particular, most of the indicators considered in this study are synchronized based on the reference date, which is the day on which a cumulative number of reported cases surpassed 100. Supplementary table S1 shows the reference date for individual countries and also lists the date when the first nation-wide large scale non-pharmaceutical interventions (NPIs) affecting community spread (e.g. social distancing measures, masking rules, school closures, limitations on large gatherings) were introduced (European Centre for Disease Prevention and Control, 2021; Cheng et al., 2020) (here, we did not consider international travel restrictions, since such restrictions mainly delay the start of the epidemics). In most countries, NPIs were introduced approximately around the reference date; in only a handful of countries the NPIs preceded the reference day by more than a week. Since the effects of NPIs are typically delayed by 10-20 days, and the full effect takes even more time (Nader et al., 2021) , their impact on the indicators used in this study is likely minimal. For the earlystage indicators, we did not normalize for the population size, since the numbers of cases at this stage were very low and the population size was unlikely to pose limitations to the unmitigated disease spread at the time.

Estimates of early stage R are lower in countries with recent BCG vaccination policies. The reproduction number R, the average number of secondary cases of disease caused by a single infected individual, has been estimated using EpiEstim package (Cori et al., 2013) , based on 7-day windows, the first estimate starting on the day when cumulative number of 100 reported cases have been reached (R100), the second estimate starting on 10th day afterwards (R100+10). In many countries, this time period would not reflect the effects of social distancing policies, but would also somewhat avoid the initial period when the case reporting is likely to be unreliable. In both cases, the countries with recent BCG vaccination policies show lower R estimates (Figure 1) and these shifts were statistically significant (Mann Whitney U-test, P = 0:04 for R100 and P = 0:006 for R100+10).

We have also examined the number of days between 10 and 100 reported cases (C10), 100 and 1000 reported cases (C100), 10 and 100 reported deaths (D10), and 100 and 1000 reported deaths (D100)—see details on the indicators in Methods section. These time periods reflect R in various early stages of the epidemic, longer periods meaning slower spread of the disease. Note, that C10 numbers are likely unreliable (due to initial problems in establishing testing and reporting policies in each country), and there are only a few countries that reached 1000 reported deaths before our data set cutoff. Also note that if we assume a constant case fatality rate within a specific time period (typically 6-10 days) and a specific country, and also assume exponential growth in cases within this time period, the numbers D10 and D100 do not actually reflect the death rate, but instead only depend on the underlying value of R. Death reports are likely more accurate than case reports, which are much more affected by testing and reporting policies in each country (Li et al., 2020; Flaxman et al., 2020) . On average, all of these time periods are slightly longer in countries with recent universal BCG policies, with statistically significant results for D10 4.0 3.5 3.0 0 0 1 R 2.5 2.0 1.5 3.5 3.0 2.5 0 1 + 0 0 1 R 2.0 1.5

Canada Quebec●

● Spain

US●● Italy

Germany West● Austria

Sweden● Netherlands●● Switzerland ● Israel Belgium● Denmark● ● Canada Ontario ● Iran Turkey●

● Korea, South

PForrtaungcael●● Norway

● Brazil

Ireland●● Ecuador

Germany East● ● Czechia

Dominican Republic●● Poland

United Kingdom●● Russia

InMdoornoecscioa●● Colombia

● India Mexico● RHoumngaanriay●●●● APelgreuria

● Egypt Philippines●● China Hubei Japan● yes BCG status no

US● Spain●

● Germany West

Switzerland●● Israel

Belgium●● ICtaalnyada Ontario

Austria●● Netherlands Canada Quebec●

Denmark●● Sweden

no ● United Kingdom FRruasnsciae●● Philippines Romania●● Portugal ● India

Turkey●● Peru

Poland●● Brazil Czechia●

MIreelxaincdo●● JGaepramnany East

IndoneIrsaina●●● AClogleormiabia Korea, SEoguytpht●●● EMcouroacdcoor

ChinHauHnguabreyi●●● DNoormwianyican Republic

yes

BCG status

No differences in case fatality rates. We have estimated case fatality rates on days when 100 and 1000 cumulative deaths were first reached in each country (CFR100 and CFR1000 respectively), and also used CMMID methodology (Nishiura et al., 2009; Russel et al., 2020) to correct for estimation of active cases (cCFR100 and cCFR1000) . While some small shifts were observed between countries with and without recent universal BCG vaccination policies, these shifts are not statistically significant. Significant differences in the coefficients of the Vazquez model. One of the difficulties in modelling and predicting the extent of the coronavirus spreading in a population is the divergence of the observed data (the number of confirmed active cases in individual countries) from the trends expected from the traditional SIR type models. Ziff and Ziff (2020) have observed that deaths in China did not follow the typical epidemiological curve and instead of an exponential growth they followed a combined polynomial growth with exponential decay (PGED). Polynomial growth has been also confirmed for multiple other countries (Merrin, 2020) and even though the initial spread in many countries is approximately exponential, it is followed by a steady polynomial growth and in a longer run by an exponential decay (Komarova et al., 2020) .

For a possible explanation of the transition from exponential to polynomial growth, it is natural to look into selfimposed or government-imposed social distancing measures. These measures transform the structure of virus transmitting contact networks in a population, possibly to small-world network structures or even fractal networks.

For a possible explanation of the transition from exponential to polynomial growth, it is natural to look into selfimposed or government-imposed social distancing measures. These measures modify the structure of virus transmitting contact networks in a population primarily by removing long-distance connections within the network. The impact of the growth of the pandemic on the network type was further studied by Medo (2021) on a parametrized network with two limiting cases: the random network with small average distance between the nodes and the regular networks with large average distance between the nodes. The reduction of the long-distance connections seemed to be the main cause of the power-law epidemic growth in the model. The impact of the networkbased contact reduction was also explored by Block et al. (2020) , finding that strategic network-based interactions make the contact reduction more effective.

Moreover, social networks under standard conditions contain a significant fraction of nodes with a high number of connections (that correspond to potential superspreaders). Interestingly, polynomial growth of the number of infections in time in well connected scale-free networks emerges naturally as a consequence of infection initially reaching the highly connected nodes and their neighbors, while their isolation or recovery significantly reduces the interconnectivity of the residual network ( Szabó, 2020 ). Theoretical study of the infection spread in scale-free networks by Vazquez (2006) leads to an explicit formula for the number of infected individuals in time in a form of PGED. The formula contains three key parameters: p - the coefficient of the polynomial growth (not necessary an integer), t - the rate of decay of the exponential tail (1=t is an analogue to the rate of removal of individuals from the infected class to inactive recovered class in the traditional SIR-type models), and A - the constant prefactor (scaling the total population). Based on the value of these parameters, it is straightforward to determine Nmax, the number of infected at the peak of the epidemic, which is independent of the choice of the reference time for the start of the infection. These parameters were obtained by the best fit on the linear scale to the data in each of the considered countries/regions.

Interestingly, we have found that the parameters t and Nmax significantly differ between countries split into two groups—with and without recent universal BCG vaccination policies (Figure 2). The t parameter shifts to the higher values, signifying higher recovery rate in the countries with recent universal BCG vaccination policies (Mann-Whitney U-test P = 0:04). In addition, these countries have generally lower numbers of infected cases at the peak of the epidemic (Nmax) corrected for underreporting (Mann-Whitney U-test P = 0:002).

East and West Germany. The case of Germany is interesting, since the country has been split into East and West Germany in 1949 and reunited in 1990. In East Germany, the policies regarding BCG vaccination followed Eastern Bloc practices, with universal vaccination policy in place between 1951 and 1998. In West Germany, the vaccination has been introduced in 1961, but in 1975 it was discontinued in favor of vaccinating high risk groups only. [The information has been reconstructed from the notes in BCG atlas, however we were not able to confirm this from other sources.] In the present crisis, the whole Germany follows similar practices in case reporting and treatment of the disease. Interestingly, East Germany exhibits much lower estimates of R than West Germany at the corresponding phases of the epidemic (R100 = 2.8, R100+10 = 1.55 in East Germany vs. R100 = 3.14, R100+10 = 2.76 in West Germany; see also Figure 3). Also, the death rate from COVID-19 seems to be significantly lower in East Germany, even when correcting for differences in age distribution (Table 1).

Discussion

While some of the previous studies have observed associations between BCG vaccination policy and spread of COVID-19 (Miller et al., 2020; Berg et al., 2020) , others criticized their work and showed that after corrections for various covariate factors, no statistically significant associations could be found Hensel et al. (2020) ; Fukui et al. (2020) ; Singh (2020) . Most of these studies have used indicators that were quite straightforward, such as the number of reported cases per million inhabitants on a particular date. Here, we have instead chosen a variety of indicators that reflect characteristics of various phases of the epidemics in each country, and moreover, these indicators were implicitly or explicitly adjusted according to the time from the beginning of the epidemic in each country. In fact, we hypothesize that such time adjustment is one of the key factors in such an analysis considering what we know about the spread of COVID-19. Age group A00-A04 A05-A14 A15-A34 A35-A59 A60-A79

A80+

In our data, we have observed several statistically significant associations, and we conclude that there is an association between BCG vaccination policy and spread of COVID-19. However, whether this association is causal or is merely an observed correlation due to some other common factor, is impossible to say. Moreover, most observed shifts in various coefficients are rather small and while the universal BCG vaccination policy may have had a positive impact in some of the countries, the observed impact clearly cannot replace effective policies such as lockdowns and social distancing measures which currently constitute the most effective weapon against the epidemic. At best, the existence of universal BCG vaccination policy may have provided a few days time for governments to effectively institute such policies.

One of the interesting observations is that we did not find any correlation between BCG vaccination policy and CFR. While this may suggest a hypothesis that BCG vaccination may help to limit spread, but may not be effective against difficult progression of the disease in susceptible individuals, we would be careful to draw such conclusions. This is because the estimates of CFR are clearly unreliable at this point of time, with many countries showing CFR estimates well over 10%. Likely, huge differences between countries do not reflect real differences in outcomes of the disease, but rather discrepancies in the amount and effectiveness of testing, with many light or asymptomatic cases remaining undetected. In fact, such a conclusion is partly supported by the evidence from East/West Germany, where we can assume consistent reporting of cases and outcomes, and where differences in CFR seem to be consistent with historical differences in BCG vaccination policies, even after correcting for differences in the age distribution of the population.

Methods

Obtaining case and death reports. The information on reported cases, deaths, and recoveries related to COVID-19 assembled by John Hopkins University Center for System Science and Engineering (Dong et al., 2020) has been downloaded from Humanitarian Data Exchange (Humanitarian Data Exchange, 2020) on April 14, 2020. The data set covers reports from 266 countries from January 22, 2020 until April 13, 2020. For further analysis, only 41 countries with at least 100 reported cumulative deaths have been retained. We also used the data set for Germany maintained by Robert Koch Institute, containing reported cases, deaths, and recoveries split geographically and into age groups; the data set was downloaded through ArcGIS (Robert Koch-Institut and Bundesamt für Kartographie und Geodäsie, 2020) . For our analysis, the data were split geographically into East Germany (Brandeburg, Mecklenburg-Vorpommern, Sachsen, Sachsen-Anhalt, Thüringen, and Berlin) and West Germany (Schleswig-Holstein, Hamburg, Niedersachsen, Bremen, Nordrhein-Westfalen, Hessen, Rheinland-Pfalz, Baden-Württemberg, Bayern, and Saarland).

BCG status of individual countries. For countries included in the study, we have assembled information from the BCG World Atlas (Zwerling et al., 2011) and from the WHO-UNICEF estimates of BCG coverage (World Health Organization, 2020) . Based on this information, the countries were divided into positive BCG status (the countries with current universal BCG vaccination policy and countries with past universal policies discontinued after 1990 or with recent reports of high vaccination coverage from WHO) and negative BCG status (the countries without universal BCG vaccination policy and those that discontinued universal BCG policies and did not satisfy the above conditions); see Supplementary Table S1 for details. Estimation and extraction of indicators. The indicators were extracted from the time series data sets using simple scripts, as outlined in the Results (see Supplementary Material for tables). All of the indicators are computed in time that is relative to a particular milestone, i.e. reaching a particular cumulative number of case reports or death reports. In this way, compared indicators are synchronized at a particular stage of the epidemic. Since the number of cases and deaths is highly dependent on the stage of the epidemic, using such synchronized indicators is a key in our analysis.

Case fatality rate indicators CFR100 and CFR1000 were computed on the days when the cumulative number of reported deaths surpassed 100 and 1000 respectively; lem with CFR indicators is inconsistent reporting on the number of cases in different countries, as this depends highly on testing strategy, reporting methodology, as well as testing capacities of individual countries. Thus, CFR estimates are likely dominated by these factors. We are not aware of any simple method that could overcome this problem at this point of time.

Note that indicators D10 (time from 10 death reports to 100 death reports) and D100 (time from 100 death reports to 1000 death reports), even though based on the numbers of reported deaths, are unlikely to reflect CFR, but instead simply serve as more stable estimates reflecting the underlying reproductive number R. This is because if we assume exponential growth phase and a constant CFR over this period of time, the CFR coefficient will cancel out in the computation of the expected number of days to reach 10-fold increase in the number of deaths.

Indicators R100 and R100+10 were computed using EpiEstim R package (Cori et al., 2013) . This method is based on Bayesian inference, modelling new infections as a Poisson process with rate governed by the instantaneous reproduction number and the number and total infectiousness of infected individuals at the current time interval. The instantaneous reproduction number has a gamma-distributed prior and during the inference is assumed to be constant within each seven-day sliding window to yield an estimate at the end of the window. The infectiousness is approximated by the distribution of the serial interval, which is defined as the time between the onset of symptoms of a case and the onset of symptoms of secondary cases infected by the primary case. Following previous work (Churches, 2020) , we have set the distribution of serial intervals as a discrete gamma distribution with mean of 5 days and standard deviation of 3:4 days. Here, we concentrated on monitoring early stages of the epidemic in each country, when such simple exponential growth model is relatively accurate representation of the spread of the disease. Moreover, the estimated values are used mostly in the non-parametric Mann-Whitney test, which only considers their relative ordering, not exact values.

To avoid initial uncertainty in the reproductive number estimates due to small numbers of case reports, and to adjust for the differences in the start date of epidemics in each country, the seven-day interval for the first estimate (R100) starts on the day when 100 cases have been reported and the second estimate (R100+10) is taken 10 days later. The case incidence numbers have been smoothed over a window of 7 days in order to account for differences in testing procedures on different days of the week (i.e. no or little testing over the weekend in many countries). Such smoothing will not affect the parameters of exponential growth models. It has been verified that confidence intervals at chosen points of time are not unproportionally large.

● Portugal

KoreaC,zSeocuhtiha●● Norway

● Philippines Algeria●

yes no

BCG status the cumulative number of deaths was divided by the cumulative number of reported cases 7 days prior to that date. As alternative indicators for case fatality rates, denoted as cCFR100 and cCFR1000, we have used methodology established by the Centre for the Mathematical Modelling of Infectious Diseases (Nishiura et al., 2009; Russel et al., 2020) , systematically compensating for confirmation-todeath delay using lognormal distribution with mean delay of 13 days and a standard deviation of 12.7 days (Linton et al., 2020) . Regardless of the method, the main probEpidemic curve for Germany East Estimated R for Germany East Epidemic curve for Germany West Estimated R for Germany West

Mar 15

Time Mar 01

Apr 01

Apr 15

Mar 01

Apr 01

Apr 15 7.5 where A, p, and t are parameters and t = 1 (units are days) corresponds to the first day of an infection. In practice, the available data does not report the number of infected individuals in the population due to limited testing availability and potential testing errors. Therefore we use the total number of active cases (confirmed - recovered - deaths) as a proxy for the total number of infected individuals. Bodova and Kollar (2020) (Supplemental information S3) provide a detailed explanation of why the number of active cases is a good approximation of the newly infected individuals in the Vazquez model. In countries with sufficient testing, we assume that the identified active cases represent a constant fraction of the total active cases and the formula for n(t) differs only in the constant factor A.

We used a nonlinear least squares method to infer the parameters in the above relationship from the data. The advantages of fitting the parameters in a linear scale instead of fitting logarithmically transformed data allows us to fit the data globally, as the pre-AGED regime contributes to the errors only by a small amount (Bodova and Kollar, 2020) . However, instead of directly fitting the parameters A, p, and t, we used an equivalent formulation n(t) = Nmax

t Tmax p ep(1 t=Tmax); with parameters Nmax (the maximal number of active cases during the infection), p (the power of the polynomial growth term), and Tmax = p t (the time when the peak is reached). For consistency, we have truncated the data to reduce the impact of testing irregularities during the initial onset of epidemic. Therefore we start the data from the day when a certain number (relative to the population of the country) of active cases Na was reached. The threshold Na was chosen in proportion to the population in the country to reduce effects of randomness in reporting and to account for the spreading potential. Italy served as the reference with a threshold of 200 cases (threshold chosen was always at least 10).

Acknowledgements. This work has been supported by the Slovak Research and Development Agency under the contracts no. APVV-18-0239 (T.V., B.B.), APVV-18-0308 (R.K.), APVV PP-COVID-20-0017 (R.K., K.B.), and by the Scientific Grant Agency of the Slovak Republic under the grants no. 1/0458/18 (T.V.), 1/0463/20 (B.B.), 1/0755/19 (R.K.), and 1/0521/20 (K.B.). This research was also supported by a grant ITMS:313011ATL7 “Pangenomics for personalized clinical management of infected persons based on identified viral genome and human exome” from the Operational Program Integrated Infrastructure (90%) co-financed by the European Regional Development Fund. European Centre for Disease Prevention and Control (2021). Data on country response measures to COVID-19. https://www. ecdc.europa.eu/en/publications-data/ download-data-response-measures-covid-19.

References

Berg , M. K. , Yu , Q. , Salvador , C. E. , Melani , I. , and Kitayama , S. ( 2020 ). Mandated Bacillus CalmetteGuerin (BCG) vaccination predicts flattened curves for the spread of COVID-19 . Science Advances, 6 ( 32 ): eabc1463 .

Block , P. , Hoffman , M. , Raabe , I. J. , Dowd , J. B. , Rahal , C. , Kashyap , R. , and Mills , M. C. ( 2020 ). Social network-based distancing strategies to flatten the COVID-19 curve in a post-lockdown world . Nature Human Behaviour , 4 ( 6 ): 588 - 596 .

Bodova , K. and Kollar , R. ( 2020 ). Emerging algebraic growth trends in SARS-CoV-2 pandemic data . Physical Biology , 17 ( 6 ): 065012 .

Bonten , M. J. M. ( 2020 ). Reducing health care workers absenteeism in Covid-19 pandemic through BCG vaccine (BCG-CORONA) . https://clinicaltrials. gov/ct2/show/NCT04328441.

Cheng , C. , Barcelo , J. , Hartnett , A. S. , Kubinec , R. , and Messerschmidt , L. ( 2020 ). COVID-19 Government Response Event Dataset (CoronaNet v.1.0) . Nature Human Behaviour , 4 ( 7 ): 756 - 768 .

Churches , T. ( 2020 ). COVID-19 epidemiology with R. R Views , An R community blog edited by RStudio . https://rviews.rstudio.com/ 2020 /03/05/ covid-19 - epidemiology-with-r/.

Cirillo , J. D. and DiNardo , A. ( 2020 ). BCG vaccine for health care workers as defense against COVID 19 (BADAS) . https://clinicaltrials.gov/ct2/ show/NCT04348370.

Cori , A. , Ferguson , N. M. , Fraser , C. , and Cauchemez , S. ( 2013 ). A new framework and software to estimate time-varying reproduction numbers during epidemics . American Journal of Epidemiology , 178 ( 9 ): 1505 - 1512 .

Curtis , N. ( 2020 ). BCG vaccination to protect healthcare workers against COVID-19 (BRACE) . https: //clinicaltrials.gov/ct2/show/NCT04327206.

de Vrieze , J. ( 2020 ). Can a century-old TB vaccine steel the immune system against the new coronavirus? Science . doi: 10 .1126/science.abb8297.

Dong , E. , Du , H. , and Gardner , L. ( 2020 ). An interactive web-based dashboard to track COVID-19 in real time . The Lancet Infectious Diseases , 20 ( 5 ): 533 - 534 .

Flaxman , S. , Mishra , S. , Gandy , A. , Unwin , H. J. T. , Mellan , T. A. , Coupland , H. , Whittaker , C. , Zhu , H. , Berah , T. , Eaton , J. W. , Monod , M. , Ghani , A. C. , Donnelly , C. A. , Riley , S. , Vollmer , M. A. C. , Ferguson , N. M. , Okell , L. C. , and Bhatt , S. ( 2020 ). Estimating the effects of non-pharmaceutical interventions on COVID19 in Europe . Nature , 584 ( 7820 ): 257 - 261 .

Fukui , M. , Kawaguchi , K. , and Matsuura , H. ( 2020 ). Does TB vaccination reduce COVID-19 infection?: No evidence from a regression discontinuity analysis . Technical Report doi:10 .1101/ 2020 .04.13.20064287, medRxiv.

Hensel , J. , McAndrews , K. M. , McGrail , D. J. , Dowlatshahi , D. P. , LeBleu , V. S., and Kalluri , R. ( 2020 ). Protection against SARS-CoV-2 by BCG vaccination is not supported by epidemiological analyses . Scientific Reports , 10 ( 1 ): 18377 .

Humanitarian Data Exchange ( 2020 ). Johns Hopkins University Novel Coronavirus (COVID-19) Cases Data . https://data.humdata.org/dataset/ novel-coronavirus-2019 - ncov-cases.

Khattab , A. ( 2020 ). Application of BCG vaccine for immune-prophylaxis among Egyptian healthcare workers during the pandemic of COVID-19 . https:// clinicaltrials.gov/ct2/show/NCT04350931.

Komarova , N. L. , Schang , L. M. , and Wodarz , D. ( 2020 ). Patterns of the COVID-19 pandemic spread around the world: exponential versus power laws . Journal of the Royal Society Interface , 17 ( 170 ): 20200518 .

Li , R. , Pei , S. , Chen , B. , Song , Y. , Zhang, T. , Yang , W. , and Shaman , J. ( 2020 ). Substantial undocumented infection facilitates the rapid dissemination of novel coronavirus (SARS-CoV-2 ). Science, 368 ( 6490 ): 489 - 493 .

Linton , N. M. , Kobayashi , T. , Yang , Y. , Hayashi , K. , Akhmetzhanov , A. R. , Jung , S.-m., Yuan , B. , Kinoshita , R. , and Nishiura , H. ( 2020 ). Incubation period and other epidemiological characteristics of 2019 novel coronavirus infections with right truncation: a statistical analysis of publicly available case data . Journal of Clinical Medicine , 9 ( 2 ): 538 .

Medo , M. ( 2021 ). Contact network models matching the dynamics of the COVID-19 spreading . Journal of Physics A: Mathematical and Theoretical , 54 ( 3 ): 035601 .

Merrin , J. ( 2020 ). Differences in power law growth over time and indicators of COVID-19 pandemic progression worldwide . Physical Biology , 17 ( 6 ): 065005 .

Ziff , A. L. and Ziff , R. M. ( 2020 ). Fractal kinetics of COVID-19 pandemic . International Journal of Educational Excellence , 6 ( 1 ): 43 - 69 .

Miller , A. , Reandelar , M. J. , Fasciglione , K. , Roumenova , V. , Li , Y. , and Otazu , G. H. ( 2020 ). Correlation between universal BCG vaccination policy and reduced morbidity and mortality for COVID19: an epidemiological study . Technical Report doi:10 .1101/ 2020 .03.24.20042937, medRxiv.

Moorlag , S. J. C. F. M. , Arts , R. J. W., van Crevel, R. , and Netea , M. G. ( 2019 ). Non-specific effects of BCG vaccine on viral infections . Clinical Microbiology and Infection , 25 ( 12 ): 1473 - 1478 .

Nader , I. W. , Zeilinger , E. L. , Jomar , D. , and Zauchner , C. ( 2021 ). Onset of effects of non-pharmaceutical interventions on COVID-19 infection rates in 176 countries . BMC Public Health, 21 ( 1 ): 1472 .

Nishiura , H. , Klinkenberg , D. , Roberts , M. , and Heesterbeek , J. A. ( 2009 ). Early epidemiological assessment of the virulence of emerging infectious diseases: a case study of an influenza pandemic . PLoS One , 4 ( 8 ).

Robert Koch-Institut and Bundesamt für Kartographie und Geodäsie ( 2020 ). CSV mit den aktuellen Covid-19 Infektionen pro Tag (Zeitreihe) . https://www.arcgis.com/home/item.html?id= f10774f1c63e40168479a1feb6c7ca74.

Russel , T. , Hellewell , J. , Abbot , S. , et al. ( 2020 ). Using a delay-adjusted case fatality ratio to estimate under-reporting. Available at the Centre for Mathematical Modelling of Infectious Diseases Repository https ://cmmid.github.io/topics/covid19/ global_cfr_estimates.html.

Singh , S. ( 2020 ). BCG vaccines may not reduce COVID-19 mortality rates . Technical Report doi:10 .1101/ 2020 .04.11.20062232, medRxiv.

Szabó , G. M. ( 2020 ). Propagation and mitigation of epidemics in a scale-free network . Technical Report arXiv:2004 . 00067 , arXiv.

Szigeti , R. , Kellermayer , D. , Trakimas , G. , and Kellermayer , R. ( 2020 ). BCG epidemiology supports its protection against COVID-19? A word of caution . PLoS One , 15 ( 10 ): e0240203 .

Vazquez , A. ( 2006 ). Polynomial growth in branching processes with diverging reproductive number . Physical Review Letters , 96 : 038702 .

World Health Organization ( 2020 ). WHO-UNICEF estimates of BCG coverage . https://apps.who. int/immunization_monitoring/globalsummary/ timeseries/tswucoveragebcg.html.

Zwerling , A. , Behr , M. A. , Verma , A. , Brewer , T. F. , Menzies , D. , and Pai , M. ( 2011 ). The BCG World Atlas: a database of global BCG vaccination policies and practices . PLoS Medicine , 8 ( 3 ).