Modeling the impact of SARS-CoV-2 variants and vaccines on the spread of COVID-19

The continuous mutation of SARS-CoV-2 opens the possibility of the appearance of new variants of the virus with important differences in its spreading characteristics, mortality rates, etc. On 14 December 2020, the United Kingdom reported a potentially more contagious coronavirus variant, present in that country, which is referred to as VOC 202012/01. On 18 December 2020, the South African government also announced the emergence of a new variant in a scenario similar to that of the UK, which is referred to as variant 501.V2. Another important milestone regarding this pandemic was the beginning, in December 2020, of vaccination campaigns in several countries. There are several vaccines, with different characteristics, developed by various laboratories and research centers. A natural question arises: what could be the impact of these variants and vaccines on the spread of COVID-19? Many models have been proposed to simulate the spread of COVID-19 but, to the best of our knowledge, none of them incorporates the effects of potential SARS-CoV-2 variants together with the vaccines in the spread of COVID-19. We develop here a θ−ij-SVEIHQRD mathematical model able to simulate the possible impact of this type of variants and of the vaccines, together with the main mechanisms influencing the disease spread. The model may be of interest for policy makers, as a tool to evaluate different possible future scenarios. We apply the model to the particular case of Italy (as an example of study case), showing different outcomes. We observe that the vaccines may reduce the infections, but they might not be enough for avoiding a new wave, with the current expected vaccination rates in that country, if the control measures are relaxed. Furthermore, a more contagious variant could increase significantly the cases, becoming the most common way of infection. We show how, even with the pandemic cases slowing down (with an effective reproduction number less than 1) and the disease seeming to be under control, the effective reproduction number of just the new variant may be greater than 1 and, eventually, the number of infections would increase towards a new disease wave. Therefore, a rigorous follow-up of the evolution of the number of infections with any potentially more dangerous new variant is of paramount importance at any stage of the pandemic.

θ-SIR Type model COVID-19 vaccines SARS-CoV-2 variants VOC 202012/01 501.V2 effective reproduction number a b s t r a c t The continuous mutation of SARS-CoV-2 opens the possibility of the appearance of new variants of the virus with important differences in its spreading characteristics, mortality rates, etc. On 14 December 2020, the United Kingdom reported a potentially more contagious coronavirus variant, present in that country, which is referred to as VOC 202012/01. On 18 December 2020, the South African government also announced the emergence of a new variant in a scenario similar to that of the UK, which is referred to as variant 501.V2.
Another important milestone regarding this pandemic was the beginning, in December 2020, of vaccination campaigns in several countries. There are several vaccines, with different characteristics, developed by various laboratories and research centers.
A natural question arises: what could be the impact of these variants and vaccines on the spread of COVID-19?
Many models have been proposed to simulate the spread of COVID-19 but, to the best of our knowledge, none of them incorporates the effects of potential SARS-CoV-2 variants together with the vaccines in the spread of COVID-19. We develop here a θ − i j-SVEIHQRD mathematical model able to simulate the possible impact of this type of variants and of the vaccines, together with the main mechanisms influencing the disease spread. The model may be of interest for policy makers, as a tool to evaluate different possible future scenarios. We apply the model to the particular case of Italy (as an example of study case), showing different outcomes. We observe that the vaccines may reduce the infections, but they might not be enough for avoiding a new wave, with the current expected vaccination rates in that country, if the control measures are relaxed. Furthermore, a more contagious variant could increase significantly the cases, becoming the most common way of infection. We show how, even with the pandemic cases slowing down (with an effective reproduction number less than 1) and the disease seeming to be under control, the effective reproduction number of just the new variant may be greater than 1 and, eventually, the number 2020, the U.S. Food and Drug Administration (FDA) issued the first Emergency Use Authorization (EUA) for a vaccine for the prevention of coronavirus disease 2019 (COVID-19) caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) in individuals 16 years of age and older. The EUA allowed the Comirnaty COVID-19 vaccine (from Pfizer/BioNTech; hereinafter "Comirnaty vaccine"), to be distributed in the USA (see [17] ). Vaccinations in this country began on 14 December 2020, with a critical care nurse in New York City receiving the first injection of a drug developed by Pfizer-BioNTech. On 18 December 2020, the FDA issued an EUA for the second vaccine allowing the Moderna COVID-19 Vaccine to be distributed in the USA for use in individuals 18 years of age and older (see [16] ).
The European Union (EU) reported on 21 December 2020, through the European Medicines Agency (EMA), some recommendations such as granting a conditional marketing authorization for the Comirnaty vaccine in the EU in people from 16 years of age. On 6 January 2021, the EMA recommended granting a conditional marketing authorization for the Moderna COVID-19 vaccine to prevent Coronavirus disease  in people from 18 years of age. This is the second COVID-19 vaccine that the EMA has recommended for authorization (see [13] ). At the time of writing, the EMA Committee for Medicinal Products for Human Use is performing a rolling review procedure for two more vaccines (Vaxzevria from AstraZeneca, hereinafter "Vaxzervia vaccine", and Janssen from Johnson & Johnson, hereinafter "Janssen vaccine", see [10] ). At the moment, it is not yet clear enough from Phase 3 clinical trials how effectively the upcoming COVID-19 vaccines will reduce disease severity and deaths in the population and if they will also prevent infection and decrease transmission of the disease; furthermore, the duration of protection from vaccination against COVID-19 is unknown. Additionally, the initial supply of COVID-19 vaccines is limited and, hence, it is very important for countries to develop a vaccination strategy. Mathematical modeling is essential to identify priority groups for vaccination against COVID-19 and the development of efficient and effective vaccination strategies (see [11] ).
Updated information on the vaccination process in many countries around the world is available in [43] , including the type of vaccine that is being used. Although most countries are using Comirnaty vaccines, some countries use others. For example, Argentina and Russia are using the Sputnik V Vaccine (developed in Russia), China is using the Sinopharm-CNB and Sinovac vaccines and the United Arab Emirates uses the Sinopharm-CNBG vaccine.
A natural question arises: what could be the impact of these variants and vaccines on the spread of COVID-19? It would be interesting to have a tool able to show possible changes in the behavior of the spread of the pandemic.
Many models have been proposed to simulate the spread of COVID-19 (see, e.g., [18 , 52] and the references therein), but, to the best of our knowledge, none of them incorporates the effects of potential SARS-CoV-2 variants together with the vaccines in the spread of COVID-19. In this work, based on the model developed in [29][30][31] and improved in [52] , we develop a θ -SIR type mathematical model (more precisely, a θ − i j-SVEIHQRD model) able to simulate the possible impact of these types of variants and the vaccines, together with the main mechanisms influencing the disease spread. The new model keeps the compartments and parameters used in [52] . In particular, we use the compartments S (susceptible), E (exposed), I (infectious), I u (infectious but undetected that will survive the disease), I D u (infectious but undetected that will die), Q (people in quarantine), H R (hospitalized that will recover), H D (hospitalized that will die), D (dead by COVID-19 and detected), D u (dead by COVID-19 but undetected), R d (recovered after being previously detected as infectious) and R u (recovered after being previously infectious but undetected), but we split every infected compartment (i.e., E, I, I u , I D u , Q, H R , H D ) in different subcompartments E (i ) , I (i ) , where i ∈ { 1 , · · · , M} and M ∈ N is the number of SARS-CoV-2 variants (including the reference one) that are present in the territory under study. Furthermore, we incorporate new compartments V j , for people that are immune by vaccination with the j-th vaccine, where j ∈ { 1 , · · · , J} and J ∈ N is the number of different vaccines being used in the territory under study. We note that people may need to receive several doses and to wait some days (which may depend on the vaccine) to get immunity after vaccination (we point out that it is still not clear if total immunity can be attained with these vaccines, if the vaccinated people may infect others, etc.). For instance, the Comirnaty vaccine needs a second shot 21 days after the first one and wait (in average) 28 days after the first dose to get the maximum efficacy against the virus (see [8,45] ). The Moderna vaccine is also given as two injections, 28 days apart with an efficacy of 94 . 1% at preventing symptomatic COVID-19, measured starting from 14 days after the second dose (see [15] ). We assume that people are in compartment V j once they get their immunity (not when they get one of the doses). For the sake of simplicity, we assume immunity against all the variants of SARS-CoV-2 for people immune by vaccination or recovered from an infection by any of the variants (other possibilities could be also considered in the model).
The model equations and a diagram summarizing the model (see Fig. 1 ) can be found in Section 2 , together with some other details, like formulae for the estimation of the effective reproduction number of each variant of the virus and that of the global disease. The effective reproduction number of each variant must be less than 1 in order to have the spread of the disease under control. The model developed here may be of interest for policy makers, as a tool to evaluate different possible future scenarios.
In Section 3 , we study the impact in Italy, as a particular example of study case, of the introduction of a new variant of the virus and of the vaccination campaigns. For this, we use a simplified version of the model that is explained in Section 2 .

General mathematical model
Here, we develop a mathematical model to asses the impact of SARS-CoV-2 variants and vaccination campaigns on the spread of COVID-19 in a considered territory. In Fig. 1 , we show a diagram of the general model developed here, under the assumptions made in Section 1 . We remark that, when studying a particular territory, it can be suitable to simplify the general model, according to the data available and the particularities of the case under study. For instance, Fig. 2 below illustrates the diagram of the simplified model used for the particular case of Italy studied in Section 3 .
The evolution of the compartments mentioned above is modeled by the following system of ordinary differential equations: for i ∈ { 1 , · · · , M} , j ∈ { 1 , · · · , J} , where: • M ∈ N is the number of SARS-CoV-2 variants and J ∈ N is the number of vaccines that are being used in the territory under study.
• v j (t) is the number of people that get or lose (e.g., because of changes of the vaccine efficacy through time) immunity after getting the j-th vaccine (we note that this is different from the number of people vaccinated), per unit time, at time t. Here, we have used where u j (t) is the number of people that got their first dose of the j-th vaccine per unit time, at time t, e j (s ) is the efficacy of the j-th vaccine, defined as the ratio of the number of people that is immune by vaccination s days after getting the first dose, and t j is the time when the vaccination campaign of the j-th vaccine started in the territory under study. We assume that all people receiving the first dose of a vaccine will receive, if they do not get infected, the next necessary doses recommended by the laboratories. More complex formulations can be also used, if necessary. Furthermore, the derivative of e j (i.e., ˙ e j ) can be taken in the distributional sense (see, e.g., [51] ), which generalizes the classical notion of derivation and makes it possible to differentiate functions whose derivatives do not exist in the classical sense. In particular, we can use functions e j which are piecewise constant, as those used in the case of Italy studied below.

Remark 1.
We point out that, if the vaccination is not only done for people in compartment S, then we could use more complex versions of v j . For instance, if we assume that people chosen to get the first dose of a vaccine have not been vaccinated previously and are in compartments S, where is an estimation of the number of people that can be chosen to be vaccinated at time τ , t R d > 0 is the delay (in days) that a person that enters the recovered compartment must wait before being vaccinated, is an estimation of the number of people in compartment S that can be chosen to be vaccinated at time τ and we take 4 We note that the last parentheses in the two previous formulae account for the number of vaccinated people (in the whole population or just in the susceptible compartment, respectively) that have not gotten immunity by vaccination, at time τ . Of course, we can use other possible formulae and also include stratification by age, priority people (like people working in the health system or vulnerable people), event-triggered policies (see, e.g., [3] ), etc. As explained above, in this work, for the sake of simplicity, we just assume that v j is given by (3) , since more complex formulae would need the availability of more complete data; the goal here is to get a first idea of possible changes in the behavior of the spread of the pandemic due to new variants of the virus and vaccination campaigns. Note: In the Annex we use (4) (instead of (3) ) with updated data obtained at the time of writing a revision of the initial submission of this paper, before publication.
• A detailed description of the other model parameters can be found in [52] , taking into account that now the superscript (i ) refers to each variant of the virus. For the sake of readability of this work, here we briefly recall their meaning: • N is the number of people in the territory before the start of the pandemic.
• ω (i ) (t) ∈ [0 , 1] is the instantaneous infection detected fatality ratio in the territory, at time t; i.e., the proportion of new detected people (after being in compartment I) that will die of COVID-19, per unit time, compared to the total number of new infectious people (detected or undetected), per unit time.
• ω (i ) u (t) ∈ [0 , 1] is the instantaneous infection undetected fatality ratio in the territory, at time t; i.e., the proportion of new undetected people (after being in compartment I) that will die of COVID-19, per unit time, compared to the total number of new infectious people (detected or undetected), per unit time. 1] is the proportion of the number of new infectious people that are detected and documented by the authorities (after being in compartment I), per unit time, compared to the total number of new infectious people (detected or undetected), per unit time, at time t. If we assume that ω (i ) = (1 + k (i ) ω ) ω (for a certain function ω) and, for the sake of simplicity, the proportion of detected cases is the same for all the variants, i.e., θ (i ) = θ (of course, other possibilities can be also considered) for any i ∈ { 1 , · · · , M} , we may compute θ (t) as 4 We remark that where p (i ) v (t) is an estimation of the ratio of the infections of variant i over the total number of infections (hence, where c r (t) is the number of detected cumulative cases at time t, d r is the estimated number of detected cumulative deaths at time t, s H D (t) is the mean of the duration (days) in compartment H D of a person that enters that compartment, at time t 5 , and r is the minimum natural number such that c r (t) − c r (t − r) = 0 and the previous quotient is not greater than 1 (other details can be seen in [52] ). Here, we assume that s H D is a non-decreasing time dependent function (other cases could need a more elaborated estimation).
are the disease contact rates (day −1 ) of a person in the respective compartments.
We note that some of them may change with time, due to the social distancing control measures. More precisely, we consider that Here β (i ) E, 0 , β (i ) I, 0 , β (i ) I u , 0 and β (i ) I D u , 0 are the respective disease contact rates when no control measures are taken. Additionally, m ∈ [0 , 1] represents the reduction of those specific disease contact rates, at time t, due to the efficiency of the control measures corresponding to social distancing. It is given by (see [29] ) where, for every i ∈ { 0 , ..., q } , constants m i ∈ [0 , 1] measure the intensity of the control measures (lower value implies higher intensity), κ i ∈ [0 , + ∞ ) (day −1 ) simulate the efficiency of the control measures (higher value implies faster effectiveness of the control measures) and λ i ∈ [ t 0 , ∞ ] denote the days of changes of each control strategy. • γ (i ) I (t) ∈ (0 , + ∞ ) is the transition rate (day −1 ) from compartment I (i ) to compartments I (i ) and D (i ) , respectively, at time t. • p (i ) is the ratio of the number of new detected infected people that will survive the disease and are hospitalized (they go to compartment H R ) over the total number of new detected infected people that will survive the disease (and may also do quarantine, going to compartment Q (i ) without going to hospital).
• τ (i ) 1 (t) and τ (i ) 2 (t) are the people infected that arrive/leave from/to other territories, per day.
We point out that, due to the emergency of the situation and the lack of data, we have focused on trying to draw general conclusions without aiming to get very precise forecasts. Therefore, for the sake of simplicity, we have not stratified the model by age. This might result in an overestimation of the deaths, because the priority of the vaccines is being given to old people, who have a higher risk of dying from COVID- 19. We note that in system (1) -(2) we have assumed, as stated above, that people in compartment Q are not infectious. In [52] , it is explained why this assumption is reasonable for many territories (as, for example, Italy, which is the case studied here) and how to model the general case when this assumption is not suitable.
We note that system (1) is uncoupled with system (2) . Therefore, we can solve numerically the whole system (1) -(2) or we can first solve system (1) and then solve each of the equations of (2) as follows: It is important to remark that, if for some territory the model does not include some of the compartments, then these formulae must be changed accordingly.
For the numerical simulations presented in this paper, system (1) was numerically solved, together with the equations of the first line of (2) , by applying the classic 4 th order Runge Kutta method with four stages (RK4) and with four hours as time step (which was tested in [29,52] and produced suitable stable results). The reason of including the equations of the first line of (2) into the RK4 algorithm was to have the code ready for cases including the possibility of infectiousness of people in compartment Q (as explained above, for the sake of simplicity, we have not considered this case here), which would imply a disease contact rate β Q in the equations of the first two lines of (1) , coupling in that system the mentioned equations of the first line of (2) .
Simultaneously, at each of those time steps, we obtain the solution of system (2) (except for its first equation) given by (6) , using the trapezoidal rule for solving the integrals in the formulae. We also note that, in (6) , t 0 can be changed by any other time t < t (for instance, the numerical time step previous to t).
The model is able to generate several outputs of interest, like those described in [52] (taking into account that, here, we introduce the superscript (i ) , which refers to each variant of the virus). One of those outputs is the effective reproduction number of the i -th variant of SARS-CoV-2. It is defined as the number of cases that one infected person of the i -th variant of SARS-CoV-2 generates, on average, over the course of its infectious period. Let us see how to compute it.
The model given by equations (1) -(2) assumes that the duration s X (days) of a person in an infectious compartment (the case of nonconstant transition rates is discussed in Remark 3 , to avoid here more complex notation), is an exponentially distributed random variable (see [2] ), with probability density function given by f (i ) , then the number of cases generated by an infected person that enters compartment X (i ) at time z and leaves it at time z + s X , while being in that compartment, is and, therefore, the expected number of cases generated by an infected person that enters compartment X (i ) at time z, while being in that compartment, is Hence, adding up the contribution of all the infectious compartments related to variant i , the model allows to obtain the following estimation of the effective reproduction number of the i -th variant of SARS-CoV-2: This expression is complex to be applied in numerical simulations, because of the improper integrals. Then, using that is the mean of the duration s X (days) of a person in an infectious compartment X (i ) , we may approximate R (i ) which can be rewritten as where t (i ) Remark 2. Assuming that the fraction of susceptible individuals is constant throughout the calculation of the previous integrals, with and the disease contact rates are constant (which is not realistic for the COVID-19 pandemic, due, for instance, to the changes in time of the control measures), then the effective reproduction number of the i -th variant of SARS-CoV-2 can be also estimated as Of course, many other alternatives, for this or other models can be used to estimate the effective reproduction number (see, for instance, [53] ).

Remark 3. When the transition rate of an infectious compartment
D } is given by a time dependent function γ : (0 , ∞ ) → R such that γ (t) ≥ γ > 0 , for any value t ∈ (0 , ∞ ) (for the sake of simplicity of notation, in this remark we avoid in some occasions the use of subscripts and superscripts), the notation and formulae are a little more complicated. If someone enters compartment X at time z, the duration s (days) of that person in compartment X is an exponentially distributed random variable, with a probability density function given by . Therefore, the expected number of cases generated by an infected person that enters compartment X at time z, while being in that compartment, is and the model allows to obtain the following estimation of the effective reproduction number of the i -th variant of SARS-CoV-2: Now, using that s (i ) is the mean of the duration s X (days) of a person in an infectious compartment X (i ) , if entering at time z in that compartment, we may approximate R (i ) (7) , where now (which generalizes (8) ).
The effective reproduction number of the of SARS-CoV-2 (including all its variants), defined as the expected number of secondary cases produced by one person infected of SARS-CoV-2 (without focusing in a particular variant) at time t, can be estimated by is an estimation of the number of susceptible individuals entering compartment E (i ) , per unit time, at time t.
Typically, the spread of infection by each variant of SARS-CoV-2 slows down when R (i ) e (t) < 1 . Hence, in order to have the spread of the disease under control, we need that max i ∈{ 1 , ··· ,M} R (i ) e (t) < 1 . We note that there may be cases with R e (t) < 1 (and, therefore, the number of cases is expected to decrease in the following days) but max i ∈{ 1 , ··· ,M} R (i ) e (t) > 1 (and, therefore, there is, at least, one variant of the virus that is spreading, implying that the disease is not under control). Once the model is built, if we want to use it for a particular country or territory, we need to identify suitable model parameters and functions for the case under study. Some of them can be found in the literature [33] . However, due to the New cases 1 9 J a n 2 0 2 0 lack of information currently available regarding the behavior of the SARS-CoV-2 and the implicit uncertainty, others need to be estimated. In [52] , one can find the details on how to estimate some of the model parameters and particular values of them suitable for the case of Italy.

Application to the case of Italy
In order to show the possible impact of the introduction of a new variant of the virus in a territory and the vaccination campaigns, we present here, as an example of application, the results obtained with our model for the case of Italy, one of the countries that has been strongly affected by the pandemic. As explained above, Italy is one of the countries where VOC 202012/01 arrived with one case reported on 20 December 2020 (see [35] ), which led the Italian Ministry of Health to forbid flights from the UK to Italy (see [36] ).
In Section 3.1 , we present a simplified version of our model suitable for the case of Italy. Then, in Section 3.2 , we show and discuss the results obtained when we use it to simulate the spread of COVID-19 in this particular country.

Simplified mathematical model for the case of Italy
The case of the mathematical modeling and simulation of the spread of COVID-19 in Italy was studied in [52] , without considering the vaccination campaign which started on 27 December 2020 and taking into account just the reference version (or variant) of SARS-CoV-2.
Here, we study the case of Italy considering vaccination campaigns with the Comirnaty and Moderna vaccines and the introduction in the country of the VOC 202012/01 variant of SARS-CoV-2, which we assume spreads faster than the reference one. For the sake of simplicity, we focus on the case that only the disease contact rates (i.e., the βs) change in the new variant, with respect to the reference one, keeping the same values used in [52] for all the other parameters (which do not need now the superscript notation (i ) ). The only exceptions are the following: a) we will use the control measures parameters m 12 = 0 . 18 and m 13 = 0 . 17 (which correspond to the intensity of the social distancing strategies started on 23 December 2020 and 7 January 2021, respectively), which is in the range of the multiple values used in [52] in order to study different future scenarios; b) we will use the fatality ratio ω = 1 . 4555% , which is one of the three options used in [52] . As it is also done in that work, we do not use in this case a compartment for I D u . Of course, many other possibilities can be studied, considering the model above. With these assumptions, and denoting Q = Q (1) + Q (2) , the simplified model New cases 1 9 J a n 2 0 2 0  is summarized in Fig. 2 and can now be rewritten, for this case, as the following system of equations: for i ∈ { 1 , 2 } , together with for j ∈ { 1 , 2 } .
We will assign the index i = 1 to the reference SARS-CoV-2 virus, the index i = 2 to the new variant VOC 202012/01, the index j = 1 to the Comirnaty vaccine and the index j = 2 to the Moderna vaccine. For the reference virus, we use the values of the parameters given in [52] . To model that the new variant spreads faster than the reference one, we assume that β (2) X = (1 + k ) β (1) X , with k > 0 , for all X ∈ { E, I, I u , I D u , H R , H D } . In other words, the disease contact rates of the variant are 100 × k % higher than the disease contact rates of the reference version of the virus. The initial data considered to run the simulations are detailed in [52] for the reference variant. According to [12] , the first reported cases of the new variant in Italy were imported from the UK to different Italian regions. Here, we enumerate them 0 1 J a n 2 0 2 1 1 1 J a n 2 0 2 1 2 1 J a n 2 0 2 1 3 1 J a n 2 0 2 1 by date of confirmation: 20 December 2020, one case in Lazio; 21-22 December 2020, two cases in Apulia; 23-24 December 2020, two cases in Lombardy; 24 December 2020, three cases in Veneto; 26 December 2020, six cases in Naples. Those detected cases were introduced into system (9) using parameter τ (2) 1 . For the sake of simplicity, we assume that they arrived in exposed state E and we introduce them twelve days before its official date of confirmation (taking into account one day of delay in the notification and eleven days that they spend in average in states E and I before its detection). All the time series used in this work, including the one with the imported cases, are available at https://github.com/momat-ucm/T-SIR-T .
We use equation (3) for both vaccines with the following efficacy functions: • For the Comirnaty vaccine, we take which corresponds to the Comirnaty vaccine efficacy. These values, given in [49] , must be taken with caution, since, as said in that work, their "study was not designed to assess the efficacy of a single-dose regimen". Then, the e 1 derivative in the distributional sense is The Comirnaty vaccine was the first vaccine used in this country, with a first batch of 9,750 doses received on 26 December 2020 (see [38] ). The vaccination campaign continued the following days with a rate of around 470,0 0 0 doses per week (see [39] ). The data of the vaccines already administered (at the time of writing) that we use here for our New cases 1 9 J a n 2 0 2 0     simulations were taken on 17 January 2021 from the Italian official repository [9] . For the sake of simplicity, we do not consider here the possibility of using an additional dose from the vials that were supposed to contain five doses of the Pfizer-BioNtech vaccine. Taking into account this information and assuming that the second dose of the Comirnat vaccine will be administered 21 days after the first one, we use the function u 1 defined as follows: Official reported values , if t < 16 January 2021 470 , 0 0 0 / 7 , if 16 January 2021 ≤ t < 17 January 2021 , In Fig. 3 , we show the total daily amount of remaining and (first and second) administered doses, resulting from the use of this function. We note that, with this choice, from the seventh week, every day there will be a constant number of 470,0 0 0/7 vaccines administered per day (first plus second doses) and function u 1 is periodic, with period of 6 weeks, from the fourth week. Of course, other choices of function u 1 (for instance if it is assumed that fewer doses are administered on weekends, etc.) can be also used. • For the Moderna vaccine, we consider which corresponds to values for the Moderna vaccine efficacy deduced from [1] . Again, these values must be taken with caution, since this study was not designed to assess the efficacy function along different time intervals and few data are included. Cumulative incidence

14-day cumulative incidence
Est. with m i adaptative, i 19 Regarding the Moderna vaccine, according to different media, such as [19,54] , the first 47,0 0 0 doses arrived in Italy between 11 and 12 January 2021, and their administration was programmed to start one week later, on 18 January 2021. Then, the government plans to increase the amount of dispatched doses every two weeks, reaching a total of 764,0 0 0 doses on 22 February 2021. Furthermore, according to [24] , 536,0 0 0 more doses would arrive during March, adding up a total of 1.3 million doses of the Moderna vaccine. Taking into account this information and assuming that the second dose of the Moderna vaccine will be administered 28 days after the first one, we use the function u 2 defined as follows: In the previous function we assume that, ideally, on the first four weeks all the available doses are administered, and in the following ones, priority is given to second shots. We note that this implies to keep some vaccines in deposit in the week from 8 to 21 March 2021, for future second doses in the week from 22 March 2021 to 4 April 2021. In Fig. 4 , we show the total daily amount of remaining and (first and second) administered doses, resulting from the use of this function. We note that, with this choice, from the twelfth week, every day there will be a constant number of 536,0 0 0/24 vaccines administered per day (first plus second doses), and function u 2 is periodic, with period of 8 weeks, from the eighth week. Again, other choices of function u 2 are also possible.

Results for the case of Italy
For the results presented in the section, we used the simplified model, parameters and data described in Section 3.1 . Following the conclusions of the 18 December report of the UK's New and Emerging Respiratory Virus Threats Advisory Group (mentioned above) regarding the transmissibility of VOC-202012/01, we considered in the model an increase of 39% and 93% of the variant's disease contact rates (the βs), with respect to the reference's ones, since the R e value (i.e., the effective reproduction number) is proportional to these parameters. We point out that the other possible increases mentioned above are between those two values. We only report here the results corresponding to the increase of 93%, since the simulations done with an increase of 39% are quite similar to those obtained without considering a second variant, for the time intervals studied here.
We In Section 3.2.3 , we analyze and discuss the obtained results.

The control measures active in Italy as of 17 January 2021 are maintained over time.
In Fig. 5 , we show a possible evolution of the number of daily detected cases and deaths, with and without vaccination campaigns and new variant, when keeping the control measures used at the time of writing.
In Fig. 6 -(Top) , we can see the simulation of the evolution of the daily detected deaths caused by each variant, with and without vaccination campaigns. In Fig. 6 -(Bottom) , we illustrate the possible evolution of the effective reproduction numbers with and without vaccination campaigns and new variant.

The control measures active in Italy as of 17 January 2021 are relaxed on 1 March 2021.
In Figs. 7 and 8 , we show some results of the simulations, if we assume that the control measures are relaxed on 1 March 2021, with m 14 = 0 . 20 . It is a value similar to the one we estimated in [52] for the period between 13 November and 23 December 2020 (before Christmas).
Finally, in Figs. 9 and 10 , we show some results of the simulations, if we also assume that the rate of vaccination is twice the one defined by u 1 and u 2 (see 12 and 14 ), from 18 January 2021.

Discussion
From the figures above, we observe that, for the scenarios without the introduction of the new variant, the pandemic is under control if we keep the control measures as of 17 January 2021 and also if the control measures are slightly relaxed on 1 March 2021. In those cases, the effective reproduction numbers are below 1 and the daily numbers of new cases and deaths are constantly decreasing since January 2021.
On the opposite, when introducing the second variant of the virus, the behaviour of the pandemic depends on the use or not of vaccines. Indeed, without the administration of vaccines, the number of cases and deaths would start increasing in June and the situation would become out of control. Actually, in that context, the problem would have started before that, with the introduction of the second variant, which had an associated effective reproduction number ( R (2) e ) greater than 1. Therefore, even if in January the global effective reproduction number ( R e ) is less than 1 and the number of total cases and deaths is decreasing, the number of cases and deaths associated to the second variant would increase leading the pandemic to a next wave. In other words, in contrast with common thoughts, R e < 1 is not enough for having the spread of the disease under control, if a more contagious new variant is active.
The administration of the vaccine should help to avoid this dramatic situation. In that case, although R (2) e > 1 , its values are lower than in the previous scenario (i.e., without vaccine) and it exhibits a continuous decreasing trend. Thus, the numbers of daily new cases and deaths are considerably reduced in comparison with the previous case. However, although not represented here, eventually there could be a smaller new wave that would finish once R (2) e < 1 . Furthermore, maintaining the control measures to their current level seems to be necessary to keep the situation under control, if we keep the expected (at the time of writing) rate of vaccination. Actually, as can be seen in Figs. 7 and 8 , relaxing those control measures in the near future could produce a huge increase in the pandemic magnitude, with or without vaccines. In order to be able to relax these control measures (which might have a terrible economic impact if they are maintained for a long time), it would be necessary to increase significantly the rate of vaccination. For instance, we show in Figs. 9 and 10 , the simulations obtained when we relax the control measures (as before) but the rate of vaccination is doubled from 18 January 2021. It is clear that, in this scenario, the pandemic would be under control.
From all those results, we conclude that, although the situation in Italy seems to be actually under control, variant VOC 202012/01 (whose associated effective reproduction number could be greater than 1) might produce situations with an uncontrolled new wave even with the current rate of vaccination and control measures. Keeping the control measures and increasing the rate of vaccination seem necessary.

Conclusion
We have developed a mathematical model that is able to simulate the impact of SARS-CoV-2 variants and vaccines on the spread of COVID-19. It is based on the models proposed in [29,52] , which were able to simulate other main mechanisms influencing the disease spread. For the sake of simplicity, we have not included in the model stratification by age, stochasticity and other features that are relevant in this context. The reason is that the goal here is to get a first idea (as soon as possible, because of the emergency of the situation) of possible changes in the behavior of the spread of the pandemic due to new variants of the virus and vaccination campaigns. Furthermore, a more complex model would need the availability of more complete data.
The model is ready (and freely available) to be used with real data in any territory where sufficient information about the disease is gathered by the authorities.
We have applied the model to the particular case of Italy, which has allowed to draw some conclusions. The results show that the current rates of vaccination in that country might not be enough, for example if the control measures are relaxed, to avoid a next wave of the disease. Actually, a potentially more contagious new variant (like the so-called British Variant: VOC 202012/01) may lead the disease to a new wave unless the vaccination rate is increased or very strict control measures are maintained for a long time.
We also obtain formulae for the effective reproduction number ( R (i ) e ) of each variant ( i ) and also for the effective reproduction number ( R e ) of the whole disease. Interestingly, we show that, in contrast with common thoughts, R e < 1 is not enough for having the spread of the disease under control, if a more contagious variant i , with R (i ) e > 1 , is active. Note: The results presented previously were obtained with official data reported up to 17 January 2021 and submitted to the journal on 22 January 2022. In this final version of the article (submitted to the journal on 26 May 2021), following a reviewer's suggestion we have added in the Annex an update of the situation with official data reported up to 25 May 2021. This update confirms what was said in the conclusions: • A new variant has lead the disease to a new wave (the third one) in Italy (see Fig. 14 ). • Although R e < 1 when submitting the first version of this article, the pandemic was not under control, with R (2) e > 1 and eventually R e getting bigger than 1 (see Fig. 15 ).
We point out that the simulations of the scenarios studied in Section 3.2 (which correspond to the initial submission of this article) fit quite well the official data reported by the Italian authorities for one month approximately concerning new cases and even for a longer period of time regarding new deaths (see Fig. 11 ). After that month, there are differences in the dates and magnitudes of the third wave, due to several reasons (see the Annex): poor data available in the literature (when the first version of this article was submitted) about people in Italy infected with the new variant, the Italian authorities relaxed the control measures earlier than the date considered in the scenarios, more information appeared later in the literature showing a higher new variant mortality rate for the new variant, the official data/estimations regarding the vaccination campaign have changed, difficulty of the SIR type models to estimate the evolution of an epidemic/pandemic for long periods of time (due to the sensitivity of its equations to some of their parameters), etc. 1 fourteen days with respect to the scenarios studied in Section 3.2 . Furthermore, there exist suspicions that this new variant arrived to Italy before those dates (see, e.g., [34] or [55] ). Thus, we consider two additional first cases on 8 October 2020, which means τ (2) 1 (8 October 2020 ) = 2 . These values have been obtained by calibration, taking into account the percentages of the new variant officially reported by the Italian authorities in [28] and detailed below. 2. Update of the control measures : In the scenarios studied in Section 3.2 , we considered (just as a hypothetical situation) a relaxation of the control measures on 1 March 2021. Nevertheless, the Italian authorities relaxed the measures on 1 February 2021 (see [41] ) which generated, from a strictly quantitative point of view, differences between the scenarios and the observations. We have also changed the control measures on 23 December and 7 January after a calibration using more data than those available when submitting the first version of this paper. In Table 1 , we include the parameters used to model the control measures from 23 December 2020. The parameters considered before this date are detailed in our previous work [52] . The dates of application of the different control measures in Table 1 have been consulted in the following references [21,41] and [44] , and the efficiency and intensity values have been calibrated using the official data available in [46] .
Furthermore, one of the major benefits of the vaccination campaign is the reduction of the disease death rate (see, e.g., [50] ). For the sake of simplicity, here, we assume that vaccinated people will not die due to COVID-19. Furthermore, regarding the current vaccination strategy in Italy, vulnerable people are vaccinated in priority (currently, one of the main criteria is the age of a person, despite the fact that other factors are also taken into account such as the professions at risk, see [37] ). It has been observed that older people present higher death rates than younger ones (see, e.g., [32] ). Thus, the vaccination campaign implies a decrease in the overall disease death rate. To model this effect, we have considered a suitable decreasing function c ω (t) ∈ [0 , 1] that multiplies the value of ω (1) previous to the Italian vaccination campaign (i.e., ω (1) = 1 . 4555% , as stated in Section 3.1 ). More precisely, we use ω (1) where N age ∈ N is the number of groups used to classify the population by its age (here, intervals of 10 years, [37] ) and ordered from younger (i.e., i = 1 ) to older (i.e., i = N age ) people, N i (t) is the number of non-vaccinated people with an age included in group i at time t, ω CFR ,i is the CFR of people with an age in interval i reported on 16 December 2020 in [25] (which is the latest report available before the beginning of the vaccination campaign) and t vd ∈ N is the delay (in days) between the injection of the first dose of vaccine and the protection from death. To estimate t vd , we have considered the official Italian situation reports (see [25] ), and estimated that between 16 December 2021 and 14 April 2021, the CFR decreased by 28.63%. For the sake of simplicity, we assume that this decrease is mainly due to the effect of the vaccination campaign. Thus, in order to fit this value (i.e., c ω (14 April 2021 ) ≈ 71 . 37% ), we set t vd = 12 days and obtain that c ω (14 April 2021 ) = 71 . 32% . Finally, in order to estimate the evolution of N i (t) up to the end of the considered simulation, if t data is the latest date of available data, for t > t data we use the estimation given by is the average daily number of people vaccinated in group i considering the last 7 available data. In the previous formula, we assume that, for the sake of simplicity, the amount of daily vaccinated people in each group of age is similar to the latest observed data. When an age group is void, the remaining doses are distributed to people in the next younger age group.
In Fig. 12 , we plot the estimated function c ω (t) from 27 December 2020 to 31 August 2021. 4. Update of β H R and β H D : During this pandemic, healthcare workers were quite exposed to the virus. Additionally, their work was fundamental to control and fight the COVID-19. Thus, in many countries (including Italy, see [40] ), they were considered as a priority group to be vaccinated. The immediate effect of vaccinating the healthcare workers is the decrease of the cases due to contact with infected people in hospital. To simulate this in our model, we multiply β H R (t) and β H D (t) by a coefficient c β H (t) ∈ [0 , 1] given by N H , if t > t vac , and 1 elsewhere, where N H = 1 , 974 , 324 is the number of healthcare workers in Italy (according to [40] ) and N v , H (t) is the number of vaccinated healthcare workers. We estimate this function as follows Official data from [9] , if t ≤ t data , H , where t data , H is the latest date of available data and ˜ is the average daily number of healthcare workers vaccinated considering the last 7 available data. In Fig. 12 , we plot the estimated function c β H (t) from 27 December 2020 to 31 August 2021.
5. Update of γ H D : According to the official Italian situation report published on 1 March 2021 (see [26] ), the average time (considering the whole amount of reported cases) between hospitalization and death has increased with respect to the previous report published the 16 December 2020. Thus, to simulate this increase, we consider that γ H D (t) = 1 5 (days −1 ) when t is before t 1 = 18 January 2021 (we remind that we used official data reported up to 17 January 2021 for the initial submission of this paper), γ H D (t) = 1 14 (days −1 ) when t is after t 2 = 1 March 2021, and it increases linearly between those two dates. This particular function provides a good fitting between the reported and simulated deaths. Then, following the notation in Remark 3 , the mean duration in compartment H D of an individual that enters that compartment at time z is where T 1 = max (0 , t 1 − z) < T 2 = max (0 , t 2 − z) . In Fig. 13 , we show the evolution of the mean duration s H D (z) depending on z from t 0 = 1 January 2021. We point out that this is taken into account in the computation of t (i ) H D , which is used to obtain R (i ) e (t) , as explained in Remark 3 , and θ (t) , as explained in Section 2 .
6. Update of p (2) v : Since we assume in this Annex that, as said above, ω (2) = (1 + k ω ) ω (1) = ω (1) , when computing θ using (5) we need an estimation of p (2) v (t) . According to [28] , the cases due to the new variant represent 17.8% of the total detected cases in Italy on 5 February 2021, 54% on 18 February, 86.7% on 18 March and 91.6% on 15 April. Furthermore, we assume that there were no cases due to the new variant before 8 October 2020 and that the new variant will represent 99% of the new detected cases on 31 August 2021. Hence, we obtain p (2) v (t) by interpolating all those data using a piecewise cubic Hermite interpolating polynomial. 7. Update of the vaccination campaign data : We have updated the data used in the model to simulate the Italian vaccination campaign. Here, we use the new data and estimations reported by the authorities, which have changed with respect to the information available when sending the initial version of this article. More precisely: We use the official reported data of the administered doses of the four vaccines used in Italy (Comirnaty, Moderna, Vaxzevria and Janssen) until 24 May 2021. From 25 May 2021 to 31 August 2021, Italy will receive new amounts of doses of these vaccines 6 . Besides, they aim to vaccinate 0.5M people per day (see [47] ). Considering the data in [9] and [48] and the unused doses from previous months, we summarize the amount of available doses of each vaccine in Table 2 . Furthermore, to estimate the number of people that get or lose immunity after getting the j-th vaccine, in this Annex we use the more precise formula (4) instead of (3) , with t R d = 135 days (see [27] ). Therefore, we define functions u i , i ∈ { 1 , 2 , 3 , 4 } as follows: • Comirnaty vaccine: Since the second dose is administered 21 days after the first one, we define: where e 1 is given by (11) and (4) may be rewritten in this case as v * 1 (t) = 0 . 524 u 1 (t − 12) .
• Moderna vaccine: Since the second dose is administered 28 days after the first one, we define: where e 2 is given by (13) and (4) may be rewritten in this case as.
• Vaxzevria vaccine: Since the second dose is administered 84 days after the first one (see [56] ), we define: For this vaccine, according to [7] , we consider the following efficacy: Those threshold values for the cumulative incidence are inspired by the Italian Government indicators described in the Italian official bulletins reported in [25] .