Time-adjusted Analysis Shows Weak Associations Between BCG Vaccination Policy and COVID-19 Disease Progression Katarína Bod’ová1 , Vladimír Boža2 , Broňa Brejová3 , Richard Kollár4 , Katarína Mikušová5 , Tomáš Vinař2 1 Department of Mathematical Analysis and Numerical Mathematics, Faculty of Mathematics, Physics and Informatics, Comenius University, Mlynská dolina, 842 48 Bratislava, Slovakia 2 Department of Applied Informatics, Faculty of Mathematics, Physics and Informatics, Comenius University, Mlynská dolina, 842 48 Bratislava, Slovakia 3 Department of Computer Science, Faculty of Mathematics, Physics and Informatics, Comenius University, Mlynská dolina, 842 48 Bratislava, Slovakia 4 Department of Applied Mathematics and Statistics, Faculty of Mathematics, Physics and Informatics, Comenius University, Mlynská dolina, 842 48 Bratislava, Slovakia 5 Department of Biochemistry, Faculty of Natural Sciences, Comenius University, Ilkovičova 6, 842 48 Bratislava, Slovakia Abstract: In this study, we ascertain the associations be- While correcting for many covariate factors (such as tween BCG vaccination policies and the disease progres- population size, population age distribution, etc.), most sion in the initial phases of COVID-19 pandemics through of these studies, however, failed to correct for the differ- analysis of various time-adjusted indicators either directly ences in time progression of the epidemics in each coun- extracted from the incidence and death reports, or esti- try. COVID-19 epidemic in each country started from rel- mated as parameters of disease progression models. We atively few imported cases and in its initial phases spread observe weak correlation between BCG vaccination status quickly through exponential growth with high reproduc- and indicators related to disease reproduction character- tion numbers. At unchecked growth rates, a significant istics. We did not find any associations with case fatal- percentage of the country population would be infected be- ity rates (CFR), but the differences in CFR estimates were fore the disease would subside. However, this growth rate likely dominated by differences in testing and case report- only continues until effective measures, such as lockdowns ing between countries. or social distancing policies, are introduced, changing the Supplementary material is available through dynamics of the epidemics substantially, with infection GitHub at https://github.com/fmfi-compbio/ rates rarely reaching a significant percentage of the whole bcg-supplement. 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 coun- Introduction try, and found that several key indicators show weak, but statistically significant, associations with the BCG vacci- The reports on a possible use of the well-established and nation status. 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 tri- Results als have been designed to evaluate the potential of BCG for protection against the SARS-CoV-2 infection in health- To compare the COVID-19 disease progression between care workers (Bonten, 2020; Khattab, 2020; Curtis, 2020; countries with recent universal BCG vaccination policy Cirillo and DiNardo, 2020). These studies are driven by and those without, several parameters derived from the the so called non-specific effects of BCG vaccine on viral case and death reports in each country were selected. The infections, observed in animal models, as well as in hu- parameters reflect early-stage disease spread characteris- mans, although the molecular basis of this phenomenon is tics (when they are likely not yet affected by social dis- not completely understood (Moorlag et al., 2019). tancing policies), early-stage case fatality rates (before po- The associations between BCG vaccination policy and tential effects from overwhelmed health care system), and COVID-19 disease progression have also been a subject to progression of the disease after the changes characteristic contraversy in data analysis, with some studies claiming for social distancing policies take effect. significant effects on the number of cases and case fatality In particular, most of the indicators considered in this rates (Miller et al., 2020; Berg et al., 2020), while oth- study are synchronized based on the reference date, which ers criticizing weaknesses of those studies and claiming is the day on which a cumulative number of reported cases no statistically significant differences (Szigeti et al., 2020; surpassed 100. Supplementary table S1 shows the refer- Hensel et al., 2020; Fukui et al., 2020; Singh, 2020). ence date for individual countries and also lists the date ______________ Copyright ©2021 for this paper by its authors. Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0). when the first nation-wide large scale non-pharmaceutical 4.0 ● Iran interventions (NPIs) affecting community spread (e.g. so- Canada Quebec ● Turkey ● cial distancing measures, masking rules, school closures, ● Korea, South limitations on large gatherings) were introduced (Euro- 3.5 ● Spain US ● pean Centre for Disease Prevention and Control, 2021; ● Italy Cheng et al., 2020) (here, we did not consider international Germany West ● Austria Portugal ●●● Norway travel restrictions, since such restrictions mainly delay the France 3.0 ● Brazil start of the epidemics). In most countries, NPIs were in- Ireland ●● Ecuador Sweden ● troduced approximately around the reference date; in only Germany East ● R100 a handful of countries the NPIs preceded the reference day ● Czechia Netherlands ●● Switzerland Dominican Republic ●● Poland by more than a week. Since the effects of NPIs are typi- 2.5 ● Israel Belgium ● United Kingdom ● cally delayed by 10-20 days, and the full effect takes even ● Russia Indonesia ●● Colombia more time (Nader et al., 2021), their impact on the indi- ● Canada Ontario Morocco ● India Mexico ● cators used in this study is likely minimal. For the early- 2.0 Hungary ●● Peru ● Denmark ● Romania ● Algeria stage indicators, we did not normalize for the population ● Egypt size, since the numbers of cases at this stage were very low Philippines ● ● China Hubei and the population size was unlikely to pose limitations to 1.5 Japan ● the unmitigated disease spread at the time. no yes BCG status Estimates of early stage R are lower in countries with re- 3.5 US ● cent BCG vaccination policies. The reproduction num- ber 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 3.0 windows, the first estimate starting on the day when cu- mulative number of 100 reported cases have been reached ● Germany West (R100), the second estimate starting on 10th day after- 2.5 wards (R100+10). In many countries, this time period R100+10 Spain ● ● United Kingdom would not reflect the effects of social distancing policies, France ● Russia ●● Philippines Switzerland Israel ● ● but would also somewhat avoid the initial period when the Belgium ●● Canada Ontario case reporting is likely to be unreliable. In both cases, 2.0 Austria ●● Italy Netherlands Romania ● Portugal ● the countries with recent BCG vaccination policies show Turkey ●● India ● lower R estimates (Figure 1) and these shifts were statisti- Peru Poland ●● Brazil cally significant (Mann Whitney U-test, P = 0.04 for R100 Canada Quebec ● Czechia ● Ireland ●● Germany East ● Mexico ● Japan 1.5 Algeria and P = 0.006 for R100+10). Indonesia Iran ● Colombia ● We have also examined the number of days between 10 Korea, South ●● Morocco Egypt ●● Ecuador Denmark ● Sweden ● China Hubei ●●● Norway Dominican Republic and 100 reported cases (C10), 100 and 1000 reported cases Hungary (C100), 10 and 100 reported deaths (D10), and 100 and no yes BCG status 1000 reported deaths (D100)—see details on the indica- tors in Methods section. These time periods reflect R in various early stages of the epidemic, longer periods mean- Figure 1: Comparison of estimated reproduction numbers ing slower spread of the disease. Note, that C10 num- R100 (top) and R100+10 (bottom) between countries with bers are likely unreliable (due to initial problems in es- and without the universal BCG vaccination policy. tablishing testing and reporting policies in each country), and there are only a few countries that reached 1000 re- ported deaths before our data set cutoff. Also note that if (Mann-Whitney U-test, P = 0.02). 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 No differences in case fatality rates. We have estimated time period, the numbers D10 and D100 do not actually case fatality rates on days when 100 and 1000 cumulative reflect the death rate, but instead only depend on the un- deaths were first reached in each country (CFR100 and derlying value of R. Death reports are likely more accurate CFR1000 respectively), and also used CMMID methodol- than case reports, which are much more affected by test- ogy (Nishiura et al., 2009; Russel et al., 2020) to correct ing and reporting policies in each country (Li et al., 2020; for estimation of active cases (cCFR100 and cCFR1000) Flaxman et al., 2020). On average, all of these time pe- . While some small shifts were observed between coun- riods are slightly longer in countries with recent universal tries with and without recent universal BCG vaccination BCG policies, with statistically significant results for D10 policies, these shifts are not statistically significant. Significant differences in the coefficients of the Vazquez eters, it is straightforward to determine Nmax, the number model. One of the difficulties in modelling and predict- of infected at the peak of the epidemic, which is indepen- ing the extent of the coronavirus spreading in a popula- dent of the choice of the reference time for the start of the tion is the divergence of the observed data (the number of infection. These parameters were obtained by the best fit confirmed active cases in individual countries) from the on the linear scale to the data in each of the considered trends expected from the traditional SIR type models. Ziff countries/regions. and Ziff (2020) have observed that deaths in China did Interestingly, we have found that the parameters τ and not follow the typical epidemiological curve and instead Nmax significantly differ between countries split into two of an exponential growth they followed a combined poly- groups—with and without recent universal BCG vacci- nomial growth with exponential decay (PGED). Polyno- nation policies (Figure 2). The τ parameter shifts to mial growth has been also confirmed for multiple other the higher values, signifying higher recovery rate in the countries (Merrin, 2020) and even though the initial spread countries with recent universal BCG vaccination policies in many countries is approximately exponential, it is fol- (Mann-Whitney U-test P = 0.04). In addition, these coun- lowed by a steady polynomial growth and in a longer run tries have generally lower numbers of infected cases at the by an exponential decay (Komarova et al., 2020). peak of the epidemic (Nmax) corrected for underreporting For a possible explanation of the transition from expo- (Mann-Whitney U-test P = 0.002). nential to polynomial growth, it is natural to look into self- imposed or government-imposed social distancing mea- East and West Germany. The case of Germany is inter- sures. These measures transform the structure of virus esting, since the country has been split into East and West transmitting contact networks in a population, possibly to Germany in 1949 and reunited in 1990. In East Germany, small-world network structures or even fractal networks. the policies regarding BCG vaccination followed Eastern For a possible explanation of the transition from expo- Bloc practices, with universal vaccination policy in place nential to polynomial growth, it is natural to look into self- between 1951 and 1998. In West Germany, the vaccina- imposed or government-imposed social distancing mea- tion has been introduced in 1961, but in 1975 it was dis- sures. These measures modify the structure of virus continued in favor of vaccinating high risk groups only. transmitting contact networks in a population primarily [The information has been reconstructed from the notes by removing long-distance connections within the net- in BCG atlas, however we were not able to confirm this work. The impact of the growth of the pandemic on the from other sources.] In the present crisis, the whole Ger- network type was further studied by Medo (2021) on a many follows similar practices in case reporting and treat- parametrized network with two limiting cases: the random ment of the disease. Interestingly, East Germany exhibits network with small average distance between the nodes much lower estimates of R than West Germany at the cor- and the regular networks with large average distance be- responding phases of the epidemic (R100 = 2.8, R100+10 tween the nodes. The reduction of the long-distance con- = 1.55 in East Germany vs. R100 = 3.14, R100+10 = 2.76 nections seemed to be the main cause of the power-law in West Germany; see also Figure 3). Also, the death rate epidemic growth in the model. The impact of the network- from COVID-19 seems to be significantly lower in East based contact reduction was also explored by Block et al. Germany, even when correcting for differences in age dis- (2020), finding that strategic network-based interactions tribution (Table 1). make the contact reduction more effective. Moreover, social networks under standard conditions Discussion contain a significant fraction of nodes with a high number of connections (that correspond to potential superspread- While some of the previous studies have observed asso- ers). Interestingly, polynomial growth of the number of ciations between BCG vaccination policy and spread of infections in time in well connected scale-free networks COVID-19 (Miller et al., 2020; Berg et al., 2020), others emerges naturally as a consequence of infection initially criticized their work and showed that after corrections for reaching the highly connected nodes and their neighbors, various covariate factors, no statistically significant asso- while their isolation or recovery significantly reduces the ciations could be found Hensel et al. (2020); Fukui et al. interconnectivity of the residual network (Szabó, 2020). (2020); Singh (2020). Most of these studies have used in- Theoretical study of the infection spread in scale-free net- dicators that were quite straightforward, such as the num- works by Vazquez (2006) leads to an explicit formula for ber of reported cases per million inhabitants on a particu- the number of infected individuals in time in a form of lar date. Here, we have instead chosen a variety of indi- PGED. The formula contains three key parameters: p - the cators that reflect characteristics of various phases of the coefficient of the polynomial growth (not necessary an in- epidemics in each country, and moreover, these indicators teger), τ - the rate of decay of the exponential tail (1/τ is were implicitly or explicitly adjusted according to the time an analogue to the rate of removal of individuals from the from the beginning of the epidemic in each country. In infected class to inactive recovered class in the traditional fact, we hypothesize that such time adjustment is one of SIR-type models), and A - the constant prefactor (scaling the key factors in such an analysis considering what we the total population). Based on the value of these param- know about the spread of COVID-19. East Germany West Germany Age CMMDI adj. CMMDI adj. group Deaths Cases CFR Cases cCFR Deaths Cases CFR Cases cCFR A00-A04 0 126 0 79 0 1 863 0.0012 535 0.0019 A05-A14 0 303 0 196 0 0 2152 0 1396 0 A15-A34 0 3589 0 2416 0 6 26771 0.0002 17632 0.0003 A35-A59 13 6022 0.0022 4036 0.0032 123 48761 0.0025 32777 0.0038 A60-A79 77 2459 0.0313 1565 0.0492 893 21700 0.0412 13795 0.0647 A80+ 143 1148 0.1246 601 0.2378 1709 10951 0.1561 5712 0.2992 Table 1: Differences in CFR in different age groups between East Germany and West Germany. Both raw CFR values and values corrected by CMMDI methodology are presented. In our data, we have observed several statistically sig- data set covers reports from 266 countries from Jan- nificant associations, and we conclude that there is an as- uary 22, 2020 until April 13, 2020. For further anal- sociation between BCG vaccination policy and spread of ysis, only 41 countries with at least 100 reported cu- COVID-19. However, whether this association is causal or mulative deaths have been retained. We also used the is merely an observed correlation due to some other com- data set for Germany maintained by Robert Koch Insti- mon factor, is impossible to say. Moreover, most observed tute, containing reported cases, deaths, and recoveries shifts in various coefficients are rather small and while the split geographically and into age groups; the data set was universal BCG vaccination policy may have had a posi- downloaded through ArcGIS (Robert Koch-Institut and tive impact in some of the countries, the observed impact Bundesamt für Kartographie und Geodäsie, 2020). For clearly cannot replace effective policies such as lockdowns our analysis, the data were split geographically into East and social distancing measures which currently constitute Germany (Brandeburg, Mecklenburg-Vorpommern, Sach- the most effective weapon against the epidemic. At best, sen, Sachsen-Anhalt, Thüringen, and Berlin) and West the existence of universal BCG vaccination policy may Germany (Schleswig-Holstein, Hamburg, Niedersachsen, have provided a few days time for governments to effec- Bremen, Nordrhein-Westfalen, Hessen, Rheinland-Pfalz, tively institute such policies. Baden-Württemberg, Bayern, and Saarland). One of the interesting observations is that we did not find any correlation between BCG vaccination policy and BCG status of individual countries. For countries included CFR. While this may suggest a hypothesis that BCG vac- in the study, we have assembled information from the cination may help to limit spread, but may not be effec- BCG World Atlas (Zwerling et al., 2011) and from the tive against difficult progression of the disease in suscep- WHO-UNICEF estimates of BCG coverage (World Health tible individuals, we would be careful to draw such con- Organization, 2020). Based on this information, the coun- clusions. This is because the estimates of CFR are clearly tries were divided into positive BCG status (the countries unreliable at this point of time, with many countries show- with current universal BCG vaccination policy and coun- ing CFR estimates well over 10%. Likely, huge differ- tries with past universal policies discontinued after 1990 ences between countries do not reflect real differences in or with recent reports of high vaccination coverage from outcomes of the disease, but rather discrepancies in the WHO) and negative BCG status (the countries without amount and effectiveness of testing, with many light or universal BCG vaccination policy and those that discontin- asymptomatic cases remaining undetected. In fact, such a ued universal BCG policies and did not satisfy the above conclusion is partly supported by the evidence from East- conditions); see Supplementary Table S1 for details. /West Germany, where we can assume consistent report- ing of cases and outcomes, and where differences in CFR Estimation and extraction of indicators. The indicators seem to be consistent with historical differences in BCG were extracted from the time series data sets using sim- vaccination policies, even after correcting for differences ple scripts, as outlined in the Results (see Supplementary in the age distribution of the population. 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 re- Methods ports. In this way, compared indicators are synchronized at a particular stage of the epidemic. Since the number of Obtaining case and death reports. The information on re- cases and deaths is highly dependent on the stage of the ported cases, deaths, and recoveries related to COVID-19 epidemic, using such synchronized indicators is a key in assembled by John Hopkins University Center for Sys- our analysis. tem Science and Engineering (Dong et al., 2020) has been Case fatality rate indicators CFR100 and CFR1000 downloaded from Humanitarian Data Exchange (Human- were computed on the days when the cumulative number itarian Data Exchange, 2020) on April 14, 2020. The of reported deaths surpassed 100 and 1000 respectively; ● Netherlands lem with CFR indicators is inconsistent reporting on the ● Philippines number of cases in different countries, as this depends Norway ● highly on testing strategy, reporting methodology, as well as testing capacities of individual countries. Thus, CFR 14 estimates are likely dominated by these factors. We are US ● ● Czechia not aware of any simple method that could overcome this problem at this point of time. Portugal ● ● Algeria Note that indicators D10 (time from 10 death reports to tau 10 ● Italy 100 death reports) and D100 (time from 100 death reports Korea, South ● to 1000 death reports), even though based on the numbers Belgium ● of reported deaths, are unlikely to reflect CFR, but instead simply serve as more stable estimates reflecting the under- lying reproductive number R. This is because if we as- 6 ● Spain sume exponential growth phase and a constant CFR over this period of time, the CFR coefficient will cancel out in Switzerland ●● Israel the computation of the expected number of days to reach Austria ● 10-fold increase in the number of deaths. no yes BCG status 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 Italy ● as a Poisson process with rate governed by the instan- 90000 taneous reproduction number and the number and total ● Spain infectiousness of infected individuals at the current time interval. The instantaneous reproduction number has a gamma-distributed prior and during the inference is as- sumed to be constant within each seven-day sliding win- 60000 dow to yield an estimate at the end of the window. The Nmax 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 30000 Netherlands ● of secondary cases infected by the primary case. Follow- ing previous work (Churches, 2020), we have set the dis- ● Belgium ● Portugal tribution of serial intervals as a discrete gamma distribu- Switzerland ● Israel ● Austria tion with mean of 5 days and standard deviation of 3.4 Korea, South ●●● Norway Czechia ● Philippines days. Here, we concentrated on monitoring early stages 0 Algeria ● of the epidemic in each country, when such simple expo- no yes BCG status nential growth model is relatively accurate representation of the spread of the disease. Moreover, the estimated val- Figure 2: Comparison of Vazquez model parameters esti- ues are used mostly in the non-parametric Mann-Whitney mated for countries with and without universal BCG vac- test, which only considers their relative ordering, not exact cination policies. Top: Rate of decay of the exponential values. tail (τ). Bottom: Number of infected cases at the peak of To avoid initial uncertainty in the reproductive number the epidemic corrected for underreporting (Nmax). estimates due to small numbers of case reports, and to ad- just for the differences in the start date of epidemics in each country, the seven-day interval for the first estimate the cumulative number of deaths was divided by the cumu- (R100) starts on the day when 100 cases have been re- lative number of reported cases 7 days prior to that date. ported and the second estimate (R100+10) is taken 10 days As alternative indicators for case fatality rates, denoted as later. The case incidence numbers have been smoothed cCFR100 and cCFR1000, we have used methodology es- over a window of 7 days in order to account for differ- tablished by the Centre for the Mathematical Modelling ences in testing procedures on different days of the week of Infectious Diseases (Nishiura et al., 2009; Russel et al., (i.e. no or little testing over the weekend in many coun- 2020), systematically compensating for confirmation-to- tries). Such smoothing will not affect the parameters of death delay using lognormal distribution with mean delay exponential growth models. It has been verified that confi- of 13 days and a standard deviation of 12.7 days (Linton dence intervals at chosen points of time are not unpropor- et al., 2020). Regardless of the method, the main prob- tionally large. Epidemic curve for Germany East Estimated R for Germany East 5000 7.5 4000 Incidence 3000 5.0 R 2000 2.5 1000 0 0.0 Mar 01 Mar 15 Apr 01 Apr 15 Mar 01 Mar 15 Apr 01 Apr 15 Time Time Epidemic curve for Germany West Estimated R for Germany West 5000 7.5 4000 Incidence 3000 5.0 R 2000 2.5 1000 0 0.0 Mar 01 Mar 15 Apr 01 Apr 15 Mar 01 Mar 15 Apr 01 Apr 15 Time Time Figure 3: Comparison of COVID-19 epidemic progression between East Germany and West Germany. Reproduction numbers R were estimated using seven day windows using smoothed incidence numbers. Application of Vazquez model. The number of new rameters A, p, and τ, we used an equivalent formulation infected individuals at time t in the Vazquez model  t p (Vazquez, 2006; Ziff and Ziff, 2020) has the form n(t) = Nmax · · e p(1−t/Tmax) , Tmax A  t p − t with parameters Nmax (the maximal number of active n(t) = · ·e τ , cases during the infection), p (the power of the polyno- τ τ mial growth term), and Tmax = p ∗ τ (the time when the peak is reached). For consistency, we have truncated the where A, p, and τ are parameters and t = 1 (units are days) data to reduce the impact of testing irregularities during corresponds to the first day of an infection. In practice, the the initial onset of epidemic. Therefore we start the data available data does not report the number of infected indi- from the day when a certain number (relative to the popu- viduals in the population due to limited testing availability lation of the country) of active cases Na was reached. The and potential testing errors. Therefore we use the total threshold Na was chosen in proportion to the population number of active cases (confirmed - recovered - deaths) in the country to reduce effects of randomness in report- as a proxy for the total number of infected individuals. ing and to account for the spreading potential. Italy served Bodova and Kollar (2020) (Supplemental information S3) as the reference with a threshold of 200 cases (threshold provide a detailed explanation of why the number of ac- chosen was always at least 10). tive cases is a good approximation of the newly infected individuals in the Vazquez model. In countries with suf- Acknowledgements. This work has been supported by the ficient testing, we assume that the identified active cases Slovak Research and Development Agency under the con- represent a constant fraction of the total active cases and tracts no. APVV-18-0239 (T.V., B.B.), APVV-18-0308 the formula for n(t) differs only in the constant factor A. (R.K.), APVV PP-COVID-20-0017 (R.K., K.B.), and by We used a nonlinear least squares method to infer the the Scientific Grant Agency of the Slovak Republic un- parameters in the above relationship from the data. The der the grants no. 1/0458/18 (T.V.), 1/0463/20 (B.B.), advantages of fitting the parameters in a linear scale in- 1/0755/19 (R.K.), and 1/0521/20 (K.B.). This research stead of fitting logarithmically transformed data allows was also supported by a grant ITMS:313011ATL7 “Pange- us to fit the data globally, as the pre-AGED regime con- nomics for personalized clinical management of infected tributes to the errors only by a small amount (Bodova and persons based on identified viral genome and human ex- Kollar, 2020). However, instead of directly fitting the pa- ome” from the Operational Program Integrated Infrastruc- ture (90%) co-financed by the European Regional Devel- European Centre for Disease Prevention and opment Fund. Control (2021). Data on country response measures to COVID-19. https://www. ecdc.europa.eu/en/publications-data/ References download-data-response-measures-covid-19. Berg, M. K., Yu, Q., Salvador, C. E., Melani, I., and Flaxman, S., Mishra, S., Gandy, A., Unwin, H. J. T., Mel- Kitayama, S. (2020). Mandated Bacillus Calmette- lan, T. A., Coupland, H., Whittaker, C., Zhu, H., Berah, Guerin (BCG) vaccination predicts flattened curves T., Eaton, J. W., Monod, M., Ghani, A. C., Donnelly, for the spread of COVID-19. Science Advances, C. A., Riley, S., Vollmer, M. A. C., Ferguson, N. M., 6(32):eabc1463. Okell, L. C., and Bhatt, S. (2020). Estimating the ef- fects of non-pharmaceutical interventions on COVID- Block, P., Hoffman, M., Raabe, I. J., Dowd, J. B., Ra- 19 in Europe. Nature, 584(7820):257–261. hal, C., Kashyap, R., and Mills, M. C. (2020). So- cial network-based distancing strategies to flatten the Fukui, M., Kawaguchi, K., and Matsuura, H. (2020). COVID-19 curve in a post-lockdown world. Nature Hu- Does TB vaccination reduce COVID-19 infection?: man Behaviour, 4(6):588–596. No evidence from a regression discontinuity analysis. Technical Report doi:10.1101/2020.04.13.20064287, Bodova, K. and Kollar, R. (2020). Emerging algebraic medRxiv. growth trends in SARS-CoV-2 pandemic data. Physical Biology, 17(6):065012. Hensel, J., McAndrews, K. M., McGrail, D. J., Dowlat- shahi, D. P., LeBleu, V. S., and Kalluri, R. (2020). Pro- Bonten, M. J. M. (2020). Reducing health care workers tection against SARS-CoV-2 by BCG vaccination is not absenteeism in Covid-19 pandemic through BCG vac- supported by epidemiological analyses. Scientific Re- cine (BCG-CORONA). https://clinicaltrials. ports, 10(1):18377. gov/ct2/show/NCT04328441. Humanitarian Data Exchange (2020). Johns Hopkins Cheng, C., Barcelo, J., Hartnett, A. S., Kubinec, R., and University Novel Coronavirus (COVID-19) Cases Messerschmidt, L. (2020). COVID-19 Government Re- Data. https://data.humdata.org/dataset/ sponse Event Dataset (CoronaNet v.1.0). Nature Human novel-coronavirus-2019-ncov-cases. Behaviour, 4(7):756–768. Khattab, A. (2020). Application of BCG vaccine for Churches, T. (2020). COVID-19 epidemiology with immune-prophylaxis among Egyptian healthcare work- R. R Views, An R community blog edited by RStu- ers during the pandemic of COVID-19. https:// dio. https://rviews.rstudio.com/2020/03/05/ clinicaltrials.gov/ct2/show/NCT04350931. covid-19-epidemiology-with-r/. Komarova, N. L., Schang, L. M., and Wodarz, D. (2020). Cirillo, J. D. and DiNardo, A. (2020). BCG vaccine Patterns of the COVID-19 pandemic spread around the for health care workers as defense against COVID world: exponential versus power laws. Journal of the 19 (BADAS). https://clinicaltrials.gov/ct2/ Royal Society Interface, 17(170):20200518. show/NCT04348370. Li, R., Pei, S., Chen, B., Song, Y., Zhang, T., Yang, W., Cori, A., Ferguson, N. M., Fraser, C., and Cauchemez, and Shaman, J. (2020). Substantial undocumented in- S. (2013). A new framework and software to estimate fection facilitates the rapid dissemination of novel coro- time-varying reproduction numbers during epidemics. navirus (SARS-CoV-2). Science, 368(6490):489–493. American Journal of Epidemiology, 178(9):1505–1512. Linton, N. M., Kobayashi, T., Yang, Y., Hayashi, K., Curtis, N. (2020). BCG vaccination to protect health- Akhmetzhanov, A. R., Jung, S.-m., Yuan, B., Kinoshita, care workers against COVID-19 (BRACE). https: R., and Nishiura, H. (2020). Incubation period and other //clinicaltrials.gov/ct2/show/NCT04327206. epidemiological characteristics of 2019 novel coron- avirus infections with right truncation: a statistical anal- de Vrieze, J. (2020). Can a century-old TB vaccine steel ysis of publicly available case data. Journal of Clinical the immune system against the new coronavirus? Sci- Medicine, 9(2):538. ence. doi:10.1126/science.abb8297. Medo, M. (2021). Contact network models matching Dong, E., Du, H., and Gardner, L. (2020). An interactive the dynamics of the COVID-19 spreading. Jour- web-based dashboard to track COVID-19 in real time. nal of Physics A: Mathematical and Theoretical, The Lancet Infectious Diseases, 20(5):533–534. 54(3):035601. Merrin, J. (2020). Differences in power law growth over Ziff, A. L. and Ziff, R. M. (2020). Fractal kinetics of time and indicators of COVID-19 pandemic progression COVID-19 pandemic. International Journal of Educa- worldwide. Physical Biology, 17(6):065005. tional Excellence, 6(1):43–69. Miller, A., Reandelar, M. J., Fasciglione, K., Roumen- Zwerling, A., Behr, M. A., Verma, A., Brewer, T. F., Men- ova, V., Li, Y., and Otazu, G. H. (2020). Cor- zies, D., and Pai, M. (2011). The BCG World Atlas: a relation between universal BCG vaccination policy database of global BCG vaccination policies and prac- and reduced morbidity and mortality for COVID- tices. PLoS Medicine, 8(3). 19: 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 inter- ventions on COVID-19 infection rates in 176 countries. BMC Public Health, 21(1):1472. Nishiura, H., Klinkenberg, D., Roberts, M., and Heester- beek, 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 Kartogra- phie und Geodäsie (2020). CSV mit den ak- tuellen Covid-19 Infektionen pro Tag (Zeitreihe). https://www.arcgis.com/home/item.html?id= f10774f1c63e40168479a1feb6c7ca74. Russel, T., Hellewell, J., Abbot, S., et al. (2020). Us- ing a delay-adjusted case fatality ratio to estimate under-reporting. Available at the Centre for Mathe- matical Modelling of Infectious Diseases Repository https://cmmid.github.io/topics/covid19/ global_cfr_estimates.html. Singh, S. (2020). BCG vaccines may not re- duce 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 Keller- mayer, R. (2020). BCG epidemiology supports its pro- tection against COVID-19? A word of caution. PLoS One, 15(10):e0240203. Vazquez, A. (2006). Polynomial growth in branching pro- cesses 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.