Modelling the interplay of SARS-CoV-2 variants in the United Kingdom.

Many COVID-19 vaccines are proving to be highly effective to prevent severe disease and to diminish infections. Their uneven geographical distribution favors the appearance of new variants of concern, as the highly transmissible Delta variant, affecting particularly non-vaccinated people. It is important to device reliable models to analyze the spread of the different variants. A key factor is to consider the effects of vaccination as well as other measures used to contain the pandemic like social behaviour. The stochastic geographical model presented here, fulfills these requirements. It is based on an extended compartmental model that includes various strains and vaccination strategies, allowing to study the emergence and dynamics of the new COVID-19 variants. The model conveniently separates the parameters related to the disease from the ones related to social behavior and mobility restrictions. We applied the model to the United Kingdom by using available data to fit the recurrence of the currently prevalent variants. Our computer simulations allow to describe the appearance of periodic waves and the features that determine the prevalence of certain variants. They also provide useful predictions to help planning future vaccination boosters. We stress that the model could be applied to any other country of interest.

Introduction

Towards the end of 2021 the SARS-CoV-2 virus pandemic has infected more than 500 million people worldwide, with more than 6 million global deaths[1]. Global vaccination has progressed, though at very different rates throughout the world and in the context of a vaccine shortage. While half of the world population has started vaccination schemes, the distribution is highly heterogeneous. In many countries, mainly those of low income, the percentage of the population with immunity processes in progress is lower than 5%[2]. This feature favours the virus spread, and its genome replication originates mutations that lead to new virus variants[3].

Several hundred thousand new variants have been identified in the SARS-CoV-2 genome[4]. Most of them do not affect the transmissibility of the virus because they do not alter the “spike” protein shape. In spite of this fact, the large amount of infections gave rise to some clinically relevant variants (variants of concern), that are transmitted more easily and rapidly than the original ones[5,6] and that can even be more lethal[7]. The most widely spread strains include the ones first detected in: the United Kingdom (Alpha or B.1.1.7), South Africa (Beta or B.1.351), Brazil (Gamma or P.1), and India (Delta or B.1.617.2)[8]. The emergence of new COVID-19 variants raises doubts about the effectiveness of vaccination campaigns[9]. Furthermore, these strains allowed the appearance of fast reinfections[10] and co-infections[11]. The available information suggests that new variants could draw out the pandemic if measures are not taken[12]. Further research is needed to understand the impact of their spread. For instance, Ramos et al.[13] studied the impact of the introduction of a new variant of the virus on a territory, including vaccination campaigns, for the case or Italy. Another recent work by Okabe and Shudo[14] focused on the spread of a new variant on top of the original strain, using two contact network topologies, observing no relevant differences on the effect of the new variant. There is also a study based on Ontario[15] considering the interplay of three different strains, and vaccination, without taking into account the spatial heterogeneity. Accordingly, we propose the development of a stochastic geographic expansion model, that includes several different strains, to study the emergence and dynamics of these new variants. The model is based on a previously used discrete mathematical map that starts from an extended SEIRS model[16-19]. The model is applied to the case of the UK by fitting the actual properties of each strain with the available data up to the appearance of the Delta variant. Vaccination effects are also taken into account in order to study its effectiveness under scenarios with different strains.

Model

The model used in this work is an extension of the one presented by Barreiro et al.[19] that includes the emergence of new virus variants at different times. The map of the country being studied is divided into a two-dimensional grid in which the dynamics occurs. We use a modified SEIRS-V model to simulate the local dynamics (see Fig. 1A) of infection spreading, taking place in each cell (i, j) of an area of a few , with average population density . Global dynamics (Fig. 1B), comprising the entire country, considers geographical spread to neighboring cells and along terrestrial roads or aerial routes. Governmental non-pharmaceutical interventions, and social behavior are reflected in the mobility parameters used (v).

Figure 1

Diagram showing the geo-stochastic model scheme. (A) Represents the local dynamics of a SEIRS-V model. The compartments taken into account are: Susceptible, Recovered, Vaccinated, Exposed and Infected with different strains Vaccinated people can get infected with a certain variant k with a probability . The exposed, infectious and immune periods for the kth strain are , and , respectively. The variables and stand for the vaccination rate and the immunity period conferred by the vaccine. (B) The global dynamics on a geographical area, divided into a grid of cells, is followed by placing a SEIRS-V model on each one and allowing contagions between them. Three mobility processes are considered: movement to neighbor and far cells, and thermal noise. (Diagram was made by N.L. Barreiro using standard free software).

Description of the local dynamics

The compartmental model includes: susceptible (S), exposed to the different variants present in the location, yet not infectious (, ,...,N), after days people become infectious and can transmit the variant with which they got infected (, ,...,N), they stay in this compartment days and then recover (R). Acquired immunity lasts for days and then, a proportion of people who survived, become susceptible again (S). At a specific time, vaccination designed to protect from the initial variant started, being applied at a rate ; vaccinated people (V) stay immune for days but they may be infected with variant k, with a probability (Fig. 1A).

Description of the global dynamics

Geographical spread (Fig. 1B) is simulated by using mobility parameters: for movement between neighbor cells, and for long distance trips by airplane, car, or train. We consider that flows between large cities are greater than flows between small ones, therefore, is weighted by i.e. by the origin and destination population densities. All mobility parameters , interpreted as the probability of traveling from one cell, in which , to another cell, are modeled by a Metropolis Monte-Carlo algorithm (for details see Supplementary Information). When infection is allowed on cell (m, n), one sets and .

Furthermore, there is the possibility of outbreaks in remote and isolated places, because people sporadically travel to unexpected places, regardless of their population. These random movements may be considered as “kinetic energy” kT of the system. In this case we use a Monte-Carlo procedure comparing a random number with the quantity .

Interaction between variants can be included in our model, but until now, very few cases of co-infection with two COVID-19 variants have been reported, so we decided not to include possible interactions in the present work.

As in any model depicting a complex system, it is important to consider its limitations. In this case, one limitation is that people become part of the vaccinated compartment only after receiving the shots recommended by pharmaceutics and having developed immunity. Thus, we are not considering the partial immunity acquired while vaccination is in progress. However, the model could be adapted to include more details on individual dose administration if needed. Another important limitation is that predictions depend on the available information about the mobility parameters. Those parameters depend on population compliance with government interventions and recommendations. Any change in people’s behaviour may swiftly change mobility parameters and consequently modify the virus spreading rate. Finally, most epidemiological parameters, as the periods of natural immunity , infection and vaccine immunity could change depending on the strain. Due to the lack of reliable and complete data about these quantities and for definiteness we used the same values for all the virus variants. This could also limit the prediction capacity of the model.

Results

The model was applied to the case of the United Kingdom (UK). As stated above, to apply the model to a specific country, data about its population density are required. The information to generate the density map was extracted from the “Global High Resolution Population Denominators Project”[20]. The main routes were obtained from Google Maps. In order to compute the geo-stochastic spreading of the virus, the map was divided into a grid of squares, with an area of 25 km (see Fig. 2).

Figure 2

UK normalized population density and main routes map. The map is divided in a grid of squares of 25 km. The map was generated using custom code[21].

The dynamics of different strains of the virus is very well documented in the UK, consequently the available information suffices to fit the parameters of the model for the different variants. For simplicity, only the most abundant strains over time were considered in the adjustment. These are: the EU1 strain (B. 1.177 originated in Europe), the Alpha (B.1.1.7) and the Delta (B.1.617.2). The rest of the variants, less abundant and of less interest, are grouped together under the name “other strains”. Information about the variants was obtained from CoVariants[22] enabled by data from GISAID[23]. We used the same values of the time delay parameters , and for all the variants, since, despite this simplification, the model produces very good results. The upper panel of Table 1 shows the values of the parameters used in the calculations.

Table 1

Model parameters.

Parameter	Meaning	Value
ϵ_k	Period before being infectious	1 day ^(a) ^(b)
σ_k	Infectious period	14 days ^(b)
ω_k	Immunity period	140 days ^(b)
z_k	Survival parameter	0.99
δ	Vaccination immunity period	180/360 days
η	Threshold to start an infection focus	0.00005 ^(c)

The upper panel includes epidemiological parameters common to all strains. The lower panel shows the parameters that differ among strains.

This value was set to 0 for the Delta variant in order to get a better fit.

Values fitted and used with this model in previous works[17–19].

Corresponds to start of an infection focus with 1 infected persn in an averagely populated area.

Values left as zero for simplicity, considering that when massive immunization began, Alpha was the dominant variant.

The vaccination rate was fitted over time in order to adjust the model to the real immunization data (delayed 14 days to ensure immunity development). Data on vaccination and daily cases were obtained from “Our World in Data”[2,24] and Johns Hopkins University[1]. 66% of the population was fully vaccinated by the end of September. Taking into account the willingness to be vaccinated[2,24] in the United Kingdom, however, it is not expected that a large population group will be vaccinated after that. For simplicity, in the model we consider that 70% of the population will be immunized.

We used the stringency index[25] (a measure of the severity of government policies during the pandemic) as a guide to fit mobility over time. In order to reduce the amount of fitting parameters we employed the same value of for both neighbor and far mobility (the latter was weighted by the origin and destination densities). Mobility was adjusted up to August 15, 2021, leaving the same value from then on. Any change in government policies, increasing or decreasing restrictions, could produce a modification of the results. kT was left on a low level (0.1) since, from the beginning of the pandemic, there have been multiple measures to reduce mobility (internal restrictions and bans or quarantine for travellers entering the UK from high-risk regions)[2,25]. The parameter , which reflects how contagious a strain is, was changed to fit the actual curves for each variant. As we mentioned before, new variants could be more contagious[5,6] and suppress previous ones. In this sense, determines how variants will compete and eventually become the dominant strain over time. The lower panel on Table 1 shows the fitted values for each strain. At the same time, some vaccines have exhibited a lower effectiveness against some variants[26]. To take this into account we used the parameter. Given a kth strain, the higher , the more likely that a vaccinated person will get infected by it. Some parameters such as the period of immunity conferred by the vaccine and the efficiency for some variants still need further study. We emphasize here that the immunity period that the vaccines provide, together with the inherent efficiency of the vaccines are the important factors to be considered. Therefore, we decided to propose three possible different scenarios:When computing the dynamics of the model, we adjusted and to fit the new cases of each variant and the real vaccination rate, respectively. Figure 3A shows the fitting for each variant until September 10. Adjusted values of are shown in Table 1. The inset in the figure shows the fraction of infections with each variant. In order to fit the fast case growth of the delta variant (which cannot be explained only in terms of mobility changes) we had to use a very high value and a lower value of (0) (see Table 1). In order to emphasize the evolution of the Delta variant, in the inset we show a longer period of time. From Fig. 3A it is clear that at the current transmission rate this strain will become dominant and, unless a new, more contagious, variant appears, it will prevail in the near future. Figure 3B shows immunization as a function of time. In this case the vaccination rate was changed on the model to fit the actual data with a delay of 14 days. This delay was added to ensure that fully vaccinated people reach actual immunization. The dashed orange curve represent official data while the red solid line corresponds to the model immunization progress. As mentioned before, vaccination stops when 70% population is reached

Figure 3

Model fitting to strain and vaccination data. (A) Daily cases differentiated by strain. Dashed lines represent actual cases of each strain (scaled). Solid lines and shaded areas represent the average and standard deviations obtained from 100 runs of the model, respectively. The inset shows the fraction of cases for each variant for a longer period of time. The Delta variant is expected to prevail if no new more contagious variants appear. (B) Immunization over time. Yellow dash-dot line signals the day immunization started. Red solid line represents the percentage of fully immunized people in the model in time. The dashed orange line depicts the actual two-dose vaccination data (with 14 days delay, which is the average time needed to acquire full immunity after the second dose).

Short vaccine immunity period ( days) and relatively high vaccine efficiency against delta variant ()

Long vaccine immunity period ( days) and low vaccine efficiency against delta variant ()

Long vaccine immunity period ( days) and relatively high vaccine efficiency against delta variant ()

Figure 4 shows the three scenarios mentioned above. The orange solid line represents a 7-day moving average of the new daily cases which are depicted by the blue bars. The start of the vaccination in the model is indicated with a yellow dash-dot line. The model evolution under the first scenario is represented by the green solid line. With a high efficiency vaccine, a low amount of breakthrough infections (infected while fully vaccinated) is expected in the short term, at the current mobility rate. On the other hand, a short immunity period (only 180 days) leads to a massive COVID-19 outbreak for the beginning of 2022, if no restriction measures are imposed in the future.

Figure 4

Model results for three different scenarios. Red, green and violet solid lines are the model simulation results under scenarios 1, 2 and 3 (see text), respectively. Shaded areas on each color represent the standard deviation from 100 model runs for each scenario. Yellow dash-dot vertical line indicates the day when immunization started. Blue bars and solid orange lines are the actual daily cases and their 7 day rolling average, respectively. In scenario 1, a short vaccination immunity period implies a growth of the daily cases in the near future, if current mobility is sustained. Under scenario 2, people are expected to be immune for a longer time and breakthrough infections will act as antibodies boosters, prolonging the defence against the pandemic. In scenario 3 most people will be immune in the near future, lowering the number of cases, and a new wave would appear when vaccine immunity ends.

The model applied to the second scenario is shown in red. As we mentioned before, the model was conceived with the assumption that if a person presents a breakthrough infection, his/her immune system is boosted and therefore he/she will remain immune additional days after infection. This assumption is sensible, given the abundant evidence of increased antibody levels when a previously infected person is vaccinated[27-29]. This is clearly shown by the model dynamics in Fig. 4. In this scenario, while a quite large number of people have a breakthrough infection (given a relatively low effectiveness of the vaccine against delta variant) their immune system is boosted and they stay immune for a longer period of time. Therefore, we expect a larger amount of cases in the short term, while in the long term, more people will remain immune, consistently reducing the number of cases over time. In this scenario, if the majority of the population loses their immunity, a new wave is expected in July 2022. This suggests the need of new vaccines, effective against new variants, in the years to come. Additionally, breakthrough infections could help reducing the number of cases in the long term while not increasing the amount of deaths.

The third scenario is depicted by the violet solid line in the figure. A smaller amount of cases is expected in the near future since most people are fully immunized. Afterwards, however, we expect a new growth of daily cases by the spring 2022. This growth is due to the loss of most people’s immunity. It is interesting to notice that, under our assumptions, the last peak in the 3rd scenario is predicted to be higher than in the 2nd one. This occurs because a lower efficacy of the vaccine implies a greater number of infections and, thus, re-immunized people. This antibody boost reduces the number of susceptible people by July 2022 lowering the height of the cases peak.

More data are still needed in order to decide which scenario is closer to the real one. However, our predictions raise questions about the proper management of the pandemic. How should vaccines be administered? Are non-pharmaceutical interventions a long-term solution? Is the pandemic here to stay? Our results suggest that new COVID-19 waves will come if highly transmissible, vaccine-resistant variants, are present. This would imply the need of new shots including strain-specific proteins. However, mass vaccination is clearly the key to reduce the appearance of new variants of concern and, consequently, the need of those specific shots. There is currently a huge debate in the world on this matter[30-32]. If vaccinated people get a mild version of the infection, COVID-19 will turn out to be a common, relatively harmless, influenza-like disease and (as shown in our model) there will be less susceptible people over time. Under these assumptions, mobility restrictions should be applied mainly to contain the spread of new variants of concern and boosters should only be given to high risk groups in the near future. Long-term booster vaccination is a different discussion to be considered, as there is no information on the duration of its immunity yet. The results of the present study seem to imply that the world’s population would probably need to be vaccinated periodically in the years to come. (See the Supplementary Materials for additional results.

In a changing scenario, during the publication process of this manuscript, the Omicron variant emerged, having a higher transmission coefficient, and becoming the prevalent variant worldwide. Updated results including this strain can be found in the Supplementary Materials.

Discussion

Given the global impact of the COVID-19 pandemic, reliable models to analyze the virus propagation, including the variants of concern, are crucial to explore effective mitigation strategies. Our geo-stochastic multi-variant SEIRS-V model provides means of accurately describing the dynamics of the pandemic, as shown here. The model clearly separates the biological and social parameters. This property enables it to explain some global features observed in the daily number of worldwide cases. For instance, the surge of periodical waves is not strictly related to the appearance of new variants.These unexpected periods, have been seen in various countries with varying vaccination and social-distancing schemes. Our model can explain this global pattern, as naturally related to the biological parameters it includes, such as the immunity time of the recovered patients and the immunity time conferred by vaccines.

In the present work, we applied the model to the UK, and found that it successfully describes the dynamics of the pandemic. Notice that it was essential to fit variant-dependant epidemiological parameters like the transmission coefficients () for the Alpha and Delta variants of the disease. We identified the epidemiological parameters which determine the dominant variant: the one possessing larger transmission coefficient () or smaller exposed period () will eventually become dominant. Interestingly, as new variants with higher transmission coefficient appear, they quickly become the prevalent strain.

We emphasise that the model is applicable to any country with reliable data, as previously shown for Finland, Iceland, Mexico, Spain and Argentina[17-19]. Here we included the main variants of concern for the UK. The social parameters included in the model, as the different mobility types, allow to distinguish the effects of behavioural and cultural differences. Notice that in our model the biological parameters describe the properties of the disease, thus the same parameters could be used for any country, as in particular for Spain and Argentina is the case (see Supplementary Materials).

Furthermore, this model is useful to predict future scenarios, testing pharmaceutical and non-pharmaceutical interventions, and in particular to optimize the timing of vaccination boosters in order to minimize the appearance of new waves of the disease.

Supplementary Information.