Socio-demographic and health factors drive the epidemic progression and should guide vaccination strategies for best COVID-19 containment.

We propose and study an epidemiological model on a social network that takes into account heterogeneity of the population and different vaccination strategies. In particular, we study how the COVID-19 epidemics evolves and how it is contained by different vaccination scenarios by taking into account data showing that older people, as well as individuals with comorbidities and poor metabolic health, and people coming from economically depressed areas with lower quality of life in general, are more likely to develop severe COVID-19 symptoms, and quicker loss of immunity and are therefore more prone to reinfection. Our results reveal that the structure and the spatial arrangement of subpopulations are important epidemiological determinants. In a healthier society the disease spreads more rapidly but the consequences are less disastrous as in a society with more prevalent chronic comorbidities. If individuals with poor health are segregated within one community, the epidemic outcome is less favorable. Moreover, we show that, contrary to currently widely adopted vaccination policies, prioritizing elderly and other higher-risk groups is beneficial only if the supply of vaccine is high. If, however, the vaccination availability is limited, and if the demographic distribution across the social network is homogeneous, better epidemic outcomes are achieved if healthy people are vaccinated first. Only when higher-risk groups are segregated, like in elderly homes, their prioritization will lead to lower COVID-19 related deaths. Accordingly, young and healthy individuals should view vaccine uptake as not only protecting them, but perhaps even more so protecting the more vulnerable socio-demographic groups.

Introduction

Currently we are experiencing a global COVID-19 epidemic affecting almost all countries in the world. The high infection rate and relatively long hospitalization of risk groups presents a great concern for governments, and a challenge for researchers trying to prevent an economic, social and political crisis [1], [2]. This has stimulated an unprecedented interest in epidemiological models, especially predicting the outcomes of scenarios considering counter- and prevention measures, in particular the forthcoming vaccination. Complex systems and network science approaches, along with technological advances and data availability, are becoming instrumental for the design of effective containment strategies [3], [4], [5], [6], [7]. Numerous recent works are devoted to fitting of the available data [8], [9], inferring the key epidemiological processes [10], [11], [12], identifying the control knobs for non-pharmaceutical interventions [13], [14], aiding decisions on emergency management [15], [16], checking the effectiveness and importance of lockdowns [17] and making predictions about the further epidemic progression [18], [19], [20], with the aim to support the decision-making process amid the crisis.

For COVID-19 it is widely believed that the society will not be able to return to normal pre-pandemic life until a global vaccination program is successfully implemented [21]. Finding efficient vaccine administration strategies is crucial, not only because the pandemic has generated simultaneous demand all around the globe, which rises technical as well as moral considerations in decision making, but also due to anti-vaccination movements and conspiracy theories that have the potential to reduce vaccine uptake [22], [23]. Moreover, the epidemic progression depends largely on the heterogeneity of the population [24], [25]. Firstly, different age and social groups can have different frequencies of interactions, which has been shown to substantially affect epidemiological parameters of the new coronavirus disease [26], [27]. Secondly, infected individuals express highly heterogeneous clinical and immunological manifestations of COVID-19 and high variability is also expected in responses to vaccination [28]. For these reasons, assessing the pandemic evolution and the corresponding ongoing immunization processes should also take notice of the broader socio-economic, cultural as well as immunological and other biological factors.

The nature of protective immune responses is one of the main questions in COVID-19 research [21]. Besides the age, it is not clear how different sub-populations respond to natural infection or vaccination. One of the biggest concerns nowadays is how metabolic disorders and other related diseases may hinder COVID-19 vaccine efficacy [29]. Previous studies for other viral diseases, e.g., influenza, have shown a reduced effectiveness of vaccinations in people with obesity [30], [31]. It has been shown that obese individuals 12 months after vaccination had a decreased CD8 + T cell activity and a more pronounced decline in influenza antibody titers [32]. For individuals with metabolic-related diseases, a shorter post-vaccination immunity has also been confirmed for hepatitis B and tetanus [33]. Therefore, different types of comorbidities affecting perturbations in the immune system are also likely to have implications for COVID-19 vaccination [31].

In the last year, a significant number of different types of models have been designed to simulate and predict the trajectories of the COVID-19 pandemic and to explore how herd immunity will contain the pandemic once a sufficiently high proportion of the population acquires immunity [34], [35]. Mathematical models have been used to explore the possible epidemic progressions with regard to variations in immune responses following SARS-CoV-2 infection and vaccination [25], [36] and to evaluate and derive optimal prioritization strategies of vaccine administration [37], [38], [39]. It has been argued that models of COVID-19 transmission should aim to capture the complex interactions between individuals in the heterogeneous and geographically dispersed populations. The reason being that these factors, along with other behavioral interventions, importantly contribute to the dynamics of the disease [40], [41]. Moreover, considering that decision making will have to take place under uncertainty and imperfect information, and with only conditionally optimal outcomes, game theory and social network models could be a valuable repertoire to guide decisions pertaining to vaccination programs [42]. In studies of disease transmission in heterogeneous populations, where individual characteristics or characteristics of subgroups may play a decisive role, network models and agent-based modeling formulations have already proven beneficial [43], [44], [45], [46], [47], including in simulations of vaccination [36], [48], [49], [50]. These models represent a clear departure from the homogeneity assumptions of the traditional well-mixed compartment models and have already been used in the context of the COVID-19 epidemic. Specifically, it has been shown that utilizing accurate social network models can provide better insights into the spread of the new coronavirus disease [51], [52], [53], [54], improve the forecast of the evolution of the epidemic [55], [56] and help to design efficient strategies to control COVID-19 outbreaks [57].

In this work, we explore how heterogeneity in health status along with spatial distribution patterns within the population influences the progression of the epidemic. Subpopulations with different grades of comorbidities have different responses to virus infection and to vaccination. In subpopulations with a lower health status, it is expected that the probability for a severe form of the disease is higher and the immune periods are shorter, which might considerably affect the spread of COVID-19 in a network of interacting individuals. To address these issues, we utilize a stochastic multi-compartment epidemic model and analyze how a variable transient immunity in a heterogeneous population affects the dynamics and outcome of the epidemic.

Computational model

We have developed an agent-based computational model to simulate the spread of the new coronavirus disease and vaccination in a heterogeneous population. In the model we consider that individuals have diverse medical conditions by arranging them into three subgroups. We utilize a stochastic extended SEIRS formalism and explore how these aspects affect the progression and outcome of the COVID-19 epidemic. Our framework incorporates the epidemiological compartments characteristic for COVID-19 and in which individuals interact via a realistic social network scheme. The outline of the model is schematically presented in Fig. 1 , whereas a detailed description of individual segments of the model is given in continuation. Moreover, we have performed a detailed literature overview to estimate the implemented model parameters on published data, which are provided in the Supplementary text S1.

Fig. 1

Setup of the multi-compartmental epidemiological model. (A) Distribution of underlying medical conditions within a population (green, healthy (RG1); orange, moderate underlying medical conditions (RG2); red, severe chronic medical conditions (RG3)). (B) Individuals within the population are represented as nodes, and social interactions are represented as edges between the nodes. (C) The scheme of the extended SEIRS epidemiological model. Individual model compartments stand for susceptible to infection (S), exposed to the infection (E), presymptomatic state (P), an asymptomatic form of the disease (IA), symptomatic form of the disease (IS), hospitalized (H), recovered (R), deceased (D), vaccinated (V) and vaccinated individual who developed antibodies (Im). See main text and Tables S1 and S2 in the Supplementary text S1 for a detailed description of the compartments, residence times, transmission probabilities, and other parameters. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Heterogeneity of the population

In our simulations, N individuals are initially distributed into three categories, characterizing the individual’s health status. Fractions p RG1, p RG2, and p RG3 of individuals are assigned to the categories risk group 1 (RG1), risk group 2 (RG2), and risk group 3 (RG3), respectively, so that p RG1 + p RG2 + p RG3 = 1 (see Fig. 1A). In the subcategory RG1, we consider healthy Individuals with no underlying medical condition. Individuals assigned to RG2 have moderate underlying medical conditions, whereas individuals assigned to RG3 have full-blown metabolic disorders and other related severe diseases. These risk groups are used to model COVID-19 progression specific to the corresponding underlying medical conditions. Individuals from RG2 and RG3 exhibit a higher risk of developing severe symptoms, they are more likely to get hospitalized and their immune periods are shorter. Their specific parameters in terms of transition probabilities between compartments, compartments specific residence durations, and length of immunity after vaccination and recovery from COVID-19 are given in Table S3. It should be noted that in our model, we do not explicitly consider different age groups. However, the distribution of individuals within the categories can be regarded to reflect the population's age structure, but a lower health status can, in general, also be a consequence of the presence of disadvantaged socio-economic groups.

Social network model

To simulate the human interaction patterns, we utilize the random geometric graphs in hyperbolic spaces. This network model produces interaction structures with genuine characteristics of social networks [58], [59], [60]. Networks embedded in hyperbolic spaces exhibit a small diameter, strong clustering, community structure, and a heterogeneous degree distribution and are therefore commonly used to simulate various social phenomena, including epidemics [61], [62], [63]. In our study, we utilized the geometric preferential attachment model, where each newly added node i is mapped into the hyperbolic disc with randomly assigned polar coordinates: where, α is the internal growth parameter, R hd = 1 is the radius of the hyperbolic disc, and u 1 and u 2 are independent random variables sampled from the uniform distribution on a unit interval. The new i-th node connects with n = 5 existing nodes with a probability that is proportional to the distance between the i-th and the j-th node:where is the angular distance [58]. Notably, the polar coordinate r reflects individuals' popularity, as the nodes at the disc's periphery exhibit fewer connections than the nodes that are close to the center. Moreover, nodes separated by a small angular distance correspond to individuals that can be regarded to live and operate in the same geographic region or belong to the same geo-political or socio-economic group [64]. We used parameters , , , and nodes, which yielded realistic social network architectures.

Arrangement of individuals for different risk groups within the network

In our study, we additionally investigate how the epidemic trajectories are affected by the spatial distribution of the risk groups within the population. We consider two extreme scenarios. In the first case, we homogeneously distribute the individuals from all 3 risk groups throughout the whole network. In the second case, we spatially segregate individuals belonging to RG1 from individuals belonging to RG2 and RG3, so that they are all in one district. The distributions of individuals with lower health status for both types of populations and how they are arranged within the social network model are visualized in Fig. 2 .

Fig. 2

Schematic representation of risk group arrangement within the social network model. Scheme of social interactions based on a homogeneous risk group distribution (A), the probability distribution of different risk groups across the population (B), social interaction network where risk groups are in different districts, i.e., a segregated population structure (C), and probability of risk group share in the individual compartments of the network (D). Nodes in the network are color-coded based on their corresponding risk group (green – risk group 1, orange – risk group 2 and red – risk group 3). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Network-based multi-compartment epidemiological model

Infection-disease spreading in a network is modeled with the extended SEIRS compartment model. In addition to the standard compartments susceptible (S), exposed (E), infected (I), and recovered (R), we incorporated the following compartments: presymptomatic (P), hospitalized (H), vaccinated (V), vaccine-induced immunity (I m) and deceased (D). The compartment labeled as infected was divided into two separate compartments, namely the asymptomatically infected (I A) and symptomatically infected (I S), as schematically presented in Fig. 1C. It should be noted that for practical reasons only nodes which will certainly become infected enter the state transient E, whereas nodes who were in contact with an infectious node and did not become infected, remain formally in state S. In the simulation, each node of the network can be in one of these exclusive states. Nodes interact in accordance with the underlying network structure and once a susceptible node interacts with an infected node, it is transferred to state E with a given probability. Afterwards, the state of the nodes changes with time in accordance with the model assumptions, which are explained in more detail in continuation. Description of the transition probabilities and residence times, separately for different population groups, along with the corresponding references, is provided in the Supplementary text S1. It should be noted that the data from the literature that we used to calibrate our model exhibits large variations. As pointed out by Català et al. [65], any analysis based on diagnosed cases is biased by diagnosis protocols. Moreover, several other factors contribute to the observed variation in country specific COVID-19 statistics, like for example the pool of asymptomatic cases and the health system structure, which should be kept in mind when designing COVID-19 epidemic models.

Finally, we incorporate restrictions q, which simulate the effect of non-pharmaceutical interventions to control virus spreading and affect infection rates as follows: where , and are the probabilities of a susceptible individual becoming infected from a presymptomatic, asymptomatic, or symptomatic individual when no restrictions are applied. These simplified terms, given in Eqs. (4–7) mimic both external (containment policies) and internal (endogenous) social distancing principles and contact reductions, which are a function of symptomatic infected individuals and hospitalized individuals. To consider real-life implications of restrictions governed by those two variables, we consider the 14-day average of the share of symptomatic infected individuals and share of hospitalized individuals , whereby the share of hospitalized individuals is given a larger weight (see Table S3). We considered two intensities of prevention measurements in our simulations: one strict (Q50) and the other more permissive (Q25).

Monte Carlo simulations

We perform Monte Carlo simulations to simulate the spreading of the COVID-19 through a complex network of interacting individuals, some of which have different grades of comorbidities. In accordance with the underlying health status category, we assign every node a specific set of parameters that govern the disease's progression. We initially select 0.1% of nodes at random in all our simulations and designate them as infected/presymptomatic (P), whereas all other nodes are designated as susceptible (S). In each full Monte Carlo step (MCS), the following elementary step is repeated N times, with N being the size of the population. A node i is selected uniformly at random from the whole network, and if the selected node i is in state S, we randomly choose one of its neighbors j. If node j is in an infectious state (compartments P, I A, or I S), node i is flipped to state E with probabilities p P, p A, or p S, depending on whether node j is in states P, I A or I S. If, however, the neighbor j is in any of the other states, nothing happens.

Newly infected nodes are transferred to the exposed compartment (E). In this latent state, nodes are not treated as contagious and cannot spread the infection within the population. After T E full MCS the node is shifted to the presymptomatic compartment P, and it can infect other nodes of the network. From there, nodes are transferred either to the asymptomatic compartment (I A) with a probability p PA , or the symptomatic compartment (I S) with a probability p PS, , All nodes in compartment I A, remain there for T A full MCS and are relocated to the recovered compartment (R) afterwards. Nodes transferred to the compartment I S, stay there for T S, full MCS. Afterward, nodes are transferred either to the hospitalized compartment (H) with a probability p SH, or to compartment R with a probability p SR, . Nodes in compartment H, remain there for T H,i full MCS. From this compartment, nodes are moved either to the compartment R with a probability p HR, , or to the deceased compartment (D). It should be noted that we did not incorporated natural birth and death processes in our model, as the temporal scale at which COVID-19 evolves is much faster [66]. Moreover, the length of our simulations and the total mortality rates are rather low, so that the final number of individuals in compartment D remains much lower compared to the size of the whole population. Deceased individuals therefore do not noticeably affect the spreading dynamics, even though they do not further interact with other nodes and are not replaced.

Nodes in the compartment R are treated as immune, but after a transient period T R, they switch back to the susceptible compartment S. The residence times T P, H, and T R, are drawn from a log-normal distribution, whereas other residence times are considered constant. It should be noted that transition probabilities p PA, , p SR, , p SH, , p HD, , and p VA, as well as the residence times T R, and T H, depend on the health status of individuals. The former are considered constant for nodes of the same risk group and are provided in Supplementary text S1 (Table S3), the latter, on the other hand, are drawn from log-normal distributions described in Supplementary text S1 (Table S2) and weighted for different risk groups according to parameters f Im, and f H, described in Supplementary text S1 (Table S3), respectively.

In parallel, nodes in the susceptible state are being vaccinated with a given vaccination rate p V, whereby a vaccinated individual remains in the susceptible compartment for T V full MCS. Afterward, the nodes switch to the immunity compartment I m with a probability p VA, , while otherwise, they remain susceptible. Vaccinated individuals who develop immunity remain in the compartment I m for T Im, full MCS. This immunity period is drawn from a log-normal distribution described in Supplementary text S1 (Table S2), and the final value is weighted in accordance with the node's health status with parameter f Im, (see Supplementary text S1, Table S3).

Results

We developed an extended SEIRS network model to simulate the dynamics of the SARS-CoV-2 epidemics in a heterogeneous population subjected to infected-dependent restrictions in social interactions and to explore how different vaccination strategies contribute to the confinement of disease progression. In our model, we incorporated the presence of subpopulations with different grades of comorbidities. We integrated different simulation scenarios, where individuals with these systemic diseases were either distributed homogeneously among the population or all located in one district, as presented in Fig. 2. We investigated the evolution of the epidemics in two types of populations. The first one resembles a younger and healthier society (MHP) with a low fraction of individuals with poor medical conditions, and the second one reflects an older society (LHP) with a relatively high fraction of people with systemic metabolic disturbances and associated comorbidities. Besides, we considered two strengths of enforced prevention measurements, as defined by Eq. (4). One strict (Q50) and the other more permissive (Q25). Initially, we focused on simulating the evolution of the COVID-19 epidemic in the heterogeneous network without the introduction of vaccination. The aim was to determine how variability in metabolic health affects the first 9 months of epidemic progression with regard to the structure of the society (homogeneous and segregated) and disease containment strategies (Q25 and Q50). In continuation, the state of the system after 9 months was the initial state for additional simulations, wherein the future epidemic trajectories were studied by incorporating vaccination. Special attention was given to the interplay between vaccination strategies and the heterogeneity of the population.

Simulation of the progression of COVID-19 epidemics without vaccination

We start with showing the time curves of the simulations from the extended SEIRS model without vaccination for the first 9 months of the epidemic. In Fig. 3 A-F, we show the temporal evolution of the relative number of infected, hospitalized, and the cumulative share of deceased individuals. The simulations were performed separately for the homogeneous (Fig. 3A-C) and segregated (Fig. 3D-F) population distributions within the network, for both types of societies (blue and brown lines), and two levels of non-pharmaceutical intervention measures (full and dotted lines). In all scenarios, the number of infected individuals exhibits an initial peak, followed by a rather steady and prolonged plateau phase. The curve showing the current number of hospitalizations has a very similar shape but is delayed by 15–20 days. As expected, the corresponding cumulative share of deaths increases the fastest shortly after the initial peak of infections and stabilizes to a nearly linear trend afterward in the plateau phase. In both types of societies, i.e., MHP and LHP, we can observe that all the observed measures (infected, hospitalized, and deceased) were lower in the case of more intense restrictions. Moreover, our results show that in the overall healthier society (MHP), the share of infected individuals is higher than in the society with more abundant chronic health problems (LHP). However, in the LHP society, a higher share of hospitalization and deaths can be observed regardless of the intensity of the enforced restriction.

Fig. 3

Simulating the first 280 days of disease dynamics in heterogeneous populations. Time series of actively infected individuals (A, D), the share of hospitalized individuals (B, E), and share of cumulatively deceased individuals (C, F) when risk groups are homogeneously distributed throughout the social network (A-C) and when individuals with lower health status are located in one district, i.e., a segregated structure (D-F). The temporal evolution of the epidemic is presented by the mean value and the 95% confidence interval around the mean value determined based on 50 independent simulation runs. Cumulative fractions of infections (G), hospitalizations (H), and deaths (I) on the 280th day for both population types (MHP, blue; LHP, orange), restriction strengths (Q25 and Q50), and homogeneous and segregated risk group distribution within the social network. Error bars denote standard deviations. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

To evaluate these differences quantitatively, we present in Fig. 3G-I the overall fractions of infected, hospitalized, and deceased individuals. It should be noted that the overall fraction represents the share of individuals who got infected, were hospitalized, or died in the first 9 months of the simulation. It can be observed that a higher fraction of the population gets infected in a healthier society (MHP). This is a consequence of a relative high number of individuals with mild symptoms, which are more likely to further engage in social interactions. In contrast, in a society with more abundant chronic health problems (LHP), individuals are more likely to develop severe symptoms and are less likely to transmit the disease further. Consequently, in the LHP society, there are more hospitalizations and deaths. Interestingly, it appears that how individuals with poor health are distributed within the population affects the epidemic outcome as well. In a segregated society, where these individuals are located all in one district, mortality is higher when compared to the homogeneously mixed population. This holds true for both types of societies (LHP and MHP) and is more pronounced if the containment measures are mild. Otherwise, the impact of restrictive measures is quite trivial and affects only the absolute values. In continuation, we explore these issues in further detail by examining the disease's progression among individual sub-groups.

In Fig. 4 we show the total number of infections, hospitalizations, and deaths during the first 280 MCS of simulation, separately for each risk group. The results are shown for both types of populations, for both types of distributions of individuals with poor health within the network, and for two different levels of containment policies. We can observe that in the MHS, there are more hospitalizations and deceased in RG1 when compared to the LHS (Fig. 4B and C), which is not only a consequence of a higher number of individuals in this group but also due to higher numbers of infections (Fig. 3G and Fig. 4A). However, it should be noted that the absolute numbers of casualties are in this risk group (RG1) low. On the contrary, the numbers of infections, hospitalizations, and deaths in the subpopulations with chronic comorbidities (RG2 and RG3) were significantly higher in simulations of epidemic spread for the LHS society (Fig. 4D-I). Notably, we can infer some further insights into the role of the population structure. If the population is segregated so that individuals with lower health status are all in one district, there are fewer infections in RG1, whereas, in RG2 and RG3, we observe the opposite. The differences are the most pronounced in RG3, where the number of casualties is approximately 20–30% higher when compared to the homogeneous society. This holds true for both types of society, LHS and MHS. Finally, the effect of containment measures strength influences all the results in a rather expected manner.

Fig. 4

Epidemic outcomes after 280 days within individual risk groups. Shown are the absolute total numbers of infections, hospitalizations and deaths for RG1 (A-C), RG2 (D-F), and RG3 (G-I). Results are presented for two populations (MHP, blue; LHP, orange), for two restriction strengths (Q25 and Q50) and for homogeneous and segregated risk group distributions. Error bars denote standard deviations and bar height, the mean value. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Simulating the effect of vaccination on the evolution of COVID-19 epidemics

In what follows, we examine the epidemic trajectories after the population is progressively vaccinated. We consider two vaccination strategies, i) risky strategy, where individuals in RG3 (first) and RG2 (second) have absolute priority for vaccination, and ii) random strategy, where individuals are vaccinated at random, irrespective of their status. We also consider two vaccination rates, i) 300 individuals/day and ii) 50 individuals/day, which annually amount to 109.5% (full coverage) and 18.25% (partial coverage) of the population, respectively. The first value roughly mimics the rates that are realized in developed countries, whereas the second represents the rates in less developed countries with a poor supply of vaccines. We introduce vaccination on the 281st iteration of our simulation, and in each MCS individuals are vaccinated following the selected strategy. Specifically, in each iteration, 300 (or 50) susceptible individuals are selected either entirely at random or at random first from RG3, then from RG2, and finally from RG1, when there are no more susceptible individuals in RG3 and RG2. By this means, we explore how the distribution of the risk groups in the social network combined with vaccination strategy and the population's overall health affects the mitigation of the disease.

In Fig. 5 we show the progress of the epidemic after the start of vaccination (281st iteration step) with a vaccination rate of 300 individuals/day in the networks with a homogeneous (Fig. 5A-C) and a segregated (Fig. 5D-F) distribution of heterogeneous individuals, for both population types (MHP, blue; LHP, orange), distributions (homogeneous and segregated), and both vaccination strategies (“random” and “risky”). In all simulations, we use the lower value of the restriction strength parameter (Q25), as usually, the non-pharmacological interventions in most countries get more permissive once the vaccination program is initiated. In the social network with homogeneously distributed individuals from different risk groups, both vaccination strategies lead to a mitigation of the disease after around 6 months, whereby the number of active infections drops more quickly if the population is vaccinated at random. Conversely, the number of hospitalizations and casualties is lower if the RG3 and RG2 are vaccinated first. This holds true for both types of populations, MHP and LHP. Interestingly, in simulations with the segregated population model, the disease persisted for longer when the risky individuals were vaccinated first because it took several months before vaccination started in the district with healthier individuals. However, the overall death toll is still lower, as individuals from RG1 mostly contribute to the high numbers of infections, which are very likely to recover from the disease. These results are summarized in the lowermost row of Fig. 5, where we show the change in the total number of additional infected (Fig. 5G), hospitalized (Fig. 5H), and deceased (Fig. 5I) individuals at the end of simulations compared to the starting point of vaccination. Indeed, if the priority is to vaccinate risky individuals first, the number of casualties is lower, irrespective of the type and spatial structure of the population. However, in this case, the disease persists longer, and the cumulative number of infected individuals is higher. As a result of this interplay between higher numbers of infections and protected individuals, the number of hospitalizations remains similar for both types of vaccination strategies.

Fig. 5

Simulating vaccination control for the COVID-19 epidemic with a rate 300 individuals per day. (A-C) Results for the homogeneous distribution of risk groups in the social network: (A) share of actively infected individuals, (B) share of hospitalized individuals, (C) cumulative share of deaths. (D-F) Results for the segregated distribution of risk groups in the social network: (E) share of actively infected individuals, (F) share of hospitalized individuals, (G) cumulative share of deaths. (G-I) The total increase with respect to the onset of vaccination in the shares of (I) infected, (J) hospitalized, and (K) deceased individuals for a given population type (MHP, blue; LHP, orange), distribution (homogeneous and segregated), and vaccination strategy (random and risky). Restriction strength was set to Q25 and vaccination rate to 300 individuals/day. Highlighted area around the curves in panels (A-F) gives the 95% confidence interval of the 50 realized simulations. Error bars in panels (G-I) are the standard deviations of the corresponding feature among different simulation runs. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

We next explored how the epidemic behaves if a lower vaccination rate, i.e., 50 individuals/day, is implemented. The simulations were again performed for both types of populations (MHP and LHP) and for both vaccination strategies (random or priority for risky groups), with a restriction strength set as before (Q25). In Fig. 6 we show the progress of the epidemic after the start of vaccination (281st iteration step) for either homogeneous (Fig. 6A-C) or segregated (Fig. 6D-F) distribution of risk groups within the social network. It can be seen that if the population is vaccinated at random, the epidemics is progressively suppressed in all scenarios, but much slower when compared to the higher vaccination rate (see Fig. 5). Notably, when priority for vaccination is given to individuals from RG3 and RG2, the epidemic is not mitigated. In a society where individuals are distributed homogeneously across the population, the epidemic persists, but the plateau of daily infected is lower than the situation before vaccination. On the other hand, in a segregated society, the daily share of infected individuals even increases. This increase is predominately driven by the increasing numbers of infections among individuals in the RG1 subgroup, who are not vaccinated at all, and, in addition, they progressively lose immunity after the first infection. Furthermore, for the LHP and the homogeneous distribution of individuals within the network, somewhat surprisingly, the cumulative number of deceased individuals keeps rising when priority is given to individuals from risky groups, and the outcome is worse when compared to the random vaccination strategy. Namely, it is not possible to protect all individuals with a low health status at this vaccination rate, and since they interact with individuals from RG1, among which the level of infections is high, they often get infected. Remarkably, in the segregated population structure, the opposite is noticed, as the more vulnerable individuals are, in this case, surrounded mostly by vaccinated individuals within their districts.

Fig. 6

Simulating vaccination control for the COVID-19 epidemic with a rate 50 per day. (A-C) Results for the homogeneous distribution of risk groups in the social network: (A) share of actively infected individuals, (B) share of hospitalized individuals, (C) cumulative share of deaths. (D-F) Results for the segregated distribution of risk groups in the social network: (E) share of actively infected individuals, (F) share of hospitalized individuals, (G) cumulative share of deaths. (G-I) The total increase with respect to the onset of vaccination in the shares of (I) infected, (J) hospitalized, and (K) deceased individuals for a given population type (MHP, blue; LHP, orange), distribution (homogeneous and segregated), and vaccination strategy (random and risky). Restriction strength was set to Q25 and vaccination rate to 50/day. Highlighted area around the curves in panels (A-F) gives the 95% confidence interval of the 50 realized simulations. Error bars in panels (G-I) are the standard deviations of the corresponding feature among different simulation runs. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

To provide a more detailed overview of the results, we plot in the lowermost row of Fig. 6 the change in the total share of infected (Fig. 6G), hospitalized (Fig. 6H), and deceased (Fig. 6I) individuals at the end of simulations (day 700) compared to the starting point of vaccination (day 281). Evidently, if the supply of vaccines is very limited, prioritizing individuals with lower health status leads to higher numbers of infections, as the epidemic cannot be suppressed. Accordingly, the number of hospitalizations is higher as well. However, the impact on mortality rates in this scenario is not trivial. In the homogeneous MHP, the outcome by either of the vaccination strategies is approximately the same, whereas in the segregated MHP, it is more favorable to vaccinate risky individuals first. This holds true also for the LHP, when the structure of the population is segregated. But if the individuals are distributed homogeneously, giving priority to individuals from RG3 and RG2, turned out to lead to more epidemic deaths.

Discussion

We have developed an extended stochastic SEIR model to simulate the epidemic spreading of COVID-19. The model incorporates a realistic scheme of social interactions and different compartments encompassing presymptomatic, asymptomatic, symptomatic, and hospitalized individuals as well as deaths, reinfections, and vaccination. Emphasis was given to population heterogeneity in terms of socio-demographic and health-related variations between individuals, which manifests themselves in different epidemiological features and therefore influence the course of the disease. Specifically, we considered that individuals with chronic comorbidities are more likely to develop severe symptoms, get hospitalized and die more often, and remain immune for shorter periods. To explore how these heterogeneities affected the course of the disease and how it is affected when vaccination is introduced, we studied the epidemic dynamics in two types of populations. In the first population type, individuals with metabolic disorders and other related diseases constitute only a small fraction of the whole population, resembling a young demographic structure. In the second population type, the fraction of individuals with chronic comorbidities was much higher, reflecting an older society or the presence of disadvantaged socio-economic groups, which are at greater risk of poor health. Two separate scenarios were considered how individuals with weaker health status were distributed among the society. In the first case, individuals belonging to different health categories were distributed homogeneously among the population, whereas in the second case, all chronically ill individuals were located in one district. The latter resembles the extreme socio-demographic pattern of all individuals with chronic health conditions living in one deprived area. On top of that, we incorporated in our model the effect of non-pharmaceutical interventions by reducing the contact rates based on the fraction of infected individuals and hospitalized individuals. Two types of regulations and containment measures were considered, one more and the other less strict.

Our simulations have revealed that the demographic structure as well as the spatial arrangement of the population affects the epidemic outcome. It turned out that in the first 280 days of the epidemics without vaccination, there were more infections in the healthier society, even though we have not explicitly considered in our model that healthy and younger individuals are expected to engage in social interactions more often. This would certainly lead to even higher differences in infections, as the social interactions were identified as one of the key parameters affecting the spread of COVID-19 [27], [67]. In contrast, the fractions of hospitalized and deceased individuals were significantly higher in the population with a higher share of individuals with chronic comorbidities and other related diseases, which is a rather expected result. Interestingly, the segregation of chronically ill individuals has led to fewer infections but also to more hospitalizations and deaths when compared to the homogeneous population. Noteworthy, the composition of the society was recognized as even more important after the onset of vaccination, particularly with regards to the vaccination strategy. If priority was given to individuals with lower health status the epidemic persisted and the number of infections in the population was much higher when compared to the random vaccination strategy. If the supply of vaccines was rather high, i.e., 300/day/105, vaccinating less healthy individuals first turned out more favorable, as the number of epidemic deaths was lower, irrespective of the type and structure of the population. If the supply of vaccines is very limited, i.e. 50/day/105, prioritizing individuals with lower health status leads to much higher numbers of infections, as the epidemics cannot be suppressed. In this case, vaccinating first chronically ill individuals leads to significantly more hospitalizations on the population level, whereas the number of casualties in the homogeneous less healthy population is even higher when compared to the random vaccination strategy. However, if individuals with chronic comorbidities are isolated in one district, the prioritizing strategy is much more beneficial. These are important aspects to consider in the design of efficient containment strategies. If the population is well-mixed and the supply of vaccine limited, so that the epidemic cannot be promptly suppressed, it might be more beneficial to vaccinate first individuals, who are more likely to spread the disease, even though they are not directly endangered, but will likely transmit the virus to individuals, who are in risky groups. Because the immune responses in less healthy individuals are more short-lived than in more healthy individuals, it might be more effective to concentrate more on the indirect protection, but for COVID-19 these issues are not yet completely resolved [68]. Nevertheless, if individuals with poor health status are rather isolated, our results indicate that priority should be given to them in either case, as they become protected due to the immunity of individuals within their district. This parallels with findings about influenza, where vaccination was found very beneficial if elderly people are living in communities, such as elderly homes [69].

In our model, we have included the essential compartments that are commonly employed in models aimed to simulate the spread of the new coronavirus disease and utilized realistic parameter values, where possible. As such, the model allows, in principle for an elementary description and projections of the evolving COVID-19 epidemics. However, we must still be aware that it contains several inherent limitations. First, high variability in country-specific COVID-19 statistics is observed [65], which can be accounted for by several factors like demography, pool of asymptomatic cases, and the health system structure. Solving this problem of unbiasing data is still an open issue and is particularly important for the development of epidemiological models as well. Second, the spreading of infectious diseases and human behavior are intertwined and is, particularly in the case with the SARS-CoV-2 virus, gradually affected by non-pharmaceutical interventions, which are calling for a coordinated international strategy [70], [71], as well as by endogenous reductions of social interactions [14], [72], [73]. Our study included this aspect in a very simplified form, i.e., by affecting contact rates by the number of symptomatic infected individuals and hospitalized individuals. This relation is, in reality, much more complicated and can result in diverse epidemic trajectories [73], [74]. Third, while our model accounts for human population heterogeneity, our representation which considers a three-tiered subpopulations scheme with different systemic metabolic disturbances, is very simplified. Not only that in reality, the immunological and health heterogeneity of the population is much more complex and connected with the socio-demographic structure, but also the extent of social engagement depends on the health status of individuals, which is not considered directly in our model [75]. Fourth, in our simulations we have used a rather small and isolated social network with 105 nodes, which warrants only a qualitative description of the COVID-19 spread but cannot account for the diverse epidemic patterns observed in different regions around the globe [76]. Fifth, we have considered only two fixed vaccination rates, one high and the other low, which is a simplification of the actual changeable delivery patterns observed in most of the countries. Adjusting an optimal control strategy for vaccine administration with respect to the limited supply in the real world situation is a very complex process [37], [38], [41], [42], [77], but this is beyond the scope of the present study. Finally, there are, of course, several other aspects to consider when modeling the behavior of spreading dynamics of SARS-CoV-2, such as more complex vaccination strategies, mutations of the virus, how strictly a population is following preventive measures, etc., but currently, the major open issues are the uncertainties of immune periods [78], [79], [80]. While all predictions are, by nature, associated with uncertainty, the reliability of the results depends to a significant extent on the duration of immunity and the probability as well as the nature of reoccurring infections. Given the scarcity of available data, it is impossible to have conclusive evidence about immunity at this stage. However, with time the data on antibody kinetics and biologic features of COVID-19 is rapidly increasing, which will facilitate the development of more accurate models.

Nevertheless, our results clearly highlight that health-related heterogeneity and the spatial social structure crucially affect the evolution of the epidemic as well as the future epidemic trajectories after the onset of vaccination. The risk groups and subpopulations are present in all societies. In the USA, for example, it has been shown that the health disparities in nutrition and obesity correlate closely with the racial and ethnic disparities related to COVID-19 severity and mortality. The age adjusted hospitalization rates for COVID-19 among Native Americans and Black Americans were approximately five and four and a half times that of White Americans, respectively. Also, Latin Americans have been hospitalized at a higher rate. i.e., approximately four times that of White Americans [81]. The proportion of COVID-19 mortality, at least in the studied regions and cities in the USA, such as Chicago and Michigan, was more than twice as high as the proportion of Black residents in their geographic area [81]. For example, looking at specific boroughs, like the Bronx, was characterized by the highest rate of hospitalizations and death related to COVID-19 among all five New York City boroughs [82]. Strong evidence exists for the link between metabolic health, in particular severe obesity, and the worse in-hospital outcomes, and higher in-hospital mortality, in this borough [83]. The Bronx has indeed considerably higher rates of obesity and chronic diseases due to the disproportionate amount of poverty and food insecurity, as compared with the other boroughs. Obviously, these disparities make the borough’s predominantly Black and Latin residents more vulnerable to the devastating effects of COVID-19 [81]. These correlations between the regionally specific populations with worse metabolic health, particularly obesity, and the COVID-19 severity and mortality were also extensively studied in other countries and regions. In a very recent study, using data from 30 industrialized countries, Gardiner et al. [84] show that obesity is the factor most strongly associated with the COVID-19 death rate. Obesity, as a chronic, low-grade systemic inflammation, is a common metabolic disorder that concerns not only COVID-19 severity and mortality, but also the outcome of the vaccination. A major concern is that the COVID-19 vaccines will be less effective for individuals with obesity [31], [85]. At the moment, we still lack sufficient data to make any conclusions concerning the effectiveness of COVID-19 vaccination in these risk groups. However, less efficient COVID-19 vaccination in people with obesity and other related diseases would be unsurprising, as vaccinations for other diseases, e.g., influenza, hepatitis, and others, have shown that the vaccines are less effective in people with obesity [30], [32], [86], [87], [88].

All these data show that for a better understanding of the dynamics of virus spread and a better prediction of the outcomes, it is essential that subpopulations at risk for a higher rate of hospitalization and mortality are included in the studies. Computational models are particularly appropriate for such analyses, enabling the recognition of the critical factors for better understanding the epidemic dynamics and contributing to better decision-making processes in societies facing epidemics. Metabolic health, being particularly emphasized in our study, is not the only disadvantage in the low socio-economic brackets of our society. Influenza-related complications and hospitalization rates [89] have demonstrated that preventative and therapeutic health care, limited sick leave, and household structure might also play a role that needs to be considered.

In conclusion, a combination of the complexity of the transmission processes, limited and uncertain available data, extreme heterogeneity of humans and society, variable large-scale mitigation interventions, and limited supplies of vaccines make it impossible to predict the evolution and consequences of the COVID-19 pandemics. To improve our understanding of the transmission dynamics of this complex disease, a plethora of different kinds of models have recently been developed, each of them aimed to enlighten various aspects on how to optimize the strategies and mitigate the epidemics [11], [26], [34], [57], [90]. Our study proposes a network and agent-based epidemiological model, which offers some advantages compared to standard aggregate S(E)IR-type models, such as the inclusions of complex interaction patterns, locality of social contacts, and spatial heterogeneity of the population. Noteworthy, in our simulations, the latter has proven to be a very important factor affecting the trajectories of COVID-19 epidemics, particularly after the onset of vaccination and when the fraction of individuals with a low health status in the population is relatively high. In future studies, these findings should be further evaluated by more detailed models incorporating realistic spatial demographic data, updated information about the duration of immunity periods, and in the context of various spatiotemporal vaccine distribution strategies.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.