=Paper=
{{Paper
|id=Vol-2962/paper17
|storemode=property
|title=Time-adjusted Analysis Shows Weak Associations Between BCG Vaccination Policy and COVID-19 Disease Progression
|pdfUrl=https://ceur-ws.org/Vol-2962/paper17.pdf
|volume=Vol-2962
|authors=Katarína Boďová,Vladimír Boža,Broňa Brejová,Richard Kollár,Katarína Mikušová,Tomáš Vinař
|dblpUrl=https://dblp.org/rec/conf/itat/BodovaBBKMV21
}}
==Time-adjusted Analysis Shows Weak Associations Between BCG Vaccination Policy and COVID-19 Disease Progression==
https://ceur-ws.org/Vol-2962/paper17.pdf
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.