\documentclass[10pt]{article}
\usepackage[hmargin=1.5cm,top=2cm,bottom=2cm]{geometry}
\usepackage{mathtext}
\usepackage{multicol}
\setlength\columnsep{15pt}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{array}
\usepackage{booktabs}
\usepackage{tabularx}
\usepackage[auth-sc]{authblk}
\usepackage{longtable}
\usepackage{multirow}
\usepackage{hyperref}
\usepackage{enumerate}
\usepackage[labelfont=bf]{caption}
\usepackage[usenames,dvipsnames]{xcolor}
\usepackage{mdframed}
\usepackage{graphics}
\usepackage{multirow}
\usepackage{rotating}
\usepackage{array}
\usepackage{capt-of}
\usepackage{lscape}
\usepackage{caption}
\usepackage{breakurl}
\usepackage{todonotes}
\usepackage{hanging}
\usepackage{pagecolor}
\usepackage[final]{pdfpages}
\usepackage[leftFloats,CaptionAfterwards]{fltpage}
\usepackage[numbers,super,sort&compress]{natbib}
\setlength{\bibsep}{0pt plus 0.3ex}
\usepackage{abstract}
\usepackage{enumitem}
\usepackage{soul}
\usepackage{titlesec}
\titleformat{\section}[block]{\large\bfseries\filcenter}{\thesection.}{0.4em}{}
\titleformat{\subsection}[block]{\normalsize\sc\bfseries\filcenter}{\thesubsection.}{0.4em}{}
\titleformat{\subsubsection}[block]{\normalsize\sc\itshape\filright}{\thesubsection.}{0.4em}{}
\setcounter{secnumdepth}{5}

\usepackage[hang]{footmisc}
\setlength\footnotemargin{0em}

\setlength{\skip\footins}{0.75cm}

\makeatletter
\renewcommand\footnoterule{%
  \kern-3\p@
  \hrule\@width \textwidth height 1.5pt
  \kern2.6\p@}
\makeatother

\makeatletter
\def\@biblabel#1{\@ifnotempty{#1}{#1.}}
\makeatother

\newcommand{\filllastline}[1]{
\setlength\leftskip{0pt}
\setlength\rightskip{0pt}
\setlength\parfillskip{0pt}
#1}

\newenvironment{Figure}
{\par\medskip\noindent\minipage{\linewidth}}
{\endminipage\par\medskip}


\title{\bf The impact and progression of the COVID-19 pandemic in Bulgaria in its first two years}
\renewcommand\Authfont{\scshape\normalsize}
% \author[5]{LIST ADDITIONAL AUTHORS}
\author[1,2,$\#$]{Antoni Rangachev}
\author[3,$\#$]{Georgi K. Marinov}
\author[4]{Mladen Mladenov}
\renewcommand\Affilfont{\itshape\normalsize}
\affil[1]{Institute of Mathematics and Informatics, Bulgarian Academy of Sciences, Sofia 1113, Bulgaria}
\affil[2]{International Center for Mathematical Sciences-Sofia}
\affil[3]{Department of Genetics, Stanford University, Stanford, CA 94305, USA}
\affil[4]{Tilburg University, Warandelaan 2, 5037 AB Tilburg, Netherlands}
% \affil[5]{National Reference Laboratory of HIV, National Center of Infectious and Parasitic Diseases, 1504 Sofia, Bulgaria}
\affil[$\#$]{Corresponding author}
\affil[*]{These authors contributed equally}
\date{}

\begin{document}
\maketitle

% \centerline{}
% \centerline{}

\begin{abstract}
\noindent {\normalsize \textbf{After initially having low levels of SARS-CoV-2 infections for much of the year, at the end of 2020 Bulgaria experienced a major epidemic surge, which caused the highest recorded excess mortality in Europe and among the highest in the word (Excess Mortality Rate, or EMR $\sim$0.25\%). Two more major waves followed in 2021, followed by another one in early 2022. In this study we analyze the temporal and spatial patterns of excess mortality at the national and local levels and across different demographic groups in Bulgaria, and compare those at the European level. The country has continued to exhibit the previous pattern of extremely high excess mortality as measured both by crude mortality metrics (EMR $\sim$1.05\% up to the end of March 2022) and by standardized ones -- Potential Years of Life Lost (PYLL) and Aged-Standardized Years of life lost Rate (ASYR). Unlike Western Europe, the bulk of excess mortality in Bulgaria, as well as in several other countries in Eastern Europe, occurred in the second year of the pandemic, likely related to the differences in the levels of vaccination coverage between these regions. We also observe even more extreme levels of excess mortality at the regional level and in some subpopulations (e.g. total EMR values for males $\geq$2\% and EMR values for males aged 40-64 $\geq$1\% in certain areas). We discuss these observations in light of the estimates of infection fatality rate (IFR) and eventual population fatality rate (PFR) made early in the course of the pandemic.} }
\centerline{}
\centerline{}
\end{abstract}

\section*{Introduction}

The SARS-CoV-2 virus and the COVID-19 disease\cite{Wang2020,Zhou2020,Huang2020} that it causes have triggered the most significant acute public health crisis in more than a century. SARS-CoV-2 has spread widely in most countries around the world, and has been the driver of substantial excess mortality in many of them\cite{Karlinsky2021,Pifarre2021}. 

The pandemic took divergent trajectories in different regions of the world, initially depending on the timing of the imposition of containment measures relative to the undetected early cryptic spread of the virus, and later based on some combination of the relaxation of these measures, seasonal effects, the build up/waning of population immunity, the appearance of new variants of SARS-CoV-2 that are more contagious and/or antigenically divergent, and other factors. Some countries were heavily affected early on and then experienced further major epidemic waves, others were only hard hit at later stages of the pandemic.

By the end of 2020, Bulgaria emerged as one of the countries experiencing among the highest pandemic-related excess mortality in the world, even though it was one of the early containment success stories in the course of the pandemic, largely escaping the first major wave that affected greatly many areas in Western Europe, and the Americas. As a previous analysis of ours has shown\cite{Rangachev2021}, the EMR value for the country by January 1st 2021 stood at $\sim$0.25\% (more than twice the official death count, due to some combination of insufficient testing, registration of COVID deaths as having occurred due to other reasons, and elevated mortality from otherwise treatable other conditions due to hospital capacity being exceeded). 

Subsequently, the country experienced three more major waves, in March-April 2021, in the last few months of 2021, and early in 2022. In this study we track the development and assess the impact of the pandemic on different demographic groups and regions in Bulgaria up to the end of March 2022, using a combination of excess mortality analyses and SARS-CoV-2 genome sequencing surveillance. 

 \begin{figure*}[!ht]
 \begin{center}
 \includegraphics[width=18.5cm]{Fig1-ASYR-Europe.png}
 \captionsetup{singlelinecheck=off,justification=justified}
 \caption{
 {\bf Excess mortality in Europe and Bulgaria during the COVID-19 pandemic (up to the end of March 2022)}. 
 (A) Standardized ASYR values, total; 
 (B) Standardized ASYR values, females; 
 (C) Standardized ASYR values, males;
}
\label{Fig1}
\end{center}
 \end{figure*}

\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=18.5cm]{Fig3-sequencing-V2.png}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Evolution of the SARS-CoV-2 variant composition in the course of the COVID-19 epidemic in Bulgaria up to March 2022}. Shown is the fraction of sequenced genomes belonging to each of the listed variants for each month since the beginning of December 2020. The total number of sequenced SARS-CoV-2 genomes is shown on top.}
\label{Fig3}
\end{center}
\end{figure*}


\begin{figure*}
\begin{center}
\begin{minipage}[c]{0.70\linewidth}
\includegraphics[width=12.4cm]{Fig5-trajectory.png}
\end{minipage}\hfill
\begin{minipage}[c]{0.30\linewidth}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Temporal trajectory of the COVID-19 pandemic in Bulgaria}. 
(A) Confirmed cases over time (weekly);
(B) Officially reported and excess deaths over time (weekly);
(C) CFR over time. Data on reported cases and deaths was obtained from the Our World In Data website\cite{OWID};
(D) Evolution of the undercount ratio (excess mortality divided by official COVID deaths) over time (note that the periods between waves, for which estimates of excess mortality are uncertain, are omitted from the graph).
}
\label{Fig5}
\end{minipage}
\end{center}
\end{figure*}

These subsequent waves have dramatically increased the excess mortality burden in the country, and as a result it has become the first one (among those for which overall mortality data is available) where COVID-related excess deaths have exceeded 1\% of the total population. Furthermore, we, continuing the trend established previously\cite{Rangachev2021}, observe major discrepancies between the outcomes within the country. EMR values in some regions are now approaching 2\%, and they have exceeded that value for males in certain areas. In addition, mortality in the working age 40-64 group is approaching or has even exceeded 1\%, a surprising result considering the commonly assumed dramatic age skew of COVID-related mortality. Despite the reduced Case Fatality Ratio (CFR) associated with the newly emerged at the end of 2021 Omicron variant, considerable excess mortality, not captured by official COVID death statistics, persisted in the first months of 2022. These patterns are in stark contrast to those observed in countries in Western Europe, where excess mortality was concentrated in 2020 and decreased in 2021. They are, however, shared with most other countries in Eastern Europe, although Bulgaria still exhibits the most extreme excess mortality figures. The likely explanation for this pattern is the lower vaccination rates in Eastern Europe and particularly in Bulgaria. Finally, we discuss these findings in the context of the commonly cited figures for the infection fatality rate (IFR) of COVID-19.

\section*{Results}

\subsection*{Loss of life as a result of the COVID-19 pandemic}

In order to evaluate the impact of the COVID-19 pandemic on different countries in Europe we applied excess mortality analysis for the period from the start of the pandemic until the end of March 2022, following previously established methods\cite{Rangachev2021,Karlinsky2021} (see the Methods section for details). Excess mortality measures are more objective measures of pandemic impact as officially recorded COVID mortality is often not an accurate representation of reality, due to insufficient availability of testing, inaccurate reporting, and other factors, such as second-order impacts of COVID infections (i.e. overwhelmed healthcare systems not being able to provide adequate treatment) leading to fatalities that would not occur under normal circumstances. Specifically in Bulgaria, 95\% of the officially confirmed COVID-19 deaths occurred in hospitals, meanng that few of those who died outside hospitals entered official statistics. The view that most excess deaths are due to COVID-19 is supported by the observation that the trajectory of excess deaths generally closely tracks that of officially recorded COVID-19 cases and deaths. Considerable discrepancies can be observed between official statistics and excess deaths, with excess deaths exceeding official numbers by even an order of magnitude or more in multiple countries\cite{Karlinsky2021}, underscoring the importance of analyzing excess mortality to accurately understand the real impact of the pandemic. During its first major wave in 2020, Bulgaria exhibited not only the highest excess mortality in the European Union, but also one of the highest discrepancies between excess deaths and official COVID deaths, with an ``undercount ratio'' of 2.52$\times$\cite{Rangachev2021,Karlinsky2021}.

We previously estimated that Bulgaria had lost 19,004 lives during its first major COVID wave in 2020. The updated analysis up to the end of March 2022 reveals that this number has increased to 68,569 (95\% CI: $\pm$6,772), compared to an official COVID death count of 36,529\cite{OWID}, i.e. the current undercount ratio is 1.88$\times$ ($\pm$0.18). In 2021, results from the most recent nationwide census for Bulgaria became available, which showed a decrease of the population down to 6,520,314\cite{NSIpreliminary}. Accounting for this updated denominator estimate, the EMR value for Bulgaria has now exceeded 1\%, standing at 1.05\% circa March 31 2022. This is the highest value recorded in any country for which excess mortality data is available\cite{Karlinsky2021}.

As crude mortality measures such as the EMR and the P-score (the percentage increase in mortality relative to baseline) may not be optimal for comparisons between populations with different demographic structures, we also calculated two standardized measures that control for such variation and aim at measuring the years of life lost as a result of the pandemic: the Potential Years of Life Lost (PYLL) and Aged-Standardized Years of life lost Rate (ASYR; see the Methods section for details). Figure \ref{Fig1} shows standardized (per 100,000 population) ASYR values for European countries in the three years of the pandemic, in total and for males and females separately. Bulgaria exhibits the highest mortality by all measures among this set of countries (per-100,000 ASYR values of 11,516, 9,157 and 13,745 in total, and for females and males, respectively), followed by Lithuania and Romania (for PYLL values see Supplementary Figure \ref{FigS7}). Excess mortality in Eastern Europe countries is much higher than that in Western Europe, and, curiously, is concentrated in the year 2021 rather than 2020, while the opposite pattern is observed in severely affected early in 2020 countries such as Spain and Italy. This observation is likely explained by two factors. First, the pandemic in 2021 in Europe was dominated first by the Alpha\cite{Davies2021} and then by the Delta\cite{Mlcochova2021} SARS-CoV-2 variants, which are known to cause more severe disease than the ancestral wild-type (WT/D614G) virus\cite{COVID192022,Murison2022,Davies2021,Challen2021,Grima2022}. Second, COVID-19 vaccination rates in Eastern Europe have been consistently lower than those in Western Europe (for example, only 11.5\% of the population in Bulgaria had received two vaccine doses by July 1st 2021, and this number only increased to 29.6\% by the end of March 2022\cite{OWID}), meaning that the Alpha, and especially the Delta waves encountered a much larger proportion of completely immunologically naive individuals in populations in Eastern Europe than in Western Europe, resulting in the observed disproportionally higher mortality in the former. Indeed, we find strong inverse correlation between vaccination rates and excess mortality, in particular in 2021 (Pearson $R^2 = 0.57$, $p \leq 0.0001$, and Spearman $r = -0.69$, $p \leq 0.0001$ for ASYR values, and Pearson $R^2 = 0.56$, $p \leq 0.0001$; Spearman $r = -0.65$, $p = 0.0001$ for PYLL; Supplementary Figure \ref{FigS4}).

\begin{FPfigure}
\begin{center}
\includegraphics[width=18.5cm]{Fig2-total-mortality.png}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Regional excess mortality patterns in Bulgaria during the COVID-19 pandemic (up to the end of March 2022)}. 
(A) Excess and official mortality by region and year;
(B) Excess mortality by region, total;
(C) Excess mortality by region, males;
(D) Bulgaria, Excess mortality by region, females;
(E) Excess mortality (P-scores) per year;
(F) Excess mortality (P-scores) for males and females per year.
}
\label{Fig2}
\end{center}
\end{FPfigure}

\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=18.5cm]{Fig8-ASYR-regions.png}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Regional excess mortality patterns in Bulgaria during the COVID-19 pandemic (up to the end of March 2022)}.
(A) Total, standardized ASYR values; 
(B) Female, standardized ASYR values; 
(C) Male, standardized ASYR values; 
}
\label{Fig8}
\end{center}
\end{figure*}

\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=18.5cm]{Fig7-CFR.png}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Regional discrepancies in the extent of recording of the impact of the COVID-19 and hospital outcomes in Bulgaria}. 
(A) Percentage of the population who have tested positive for SARS-CoV-2 in Bulgarian regions; 
(B) Undercount ratio for Bulgarian regions;
(C) Total CFR values for Bulgarian regions;
(D) In-hospital CFR values for Bulgarian regions.
}
\label{Fig7}
\end{center}
\end{figure*}

%\begin{figure*}[!ht]
%\begin{center}
%\includegraphics[width=18.5cm]{Fig4-40-64.png}
%\captionsetup{singlelinecheck=off,justification=justified}
%\caption{
%{\bf Working-age excess mortality in Bulgaria during the COVID-19 pandemic (up to the end of March 2022)}. 
%(A) Excess mortality by region, females ages 40-64, EMR values; 
%(B) Excess mortality by region, males ages 40-64, EMR values; 
%(C) Excess mortality by region, females ages 40-64, P-scores; 
%(D) Excess mortality by region, males ages 40-64, P-scores.
%}
%\label{Fig4}
%\end{center}
%\end{figure*}

 \begin{figure*}
\begin{center}
 \includegraphics[width=18.5cm]{Fig4-40-64+Fig9-WYLL.png}
 \captionsetup{singlelinecheck=off,justification=justified}
 \caption{
{\bf Working-age excess mortality in Bulgaria during the COVID-19 pandemic (up to the end of March 2022)}. 
(A) Excess mortality by region, females ages 40-64, EMR values; 
 (B) Excess mortality by region, males ages 40-64, EMR values; 
 (C) Excess mortality by region, females ages 40-64, P-scores; 
 (D) Excess mortality by region, males ages 40-64, P-scores.
 (E) Standardized WYLL values, total; 
 (F) Standardized WYLL values, female; 
 (G) Standardized WYLL values, male.
 }
 \label{Fig4}
 \end{center}
 \end{figure*}

\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=18.5cm]{Fig10-WYLL-Europe.png}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Working-age excess mortality in Europian countries during the COVID-19 pandemic (up to the end of March 2022)}. 
(A) Total, standardized WYLL values; 
(B) Female, standardized WYLL values; 
(C) Male, standardized WYLL values.
}
\label{Fig10}
\end{center}
\end{figure*}

We estimate that each excess death in Bulgaria resulted in 11.70, 12.70, and 10.43  years of life lost overall, for males, and for females, respectively, based on the ASYR metric, and in 
12.57, 12.02, and 12.51 years of life lost overall, for males, and for females, respectively, based on the PYLL metric (Supplementary Figure \ref{FigS5}).

Finally, we observe that male mortality is consistently higher than that of females for all the countries examined, consistent with previous observations\cite{Bhopal2020}.

\subsection*{Temporal trajectory of the pandemic in Bulgaria}

Figure \ref{Fig3} shows the evolution of the SARS-CoV-2 variant composition in Bulgaria based on available genome sequencing data\cite{GISAID}. The first major wave, in late 2020, was driven by WT-like (i.e. with the addition of the D614G mutation\cite{Plante2021,Yurkovetskiy2020,Korber2020} but otherwise without major spike protein mutations affecting antigenic properties) B.1.x lineages. The Alpha variant came to dominate in early 2021 and drove the second wave, and was then itself replaced by the Delta variant in June-July 2021. Finally, in early 2022, the Omicron BA.1 variant\cite{Viana2022,Cele2022}  displaced Delta and triggered the fourth major wave, with the Omicron BA.2 lineage\cite{Uraki2022,Yamasoba2022} beginning the next variant displacement cycle at the end of the observation period.

Figure \ref{Fig5} shows the trajectory of the pandemic in Bulgaria in terms of recorded clinical impacts and excess mortality. We estimate that the first wave caused $\sim$19,000 excess deaths (or EMR $\sim$0.29\% of the population), the Alpha wave had a slightly lower peak and caused $\sim$15,000 (EMR $\sim$ 0.23\%), the Delta wave peaked at about the same heights as Alpha, but was much more prolonged (Figure \ref{Fig5}A-B), and thus caused the highest number of excess deaths -- $\sim$28,000 (EMR $\sim$ 0.43). The largest number of infections were recorded during the Omicron wave (Figure \ref{Fig5}A) but it caused the fewest excess deaths, at $\sim$7,000 (EMR $\sim$0.11\%). A similar pattern is observed in the evolution of case fatality rate over time, which decreased dramatically once Omicron came to dominate (Figure \ref{Fig5}C), consistent with worldwide observations of lower disease severity with the BA.1 variant than with preceding non-Omicron ones\cite{Abdullah2022,Maslo2022}. 

Finally, we examined the ``undercount ratio'' (i.e. the ratio between excess deaths and official COVID deaths). Its values were highest, in the 2.5--3$\times$ range, during the first major wave, then decreased to the 1.5--2.4 $\times$ range during the Alpha and Delta waves, and further decreased to $\sim$1.5$\times$ during Omicron (Figure \ref{Fig5}D). The most likely in our view interpretation of these patterns is that the undercount ratio is dependent on the extent to which hospital systems are overwhelmed by surges of severe COVID-19 cases; thus the Omicron wave, which caused the fewest excess deaths, was most accurately captured in official statistics, as proportionally fewer people died outside of the hospital system, which was able to accomodate a larger share of the severe cases than in previous waves. However, even with Omicron, large unaccounted for excess mortality still persisted, likely due to the aforementioned issues of lack of testing and improper recording of causes of death. 

\subsection*{Regional mortality patterns in Bulgaria}

Next, we mapped the regional patterns of excess mortality in Bulgaria (Figures \ref{Fig2}, \ref{Fig8} and \ref{Fig7} and Supplementary Figure \ref{FigS1}). Previously\cite{Rangachev2021}, we identified a stark difference between major population centers, especially the capital Sofia, and the peripheral provinces, explained by the unfavorable demographic structure and socioeconomic characteristics of the latter (where the long-term trend has been towards depopulation, resulting in a very high median age, and an attendant decline in the availability of healthcare resources). This pattern has continued in the next three waves, and thus Sofia (city) still exhibits the lowest excess mortality in Bulgaria  (EMR = 0.67\%; Figure \ref{Fig2}A-B). In contrast, excess mortality has reached as high as 1.8\% in Vidin, 1.55\% in Montana, and 1.5\% Razgrad. Overall excess mortality is below 1\% in only five Bulgarian regions, with the northeast and northwest regions showing the highest values. 

We observe even more extreme values for sex-specific excess mortality (Figure \ref{Fig2}C-D) -- male EMR is 2.1\% in Montana and 1.95\% in Vidin. Female-specific excess mortality is considerably lower in all regions, with only five of them exceeding EMR = 1\% (it is highest in Vidin at 1.28\%).

We also examined excess mortality using the P-score metric for each year of the pandemic (Figure \ref{Fig2}E-F). This analysis confirmed the previously discussed observation of very high excess mortality centered on the year 2021, but also showed that in most regions excess mortality in the first quarter of 2022 has been comparable to that in 2020, despite the less severe phenotype of the Omicron variant. This observation is explained by the successful containment measures in the first half of 2020, contrasting with the very large number of infections in 2022.

We also analyzed regional excess mortality using the standardized ASYR  metric (Figure \ref{Fig8}; for the standardized PYLL metric see Supplementary Figure \ref{FigS6}). These comparisons revealed a somewhat different picture than crude mortality comparisons -- ASYR values are not lowest in Sofia city, and according to the ASYR  metric the northeastern provinces of Razgrad and Silistra have been more heavily affected than the northwestern ones of Vidin and Montana. This is likely because of the more extreme age skew of the demographic structure of the latter, which is normalized for by the ASYR  metric but not by crude EMR estimates. As with EMR metrics, even more extreme values are observed than the already very high one for Bulgaria as a whole -- e.g. in Razgrad the ASYR value approaches 16,000 per 100,000 population, whereas the ASYR value for Bulgaria is 11,516.

Considerable region discrepancies are also present regarding the documenting of the pandemic and the hospital outcomes for COVID-19 patients. Unfortunately, no serological survey (of any kind, not just the anti-nucleocapside protein ones that could distinguish evidence for previous infections from vaccination with mRNA or adenoviral vaccines that target only the spike protein) has ever been carried out in Bulgaria, but it is highly likely that towards the end of March 2022, a majority of the population has been infected by SARS-CoV-2 (given the observed excess mortality; to be discussed further below). However, the percentage of the population who have tested positive is highest in Sofia city, at only 14.17\%, and is as low as 4.57\% in the peripheral Kardzhali region  (Figure \ref{Fig7}A). Thus, testing has been highly inadequate throughout the pandemic, with most infections remaining undocumented. 

The undercount ratio between the EMR and the officially documented population fatality rate (PFR) ranges from 1.48$\times$ in Sliven region to 2.84$\times$ in Pernik  (Figure \ref{Fig7}B). The overall CFR ranges from 2.13\% in Sofia city to $\geq$7\% in Razgrad and Smolyan (Figure \ref{Fig7}C). These discrepancies are in large part due to the inadequate testing in some of the  peripheral regions in the country, which also tend to be the ones with the lowest percentage of the population that has tested positive. 

Remarkably, when focusing on the CFR for hospitalized patients specifically, we find no region in which fewer than 10\% of COVID-19 patients died, and in Dobrich region the number exceeds 23\%  (Figure \ref{Fig7}D), underscoring the unequal and inadequate access to high-quality COVID-19 treatment across the country.


\subsection*{COVID-19-related working-age excess mortality in Bulgaria and Europe}

Finally, we mapped the regional patterns of excess mortality for working-age populations (Figures \ref{Fig4}). We focused on the age 40-64 subpopulation as COVID-19-related deaths and excess mortality are low in absolute number in the younger demographics, resulting in statistically unreliable estimates at the regional level. 

In total, excess deaths in the 40-64 group in Bulgaria amount to 11,986 (95\% CI: $\pm$693). We find that EMR values for this group exceed 0.2\% in all regions even for females, and reach as high as 0.8\% for females and 1.03\% for males in Silistra region (Figure \ref{Fig4}A-B). Working-age excess mortality has been concentrated in the northeastern and southern regions of the country  (Figure \ref{Fig4}C-D).

We also applied a standardized analysis using the Working Years of Life Lost (WYLL) metric (see the Methods section for details), which largely confirmed these regional patterns  (Figure \ref{Fig4}E) -- in Silistra region the WYLL value exceeds 2,500 per 100,000 population, followed by Razgrad and Pazardzhik. We also note that a unique feature of regions such as Razgrad and Silitra is the very high female-specific WYLL, at nearly double that observed in other areas, and also doubling the normal death rate (p-scores nearly or exceeeding 100\%).

The average working years of life lost per excess death are 8.26 for Bulgaria as a whole, and 8.18 and 8.87 for females and males, respectively.

Finally, we compared working-age excess mortality across European countries (Figure \ref{Fig10}). Bulgaria stands out in this analysis as exhibiting standardized WYLL values far in excess of those in other countries included in the comparison ($\geq$70\% higher than the next ranked country, Romania). As in the comparison of overall excess mortality, countries in Eastern Europe exhibit considerably higher working-age excess mortality than those in Western Europe, and it is concentrated in the second year of the pandemic.

\section*{Discussion}

In this study, we map out the patterns of COVID-19-related excess mortality in Bulgaria across time, space and different demographic groups. Three striking observations stand out in the available data.

First, considerable discrepancies exist in the impact of the pandemic at the regional level, with peripheral areas of the country exhibiting much higher absolute excess mortality than the capital Sofia, presumably due to the better access to healthcare resources and the more favorable demographic structure in the latter, and possibly also the less favorable health status of the population in the former. Cardiovascular disease (CVD) is a well-known risk factor for severe COVID-19 outcomes so we examined the correlation between the CVD burden in different Bulgarian regions and excess mortality during the pandemic (Supplementary Figures \ref{FigS2} and \ref{FigS3}). We find a strong positive correlation (Pearson $R^2$ = 0.4, Spearman $r = $ 0.59) for overall excess mortality and CVD burden, and weaker such correlation (Pearson $R^2$ = 0.17, Spearman $r = $ 0.43) for male-specific age (40-64) excess mortality; these observations support such a link as one of the contributing factors (we note that we also find no such correlation for female-specific working age excess mortality and for the standardized ASYR and PYLL metrics; this is likely because CVD disease burden manifests itself earlier in males and because ASYR and PYLL place less weight on excess mortality in the very elderly where CVD burden is most pronounced). The other likely major contributing factor to regional discrepancies is the unequal distribution of healthcare resources, as we previously discussed in more detail\cite{Rangachev2021}.

Second, overall excess mortality in Bulgaria is extremely high, as it is now well in excess of 1\% of the total population. This result is very important for the overall understanding of the pandemic as it finally places the early estimates of the potential impact of the SARS-CoV-2 virus in a proper context. 

Numerous estimates for SARS-CoV-2's IFR have been published, particularly early in the pandemic. A major survey of available data\cite{COVID192022} estimated the age-standardized IFR for Bulgaria to be 0.873\% in early 2020, decreasing to 0.565\% in early 2021 (likely thanks to improved treatments). An early-2020 estimate for Belgium\cite{Molenberghs2022} placed the overall IFR at $\sim$1.5\%. Published early-2020 estimates for Spain were for IFR = 1.2\%\cite{Polan2020} and 1.15\%\cite{Kenyon2020}. For Eastern Europe as whole, an IFR value $\sim$1.45\% has been published\cite{Ghisolfi2020}. Other estimates include 0.6\% for the early-pandemic IFR in China\cite{Russell2020}, 0.5\% in Switzerland and 1.4\% in Lombardy, Italy\cite{Hauser2020} during the early-2020 wave, and meta-analysis-based overall estimates of 0.68\%\cite{Meyerowitz2020} and 1--1.5\% \cite{Thomas2021}. 

In addition, several much lower values were also published during the first year of the pandemic, such as IFR at 0.04\%\cite{Mizumoto2020}, a global one at $\sim$0.15\%\cite{Ioannidis2021}, IFR at 0.17\%\cite{Bendavid2021} for Santa Clara County in California, and others.

The validity of these estimates can be evaluated in the light of the fact that Bulgaria's excess mortality stood at 1.05\% in March 2022, and that in some regions of the country it approached 2\%. This outcome is the result of a combination of the following factors. First, a majority of the population must have been infected by that point (otherwise the IFR in Bulgaria would have to exceed 2\%, which is unlikely), although how many exactly have been infected is not possible to say in the absence of an anti-nucleocapside serosurvey (and even then, seroreversion would probaly bias estimates downwards). Second, reinfections became an increasingly common phenomenon, first with the arrival of the Delta variant\cite{Marinov2022}, and especially after the appearance of Omicron. Third, the virulence of SARS-CoV-2 prior to Omicron was increasing, with the Alpha variant being more severe than the WT and the Delta variant being even more severe than Alpha; meanwhile the IFR estimates from 2020 and early 2021 were based on the WT virus. Finally, vaccination in Bulgaria remained very low throughout the examined period, meaning that the Delta and Omicron waves were met with a large population of immunologically naive individuals, resulting in much higher mortality than in countries with high vaccination coverage. While deeply regrettable as a public health outcome for the country, this fact allows the  observation of the potential full impact of the pandemic after infecting most of a population with a high median age and in the absence of vaccination, a situation that has been avoided, at least for the time being, in most other countries with similar demographic structures.

Third, we also observe extremely high excess mortality in working age populations, far higher than that in other European countries. The EMR values in the neighborhood of 1\% in males aged 40-64 that we observe for several Bulgarian regions are around or even in excess of many of the cited above IFR values for the whole population, and well in excess of most estimates for working age demographics in particular\cite{Levin2020}. Therefore, the potential impact of the SARS-CoV-2 virus for working age people may well have been underestimated previously.

\section*{Methods}

% \subsubsection*{SARS-CoV-2 sequencing data}

% Information about sequenced SARS-CoV-2 genomes was obtained from the GISAID database\cite{GISAID}.

\subsection*{Data Sources}

All-cause mortality data for European countries and for NUTS-3 (Nomenclature of Territorial Units for Statistics) regions in Bulgaria was obtained from Eurostat\cite{EurostatCNTRY,EurostatRGN}. The data featured in these datasets is sex- and age-stratified, with age groups split in increments of 5 years. 

Country-level population data was collected through Eurostat\cite{EurostatPOP}, and was further supplemented by population data from the United Nations’ UNdata Data Service\cite{UNdataPOP}. We further elaborate on this topic in the subsequent section on Potential Years of Life Lost (PYLL) and Working Years of Life Lost (WYLL) estimates.

The preliminary data from the most recent population census in Bulgaria was used for the analysis at the regional level in the country\cite{NSIpreliminary}.

Life expectancy values at different ages were obtained from three separate sources. We acquire the full life tables for Bulgaria through the country’s National Statistical Institute\cite{NSI}. Abridged life tables for all European countries were obtained from the World Health Organization’s open data platform\cite{WHOLT}. This dataset is partitioned by age, in increments of 5 years. Abridged life tables for Bulgarian regions were created using regional mortality data for 2017--2019 collected by Bulgaria's National Statistical Institute \cite{NSI} following the methodology of the ONS\cite{ONS}.

COVID-related mortality and testing data for Bulgaria was obtained from the Bulgaria's Ministry of Health. The dataset, which covers the period from the beginning of the pandemic till March 2022 , includes information about each infected individual's age, gender, region, the date of their latest COVID-19 test, their status (infected, recovered, hospitalized, deceased), their hospitalization start and end dates, if any, whether they were taken into intensive care and whether they died of COVID-19.

\subsection*{Data Availability}

All datasets and associated code can be found at \burl{https://github.com/Mlad-en/Bulgaria_Regional_Mortality} and \burl{https://github.com/Mlad-en/COV-BG}.

\subsection*{Excess mortality and P-scores}

To calculate excess mortality across countries as well as across Bulgarian regions, we analyze the mortality observed between week 10 of 2020 and week 13 of 2022, and compare it to expected (baseline) mortality using the historical data for the five pre-pandemic years (2015--2019). The model we used is the Karlinsky--Kobak regression model\cite{Karlinsky2021}:

\begin{equation*}
D_{t,Y}=\alpha_t +\beta \cdot Y + \epsilon
\end{equation*}

Where $D_{t,Y}$ is the number of deaths observed in week or month $t$ in year $Y$, $\beta$ is a linear slope across years, and $\alpha_t$ are separate intercepts (fixed effects) for each week or month, and $\epsilon \sim \mathcal{N}(0,\sigma^2)$ is Gaussian noise. The model prediction for a year $Y$, where $Y=2020,2021 \ \text{or} \ 2022$ is $\mathrm{Expected \ Mortality}_{t,Y}=\hat{\alpha_t}+\hat{\beta}\cdot Y.$

We then establish a 95\% confidence interval for the expected mortality. This range is used to calculate the excess mortality $\Delta_t$ for a week or a month $t$ and a year $Y$ as:

%We derive this mean both on a weekly and aggregate total basis. 
\begin{equation*}
 \Delta_{t,Y}= \mbox{Mortality}_{t,Y} - \mbox{Expected Mortality}_{t,Y}.
% $\mathrm{Excess \ Mortality} = \mathrm{Mortality}_{2020} - \mathrm{Mean \ Mortality}$
\end{equation*}

This calculation is done both as a sex- and age-stratified metric, as well as an aggregated total excess mortality for a year $Y$, which we denote by $\Delta_Y$. To normalize excess mortality across countries, we calculate excess mortality per total population. To do this, we use population data from Eurostat for 2020.

Set $z_Y:= |\Delta_Y|/\sqrt{\mathrm{Var}[\Delta_Y]}$, where $\mathrm{Var}[\Delta_Y]$ is computed in \cite{Karlinsky2021}. If $z_Y$ is significantly below $2$ for a given country, we consider the excess mortality for this country to be not significantly different from zero. In the computations related to the years-of-life lost metrics considered in the paper, we excluded a few countries having both $z_Y$-values significantly below $2$ (typically less than $1$) for each age interval and wide confidence intervals that included $0$ for the excess mortality associated with each of these age intervals.

Based on the excess mortality ranges we also compute a P-score value for each country/region. A P-score value is defined as the ratio or percentage of excess deaths over certain period relative to the expected deaths for the same period based on historical data from the years 2015--2019 (see \cite{OURWORLD}).
We calculate the P-score for a year $Y$ as follows: 

$$P_Y:=\frac{\mathrm{Mortality_{Y}} - \mathrm{Expected \ Mortality}_{Y}}  {\mathrm{Expected \ Mortality}_{Y}} * 100$$

To calculate a total $P$-score we replace each term in  in the right-hand side in the formula above by the corresponding summation over the three-year period considered in our analysis. We also calculate the ratio between excess mortality and official COVID-19-attributed mortality. Due to the demonstrably low testing in Bulgaria\cite{ref5} and other countries, this allows us to estimate under-reported COVID-19 fatalities. We also use the total positive tests per region to compute a Case Fatality Ratio (CFR) which estimates the proportion of COVID-19 fatalities among confirmed cases.

\subsection*{Potential Years of Life Lost (PYLL), Aged-Standardized Years of life lost Rate (ASYR), and Working Years of Life Lost (WYLL) estimates}

Potential Years of Life Lost (PYLL) is a metric that estimates the burden of disease on a given population by looking at premature mortality. It is derived as the difference between a person's age at the time died and the expected years of life for people at that age in a given country. As such, the metric attributes more weight to people that have died at a younger age. 

We compute the PYLL across countries and Bulgarian regions by taking the positive all-cause excess mortality for all ages groups (in Eurostat they are aggregated at 5 year intervals). For the European countries considered in our paper we use the abridged life expectancy tables by the WHO (also aggregated at 5 year intervals) and for the Bulgarian regions we create abridged life expectancy tables following the ONS methodology\cite{ONS} to calculate a total and average PYLL value for all countries and Bulgarian regions. To be more precise, for an age interval $[x,x+4]$ and sex $s$ (if no sex is specified we assume it's for both sexes) define by $\mathrm{ED}([x,x+4],s)$ the excess deaths and by $\mathrm{LE}([x,x+4],s)$ the life expectancy. Then the potential years of life lost are computed as $$\mathrm{PYLL}([x,x+4],s)=\mathrm{ED}([x,x+4],s)*\mathrm{LE}([x,x+4],s).$$
The total PYLL is computed by summing over all age intervals. In our computations we take into account the margin of error for each $\mathrm{ED}([x,x+4],s)$.

A limitation on this approach is the upper-boundary aggregation value for the two datasets. The all-cause mortality dataset's upper boundary is 90+, while the WHO's abridged life tables only go up to the 85+ age bracket. To account for this, we attribute the life expectancy of the 85+ age group to the 85-89 mortality group. We have further excluded the 90+ mortality group from our analysis. % This is further elaborated on in the Limitations subsection, where we also provide a way of correcting for this exclusion.

Finally, we standardize PYLL values across countries by dividing the total sum value by the population and normalizing it per 100,000 people: 

$$\mathrm{PYLL_{std}}:=\frac{\mathrm{PYLL_{total}}}  {\mathrm{Total \ Country \ Population_{0-89}}} * 100,000$$

The data for country-level populations in Eurostat has a similar limitation in the upper boundary of the age distribution (a cut-off at 85+). To mitigate this limitation, we supplement the population data from Eurostat for ages 0-84 with population size data for the 85-89 age group from the UNdata Data Service. 

To compare the impact of the pandemic across European populations and Bulgarian regions with different age structures we compute the Age-Standardized Years of Life Lost Rate (ASYR)\cite{DBurden,Asyr}. Let $([x,x+4],s)$ be an age interval for a sex $s$ in a standard life expectancy table for a given population. Denote by $P([x,x+4],s)$ the population size of $([x,x+4],s)$. Define the PYLL rate for $([x,x+4],s)$ as $$\mathrm{PYLL_{rate}}([x,x+4],s):=\frac{\mathrm{PYLL}([x,x+4])}{P([x,x+4],s)}*100,000.$$

For the 2013 European Standard Population (ESP) denote by $W([x,x+4],s)$ the weight of $([x,x+4],s)$ in the standard population. Define $$\mathrm{ASYR}(s):=\sum \mathrm{PYLL_{rate}}([x,x+4],s) * W([x,x+4],s)$$ where the sum is taken over all age intervals. For a given population of sex $s$ this measure is interpreted as the years of life lost per $100,000$ people (of sex $s$) if the population has the same age distribution as the ESP. We do the same for the Bulgarian regions using a standardized population for Bulgaria based on $2019$ census estimates by the Bulgarian NSI. $\mathrm{ASYR}$ allows for comparison of the pandemic impact on EU countries and Bulgarian regions having different age distributions. 

Finally, we derive total, average and total standardized WYLL value approximations. To accomplish this, we first assume people to be in the working age group if they are $15$ to $64$ years old, and thus exclude excess mortality for all age groups over $65$. To calculate the remaining years of working life, we further assume a mean age for each age group, e.g. for the age interval $60-64$ we assume a mean age at $62.5$ years. This would leave this group with approximately $2.5$ years until retirement. $95\%$ CI for EMRs, P-scores and the values of all years of life lost functions can be found in our GitHub repository.
% Limitations on this approach are discussed in the subsequent Limitations subsection.

% \subsection*{Stringency index and Mobility data}

% \hl{XXX}

% \subsection*{XXX Hospitalization statistics XXX}

% \hl{XXX}

\subsection*{Limitations}

Each of the presented data sources and approaches to analysis have their own limitations. Below we discuss each one in detail.

\subsubsection*{Limitation of Excess mortality measures}

Influenza outbreaks in the period 2015--2019 contribute to the estimation for the expected mortality for 2020--2022. Thus the expected mortality is an estimate of the ``normal" death rate in the presence of seasonal influenza.
%Our mean mortality calculations do not account for population changes in the 2015-2019 period. For example, Bulgaria has one of the highest negative growth rates in the world. Between 2015 and 2020 Bulgaria's population has decreased by around 250,000 or approximately 6.5\% per year\cite{ref1}. The country also has a negative net migration and has seen approximately 25,000 people leave its borders in the last five years\cite{ref2}
% Compare this to the relatively net neutral population change in Czechia for the last five years or the positive population growth of Iceland with approximately 5.5\% percent per year.

\subsubsection*{Limitations of PYLL/ASYR/WYLL}

Since PYLL, ASYR and WYLL data only take into account fatalities, these metrics do not provide information about any worsened quality of life of surviving individuals, reduced life expectancy of these individuals and working capacity. Metrics such as Disability-Adjusted Life Years (DALY), Quality-adjusted life year (QALY) and Healthy Years of Life (HALE) metrics may illuminate further the total disease burden on the European population, however, obtaining the necessary information for these measurements is not yet possible. % for all countries analyzed.

% \hl{NOT SURE IF THIS PARAGRAPH MAKES SENSE}
% Furthermore, neither metric takes into account the potential comorbidities of the people that have died. \hl{Information about reduced life expectancy of people with diabetes/hypertension/etc.}
% Such research, however, is a massive undertaking and would require a lot more data being made available.

As mentioned before, due to data availability limitations from Eurostat in our computations of PYLLs and ASYRs we excluded the $90+$ group. Given that countries like France, Italy and Spain have significant excess mortality in this age group, we also present a computation of the ASYRs including the $90+$ age by assuming $4$ years of life expectancy (the average life expectancy for the $90+$ age group for the European population is $4.74$, according to the UNdata Data Service)% in Supplementary Tables \ref{TableS16} and \ref{TableS17}
. By linear interpolation, $4$ is approximately the life expectancy for the age interval $[90,94]$ assuming that the average age of deaths for this interval is in the range $92.7-93$ (the average age of the COVID-19 fatalities above 90 years of age is $92.7$ in Czechia and $92.1$ in Bulgaria, and likely it's higher in Western countries which overall have higher life expectancy). 


This rough approximation gives an upper bound of how large the ASYRs can go. It leads to $5\%-14\%$ and $14\%-22\%$ increase in the ASYRs for the $(0-89)$ population of Eastern and Western European countries, respectively, but it does not yield a decrease between the inequalities of the countries from the two groups or any significant change in their ranks (see Supplementary Fig.\ \ref{FigS5}% and Supplementary Table \ref{TableS17}
).

The WYLL measure we present has some additional limitations. The first comes from the assumption that retirement age across European countries is 65. While it is most often assumed as a standard between European countries, there is actually some variation between individual member states\cite{ref3}. Furthermore, we assume that the mean age of people who have died in a given age group is the middle of the given range, e.g. for the age group 60-64 - mean age = 62.5. It may well be a fact that a majority of the fatalities are concentrated in the upper part of the age bracket. However, since we do not have data about the different causes of mortality, but rather an aggregate total, we cannot be certain that this trend will hold true for all age groups and across different countries.


\section*{Notes}

\subsection*{IRB and/or ethics committee approval}

The Ethics Committee of IMI-BAS (Institute of Mathematics and Informatics, Bulgarian Academy of Sciences) gave ethical approval for this work.

\subsection*{Competing Interests} 

The authors declare no competing interests.

\subsection*{Author contributions}

A.R, G.K.M. and M.M. conceived the project, analyzed the data and wrote the manuscript.

\subsection*{Acknowledgments and Funding}
The authors would like to acknowledge the help of the Bulgarian Ministry of Health and Information Services for providing us with raw data about COVID-19 mortality, hospitalizations and testing. A.R. would like to acknowledge the financial support of a ``Petar Beron i NIE" fellowship [KP-06-D15-1] from the Bulgarian Science Fund. 

\begin{thebibliography}{100}

% \section*{References}

\input{references}

\end{thebibliography}

\clearpage

\setcounter{table}{0}
\renewcommand{\tablename}{Supplementary Table}
\setcounter{figure}{0}
\renewcommand{\figurename}{Supplementary Figure}

\setcounter{page}{1}
\renewcommand\thepage{{SM }\arabic{page}}

\begin{center}
% {\LARGE \textbf{\begin{spacing}{1.1}XXXX. \\ Supplementary Materials\end{spacing} }}
{\LARGE \textbf{Supplementary Materials}}
\end{center}

\section*{Supplementary Figures}\label{SM}

\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=18.5cm]{FigS7-PYLL-Europe.png}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Excess mortality in Europe and Bulgaria during the COVID-19 pandemic (up to the end of March 2022)}. 
(A) Standardized PYLL values, total; 
(B) Standardized PYLL values, females; 
(C) Standardized PYLL values, males.
}
\label{FigS7}
\end{center}
\end{figure*}

\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=18.5cm]{FigS4-ASYR-PYLL-vs-vaccination.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Correlation between rates of full vaccination and ASYR and PYLL excess mortality measures in European countries}. \\
(A) ASYR, total (Pearson $R^2 = 0.5$, $p \leq 0.0001$; Spearman $r = -0.64$, $p = 0.0002$); \\
(B) ASYR, 2020 (Pearson $R^2 = 0.22$, $p = 0.0105$; Spearman $r = -0.36$, $p = 0.0549$); \\
(C) ASYR, 2021 (Pearson $R^2 = 0.57$, $p \leq 0.0001$; Spearman $r = -0.69$, $p \leq 0.0001$); \\
(D) ASYR, 2022 (Pearson $R^2 = 0.08$, $p = 0.1286$; Spearman $r = -0.2$, $p = 0.2982$); \\
(E) PYLL, total (Pearson $R^2 = 0.47$, $p \leq 0.0001$; Spearman $r = -0.62$, $p = 0.0003$); \\
(F) PYLL, 2020 (Pearson $R^2 = 0.19$, $p = 0.0169$; Spearman $r = -0.32$, $p = 0.08$); \\
(G) PYLL, 2021 (Pearson $R^2 = 0.56$, $p \leq 0.0001$; Spearman $r = -0.65$, $p = 0.0001$); \\
(H) PYLL, 2022 (Pearson $R^2 = 0.13$, $p = 0.0523$; Spearman $r = -0.19$, $p = 0.3106$).
}
\label{FigS4}
\end{figure*}

\clearpage

\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=15cm]{FigS5-PYLL-ASYR-per-death-Europe.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Year of life lost per excess death in European countries}. Note that the very high values for a few countries (e.g. Cyprus, Iceland, Luxembourg, Finland, Malta) might be artifacts resulting from significantly insignificant excess deaths (z-score significantly below $2$).
}
\label{FigS5}
\end{figure*}

\clearpage

\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=18.5cm]{FigS1-total-p-scores.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Excess mortality in Bulgarian regions (P-scores)}. 
}
\label{FigS1}
\end{figure*}

\clearpage

\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=18.5cm]{FigS6-PYLL-regions.png}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Regional excess mortality patterns in Bulgaria during the COVID-19 pandemic (up to the end of March 2022)}.
(A) Total, standardized PYLL values; 
(B) Female, standardized PYLL values; 
(C) Male, standardized PYLL values.
}
\label{FigS6}
\end{center}
\end{figure*}

\clearpage


\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=18.5cm]{FigS2-CVD-vs-EMR.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Correlation between cardiovascular disease (CVD) prevalence (measured as the deaths from CVD per 100,000 people in 2019) and COVID-related excess mortality in Bulgarian regions (as measured by EMR)}. Shown is the Pearson $R^2$ correlation coefficient. \\
(A) total (Pearson $R^2 = 0.4$, $p = 0.0003$; Spearman $r =  0.59$, $p = 0.0008$); \\
(B) males (Pearson $R^2 = 0.35$, $p = 0.0009$; Spearman $r =  0.57$, $p = 0.0014$); \\
(C) females (Pearson $R^2 = 0.21$, $p = 0.0129$; Spearman $r =  0.43$, $p = 0.0196$); \\
(D) males ages 40-64 (Pearson $R^2 = 0.17$, $p = 0.026$; Spearman $r = 0.43$, $p = 0.0214$); \\
(E) females ages 40-64 (Pearson $R^2 = 0.02$, $p = 0.41$; Spearman $r = 0.12$, $p = 0.51$). 
}
\label{FigS2}
\end{figure*}

\clearpage

\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=18.5cm]{FigS3-CVD-vs-ASYR.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{
{\bf Correlation between cardiovascular disease (CVD) prevalence (measured as the deaths from CVD per 100,000 people in 2019) and COVID-related excess mortality in Bulgarian regions (as measured by ASYR)}. Shown is the Pearson $R^2$ correlation coefficient. \\
(A) total (Pearson $R^2 = 0.02$, $p = 0.47$; Spearman $r =  0.18$, $p = 0.35$); \\
(B) males (Pearson $R^2 = 0.01$, $p = 0.31$; Spearman $r =  0.15$, $p = 0.44$); \\
(C) females (Pearson $R^2 = 0.03$, $p = 0.79$; Spearman $r =  -0.16$, $p = 0.40$).
}
\label{FigS3}
\end{figure*}

% \begin{figure*}[!ht]
% \begin{center}
% \includegraphics[width=18.5cm]{FigS7-PYLL-Europe.png}
% \captionsetup{singlelinecheck=off,justification=justified}
% \caption{
% {\bf Excess mortality in Europe and Bulgaria during the COVID-19 pandemic (up to the end of March 2022)}. 
% (A) Standardized PYLL values, total; 
% (B) Standardized PYLL values, females; 
% (C) Standardized PYLL values, males.
% }
% \label{FigS7}
% \end{center}
% \end{figure*}

% \begin{figure*}
% \begin{center}
% \includegraphics[width=18.5cm]{FigS6-PYLL-regions.png}
% \captionsetup{singlelinecheck=off,justification=justified}
% \caption{
% {\bf Regional excess mortality patterns in Bulgaria during the COVID-19 pandemic (up to the end of March 2022)}.
% (A) Total, standardized PYLL values; 
% (B) Female, standardized PYLL values; 
% (C) Male, standardized PYLL values.
% }
% \label{FigS6}
% \end{center}
% \end{figure*}


% \section*{Supplementary Tables}\label{ST}

% \input{TableS1-Europe-excess-mortality-total}

\end{document}