key: cord-0836925-tcd3oybo
authors: Jiménez-Rodríguez, Pablo; Muñoz-Fernández, Gustavo A.; Rodrigo-Chocano, José C.; Seoane-Sepúlveda, Juan B.; Weber, Andreas
title: A population structure-sensitive mathematical model assessing the effects of vaccination during the third surge of COVID-19 in Italy
date: 2021-12-30
journal: J Math Anal Appl
DOI: 10.1016/j.jmaa.2021.125975
sha: 8d3205223cbb2591b883272c2462f004aec50d87
doc_id: 836925
cord_uid: tcd3oybo

We provide a non-autonomous mathematical model to describe some of the most relevant parameters associated to the COVID-19 pandemic, such as daily and cumulative deaths, active cases, and cumulative incidence, among others. We will take into consideration the ways in which people from four different age ranges react to the virus. Using an appropriate transmission function, we estimate the impact of the third surge of COVID-19 in Italy. Also, we assess two different vaccination programmes. In one of them, a single shot is administered to all citizens over 16 years old before second shots are available. In the second model, first and second shots are administered to each citizen within, approximately, 20 days of time-gap.

Since the break out of the COVID-19 pandemic at the beginning of 2020, a wealth of mathematical models have been published in an effort to understand the dynamics of the transmission of SARS-COV-2. By the time this manuscript has been finished, beginning of March 2021, there was already a vast amount of literature regarding models on the expansion of COVID-19. We have selected the references [4, 11, 16, 17, 19, 21, 22, 38, 39, 42, 46, 47] . The reader not completely acquainted with mathematical epidemiology may find of interest the following references devoted completely or in part to that topic [1, 2, 3, 6, 7, 9, 12, 13, 15, 18, 43, 20, 23, 44] .

The achievements reached by the mathematical community in the study of epidemiological models can be a great value in order to support controversial political decisions aimed to control the spreading of the COVID-19 pathogen. The classical SIR model is among the tools that have been proved to describe the evolution of the main epidemiological figures of the pandemic. The reader is referred to [24] , where the SIR model is adapted to describe the evolution of infections, recuperations and deaths in three countries (Italy, Spain and the USA) since the months of March and April 2020 until mid November 2020. The aforementioned model is nothing but a (properly modified) SIR model with vital forces and non constant parameters, and it is given by the equations (1.1) (1) Under total or partial lockdown, β and µ * decay exponentially to baseline values β 0 and µ * 0 according to reported data of infections, deaths and recoveries. This is what happened since mid March until the end of June 2020.

(2) Once β and µ * reach their baseline values β 0 and µ * 0 , the observed values of µ * remain stable around µ * 0 , whereas the fluctuations of β give rise to the successive surges of COVID-19 that occurred in Italy, Spain and the USA since the summer of 2020. When µ * , β, γ, µ and λ are appropriately chosen, the system (1.1) is effective in modelling some relevant epidemiological data (see [24] ). Figure 1 shows two additional examples about Italy that do not appear in [24] . Here and in the rest of this manuscript the calibrarion of the models are done on a daily basis. Despite the fact that (1.1) can be quite precise, it shows some limitations. We underline below three issues that should be taken into consideration in other to improve the response of model (1.1):

(1) The exponential decay of the observed values of β can be detected in several territories under lockdown. The USA has already been mentioned above. An additional example, not studied in [24] , is the case of Germany. A graph of β(t) for Germany constructed out of data from the Robert Koch Institute between March and April 2020 can be seen in Figure  3 (right). The same exponential decay can be observed in Figure 3 (left), where a model for β(t) in the USA is obtained out of data from Johns Hopkins University. In all cases, the exponencial decay in the effective contact ratio β can be justified as the result of lockdown and other strict social distancing measures. At the other end of the scale, mortality is associated to the inherent response of the human body to the virus. Hence, in the absence of effective treatments, an abrupt drop in mortality is more difficult to justify. Figure 2 shows that (in the case of Germany) there was not a substantial decay in mortality during the first weeks of the lockdown of March and April 2020. (2) Time and effort should be devoted to explain why mortality by COVID-19 shows such a markedly different evolution in some territories. The accumulation of evidence favoring the idea that official counts of infections and deaths were way below reality at the onset of the pandemic deserves special attention. Model (1.1) only considers documented cases, ignoring the plausible existence of thousands of infected individuals who are not even aware of the fact that they are infected and, possibly, infectious too. (3) In connection with the above two considerations, the population structure of a country may play a relevant role in the effects of the pandemic in aging populations. It is widely known that most of the deaths occur in individuals older than 65, and this is not being regarded in (1.1) (or in many other models for that matter). Issues number 1 and 3 above will be addressed in depth in the subsequent sections. As for issue number 2, in the practical implementation of our model we will not explicitly take into consideration the impact of undetected cases due to the lack of reliable data. Instead, the seroprevalence studies performed in several countries (in particular in Italy) allows us to draw an educated guess on the number of people with a certain degree of immunity (that is, in the class of the recovered) at the beginning of the period of study.

The paper is arranged as follows: In Section 2 we state a refined version of (1.1) that takes into consideration the age structure of the population and the impact of vaccines. Two different vaccination strategies will be assessed later. Before the calibration of the model defined in Section 2, a qualitative study is developed in Section 3. Here we will determine the asymptotic stability of the disease free equilibrium points in two systems closely related to the models defined in the introduction and Section 2. To achive this goal, the basic reproduction numbers of the disease free equilibrium points will be calculated. In Section 4 we calibrate the model introduced in Section 2 with reported data of Italy in the period 16th November 2020-15th December 2020. A validation based on a reduced set of data is done in Section 5. Sections 6 and 7 are devoted to assess quantitatively the consequences of the ongoing (by the beginning of March 2021) third surge of COVID-19 in Italy, and the impact in infections and deaths of two different vaccination schemes. Section 8 presents the main conclusions.

2. Formulation of a SIR-type model for a population with a given age distribution and a vaccination process

Other age-structured models have been used to study the pandemic of COVID-19 (see for instance [5, 8, 40] ). In the model we are about to formulate the whole population at a given time t will be divided into several disjoint age segments [0, a 1 ), [a 1 , a 2 ),. . . ,[a n−2 , a n−1 ), [a n , a n+1 ), where the ages 0 < a 1 < a 2 < . . . < a n < a n+1 are to be determined according to available data. If we set a 0 = 0 and a n+1 = 125 (no living human being is known older than 125), the k-th age interval is given by [a k−1 , a k ) for all k = 1, . . . , n. Let us denote by N k (t) the number of individuals in the k-th age interval. We define

We will also consider the following epidemiological compartments:

• The susceptible people. These individuals are free of the virus and have never had any contact with it. Hence they do not have immunity and are under risk of being infected. The number of suceptibles in the k-th age segment contains S k (t) people at a given time t. We define S(t) = S 1 (t) + S 2 (t) + · · · + S n (t).

• The documented infected people. These individuals have been proved to be infected with a reliable diagnose and the authorities have included them in the official records. Since they are infectious, special social distancing measures and mobility control policies may be applied to them. These people can develop serious symptoms and they might even die. The number of documented infected individuals in the k-th age segment contains I k (t) members at a given time t. We define

• The recovered people. These individuals have developed a natural immunity to the virus and cannot be infected nor are infectious. Immunity is gained either by overcoming the infection or through the action of vaccines. The number of recovered individuals in the k-th age segment contains R k (t) members at a given time t. We define

• To assess the impact of vaccines we define v k1 (t) and v k2 (t) as, respectively, the first and second doses administered on day t to people in the age group k. The following transition rates will be regarded:

• All the infected recover at a rate γ (recoveries per infected and day).

• Contacts between the susceptible in any age group and the infected have an effective transmission rate given by β. Vital dynamics is governed by the following rates:

• The rate of deaths due to other causes than COVID-19 within the age group k is represented by µ k (deaths per inhabitant and day). • The rate of deaths due to COVID-19 within the age group k is given by µ * k (deaths per infected and day).

• We assume that all the fertile population is concentrated in the first age segment, and that the population birth rate relative to this segment is λ (births per inhabitant in the first age group and day). We will also assume that babies born from all mothers are virus free. It is important to mention that vaccines may not be completely effective. A fraction of the people who have been vaccinated may not become immune. Let θ 1 and θ 2 be, respectively, the effectiveness of the first and second doses of the vaccines. We would like to remark that by the time this paper was finished (beginning of March 2021), all the vaccines scheduled in the Italian vaccination program consisted of at most two doses. The values taken by θ 1 and θ 2 range in the interval [0, 1], where 0 means that vaccines never provide immunity, and 1 that the vaccines provide immunity in all cases. The most typical values attributed to θ 1 range between 0.5 and 0.8, whereas θ 2 is commonly between 0.7 and 0.95. Since there were already several anti Covid-19 vaccines in use by the end of 2020, the values assigned to θ 1 and θ 2 could be a weighted average among all the vaccines being used. We shall also mention that some recovered people may be targeted to be vaccinated, and those would not have an impact on the number of the daily new recovered. We will denote by w k (t) the proportion of people susceptible to receive vaccination who do not belong to the recovered group. Besides, we need to consider that the vaccine may have already been effective on an individual after the first shot. Therefore, we will consider that the probability of an individual to get immunity in the second shot while not being immunized after the first one can be calculated as e 2 = θ 2 (1 − θ 1 ). Obviously, the probabilty of getting immunity after the first shot is e 1 = θ 1 . The probabilities e 1 and e 2 will be considered in the model.

For short or medium term predictions, transitions between age groups by aging do not play a relevant role. However, although births may have a minimal impact in epidemiological variables, they are important in estimating the whole population variation in medium term simulations. This is the main reason why we keep births in the model.

It is widely accepted that the official data of many countries (including Italy) about daily infections and deaths fell short from reality at the onset of the pandemic in March and April 2020. However, the massive availability of test swabs since November 2020 makes the official data more reliable. In this paper the impact of undetected cases will not be assessed. Therefore the model will be calibrated using only the official data of detected cases, recoveries and deaths.

In order to minimize the complexity of the model, the entire population of Italy will be divided into only four age groups, namely, from 0 to 39 years, from 40 to 59 years, from 60 to 79 years and 80 years and older (see Table 1 ). Therefore, the variables of the new model will be S k , I k and R k for k = 1, 2, 3, 4, where k = 1 represents the group of younger age and k = 4 the group of the elderly. Also, S = 4 k=1 S k , I = 4 k=1 I k , R = 4 k=1 R k and N = S + I + R, represents the entire population.

We will apply the same parameter β. The reason behind this simplification is the difficulty involved in the estimation of the effective contact ratio between two specific age groups due to the lack of reliable studies. On the contrary, the value of a general β can be easily estimated using the reported number of daily new infections and other epidemiological information.

Another important decision in the model design is the choice of γ. It is proved in [45] that there is a correlation between age and recuperation rates which is statistically significant, but small. Actually, the study carried out in [45] on more than 5,000 Israeli patients show a difference of about one day in the recuperation rates of patients younger than 29 with respect to patients older than 29. Therefore, we believe that it is justified to consider just one recuperation rate for all age groups.

Having in mind all the considerations commented above, the SIR-type model we will consider is thus defined by the following 12 equations:

with k = 1, 2, 3, 4. Observe that, for convenience, we have put λ 1 = λ and λ 2 = λ 3 = λ 4 = 0.

Once the parameters appearing in equations (2.1) have been properly calibrated, the model's outcome for the variables I and R produce excellent approximations on the number of reported active cases and people with immunity respectively. Additionally, other relevant parameters in any pandemic can also be approximated. If we define

then d(t), D(t), i(t), r(t) and Ci(t) are good approximations of, respectively, the daily number of deaths, the aggregated number of deaths, the daily number of new infections, the daily number of recoveries and the cumulative incidence in 1 day per 100,000 inhabitants on day t. An alternative way to express (2.1) is to consider the equivalent systems of difference equations, namely:

with k = 1, 2, 3, 4. Since the data of infections and deaths are usually broadcasted on a daily or weekly basis, t can be measured either in days or weeks. Throughout this manuscript, time will be given in days unless otherwise specified. This has already been considered implicitly in the definition of all the parameters. Observe that the discrete versions of, respectively, d(t), D(t), i(t), r(t) and Ci(t), are

The numerical integration of the equations (2.1) requires specialized software such as Matlab and Mathematica, using the corresponding numerical integration packages for ODEs. Additionally, programming can be used to run the discretized version of the model in the previous equations. Also, Mathematica has been used in some routine but messy calculations in the qualitative study of the next section.

The parameters λ, µ * k , µ k , γ and β will be determined in Section 4 using reliable sources. First we will stop for a while to study the systems (1.1) and (2.1) qualitatively.

This section is devoted to study the stability of the disease free equilibrium points of the models introduced in the previous sections. The study of model (2.1) is complex even considering a small number of age segments. However, several plausible arguments, which will be discussed in this section, can be used to reduce (2.1) to (1.1) for long term studies. Some effort will be devoted to calculate the basic reproduction number of the models introduced so far. In order to accomplish this calculation with mathematical rigor it will be needed to compute the spectral radius of the so-called next generation matrix following the procedure employed in [14] which, in its turn, is fundamented on [10] .

3.1. Procedure to compute the basic reproduction number. Let us write a given epidemiological model in the form 

is the rate at which individuals are transfered into compartment j (for negative values) or leave compartment j (for positive values). Finally, g = (g 1 , . . . , g m ) T is formed by the transition rates among non-infectious states. We will assume all the time that (3.1) has a unique disease free equilibrium point (X e , Y e ) with X e = 0. If certain conditions are accomplished (see a bit further down), then the basic reproduction number R 0 corresponding to P e = (X e , Y e ) is the spectral radius of the so-called next generation matrix. We define the next generation matrix as F V −1 , where

.

Of course, the matrix V has to be non-singular. Hence

The hypotheses that allow us to apply the previous formula are listed below:

(1) The functions F k , V k and g k are in C 2 ([0, N 0 ] n+m ; R).

(2) A disease free population remains disease free, i.e., F(0,

3.2. Stability of the disease free equilibrium point of the model with no age structure. In this section we will compute the basic reproduction number of (1.1) in order to study the stability of its unique disease free equilibrium point. Firstly it is crucial to mention that (1.1) has been conceived to assess the evolution of COVID-19 in a short period of time. For a long term evaluation the following argument suggests that some adaptations are needed. Indeed, adding the equations in (1.1) we arrive at

Therefore the population in this model is not necessarily constant. As a matter of fact, the entire population would be doomed to disappear exponentially if the net growth rate λ − µ were negative since, in that case, N ≤ (λ − µ)N , from which N ≤ N 0 e (λ−µ)t t→∞ −→ 0, where N 0 is the initial population. It could also happen that N t→∞ −→ ∞ if λ − µ > 0 and µ * is small enough, which is fairly unrealistic. A more interesting model from an epidemiological and demographical viewpoint occurs when the population N remains constant. If the birth rate is low and smaller than the death rate, which is the case of many developed countries, the loss of population due to COVID-19 and a negative net growth might be replaced by immigration. The rate at which the population decreases is given by (λ − µ)N − µ * I. Therefore the immigration rate has to be (µ − λ)N + µ * I to balance the loss of population. Let us assume that by means of strict border control, infected individuals are not allowed to enter the country. Let η ∈ [0, 1] be the ratio of immigration with a sufficient degree of immunity gained either by contact with COVID-19 or by effective vaccination.

to the first and third equations in (1.1), respectively, we end up with

We have changed the order of the variables as pointed out in (3.1) .

Observe that in the right hand side of the equations in (3.3) we have polynomial functions, and therefore they are locally Lipschitz. Hence, by virtue of the existence and uniqueness of the Cauchy-Lipschitz Theorem, system (3.3) has a unique local solution passing through every point

The system (3.3) describes now a scenario where the population is constant, say N 0 . Also, as the reader may easily check (we spare the details), the system (3.3) has a unique disease free

In order to apply formula (3.2) we need to check whether conditions (1) through (6) 

Since the coefficients µ * +µ−η(µ−λ) and µ−η(µ−λ) are strictly positive, it follows that S (t 0 ) > 0. This shows that S(t) ≥ 0 for all t ≥ 0. A similar argument proves that R(t) ≥ 0 for all t ≥ 0. Finally, if I(t 0 ) = 0 then the uniqueness of local solutions of (3.3) shows that (I(t), S(t), R(t)) must be the free of disease equilibrium point, and therefore I(t) = 0 for all t ≥ 0.

Theorem 3.2. The basic reproduction number for the disease free equilibrium point of (3.3) is given by

Proof. We just need to follow the steps introduced in Section 3.1 and check whether condition (1) through (6), shown on page 7, are satisfied. The system (3.3) can be written as in (3.1) by setting

On the other hand, the disease-free equillibrium point is

We spare the details of the easy proofs of conditions (1) through (5) on page 7. As for condition (6), the disease-free system Y = g(0, Y ) is given by

for some K ∈ R. It follows straightforwardly that (S f , R f ) is the unique equilibrium point in Ω Y and that it is globally asymptotically stable. Hence, by (3.4) we have

Applying now formula (3.2) we finally arrive at

As a direct consequence of [10, Theorem 2] we have:

where S f and R f are as in 3.4, and its basic reproduction number R 0 given by (3.5). We have:

3.3. Stability of the disease free equilibrium point of the age structured model. The stability of the disease free equilibrium point of (2.1) can be studied in a similar fashion as we did in Section 3.2. However calculations are now much more complex. The explicit calculation of R 0 is only possible for a small number of age groups. Before doing this, we will provide an informal argument that allows us to reduce the study of a model such as (2.1) to the simpler case represented by model (3.3) . This can be done under certain reasonable assumptions. First of all, we will assume that the vaccination procedure is already finished because all targeted population have already been vaccinated. By the time vaccination is completed, we may assume too that the average mortality rate due to COVID-19 has stabilized and is constant, that is

k is assumed to be constant. A similar assumption can be done for the average mortality rate µ due to other causes of death in the rest of the variables, S k , I k and R k with 1 ≤ k ≤ 4 and for the average birth rate λ. Under these hypotheses, adding conveniently the equations in (2.1) we end up with (1.1) (recall that the vaccination procedure will eventually finish, and therefore we can do without the vaccination terms), which is the starting point of the reasoning developed in Section 3.2. A formal study of the age structured model requires more work and complicated calculations. We will only consider a simplified case where only two age classes are regarded and the vaccination process is over. For long term studies, it is essential to have into consideration the aging of the population. Hence an aging rate δ will be used. With all these considerantions, system (2.1) becomes

where we have assumed that all births occur in the class of young people, that is λ 2 = 0. Adding all the equations in (3.6) yields

According to real data of most developed countries, the right hand side of (3.7) would be negative.

In the rest of this section we will assume that µ 1 − λ 1 > 0. Reasoning as we did in Section 3.2, the loss of population can be compensated with immigration. Let η ∈ [0, 1] be the fraction of immigration immune to COVID-19 and, again, let us suppose that border control minimizes the entrance in the country of infected people. Hence, a fraction 1−η of the immigrants are susceptible. Also, it will be assumed that immigrants belong to the category of young people. Once we rectify model (3.7) with the inclusion of the immigrants in the terms just described, we end up with the model

Observe that in (3.8) we have already arranged the variables as in (3.1). If we put I 1 = I 2 = 0, using elementary linear algebra we obtain the unique disease free equilibrium point of (3.8) , namely

The local existence and uniqueness of the solutions of (3.8) is an easy consequence of the fact that all the funcions on the right hand side of (3.8) are locally Lipschitz. Also, as in Theorem 3.1, it is elementary to verify that (3.8) is positively invariant in

Now we are ready to prove the following:

Theorem 3.4. The basic reproduction number for the disease free equilibrium point of (3.8) is given by

where S 1,f and S 2,f are given by (3.9) and (3.10) respectively.

Proof. Following the steps pointed out in Section 3.1, we identify

where K is given by

Therefore,

which allows us to conclude

The characteristic polynomial of the matrix F V −1 is

from which it follows right away that

It only remains to see that conditions (1) through (6) are satisfied. Properties (1) thorough (5) are easily checked. As for condition (6), we have to examine the solutions of the system

In (3.14), adding the first and third equations on the one hand, and the second and the fourth on the other yields

where A = S 1 + R 1 and B = S 2 + R 2 . Obviously, A + B = S 1 + S 2 + R 1 + R 2 = N 0 . Using the latter and elementary integration in (3.15) it follows:

Now, notice that the first and third equations in (3.14) can be written as

or equivalently as

The explicit solutions of (3.16) can be expressed in terms of exp(M t), b(t) and b 0 as

for some K 1 , K 2 ∈ R 2 . Since the eigenvalues of M , namely −δ and −(δ + µ 1 ), are both negative, it follows that

To finish we focus on the second and fourth lines in (3.14)

and the solutions for S 1 and R 1 which, according to (3.17) , are of the form 

where F and G are again linear combinations of e −(δ+µ1)t , e −(δ+µ2)t and e −δt . Hence

which concludes the proof.

As a direct consequence of [10, Theorem 2] we have:

Corollary 3.5. Consider the system (3.8), its unique disease free equilibrium point Q f , defined in (3.9) through (3.12), and its basic reproduction number R 0 given by (3.13). We have:

(1) If R 0 < 1 then Q f is locally asymptotically stable.

(2) If R 0 > 1 then Q f is unstable.

The parameters λ and µ k (k = 1, 2, 3, 4) are easily obtained with accuracy from the ISTAT databases, and are shown in Table 2 . The constant µ is a weighted average of the µ k 's. It is important to notice that on the onset of the pandemic, that is, during the first wave and in a lesser degree during the second wave too, there was some correlation between the µ k 's and the µ * k 's due to the fact that not all deaths by COVID-19 were detected. However, the excess of mortality since the end of 2020 can be attributed mostly to the detected deaths by COVID-19 as reported by ISTAT (see [29] ).

As for β and µ * we have analysed the data of daily and cumulative infections, deaths and recoveries provided by Protezione Civile Italia [33] , and downloaded from [34] . The period used to calibrate β and µ * consisted of the 30 days between the 16th of November 2020 (t = 1) and the 15th of December 2020. Time will be measured in days. Consider the definitions where, on day t, i e (t) are the new reported infections, I e (t) are the active cases, S e (t) are the susceptible and d e (t) are the new reported deaths. During the calibration period S e (t) was close to the total population of Italy, N , shown in Table 1 . We have not used though the latter approximation. On the contrary, to calculate S e (t) we subtract from N the cumulative cases between the onset of the pandemic and day t.

Observe that β e (t) and µ * e (t), as they have been defined, are natural choices for β(t) and µ * (t) respectively. This justifies that β(t) ≈ β e (t) and µ * (t) ≈ µ * e (t). It was shown in [24] that the observed values of β e (t) during the first wave of COVID-19 in Italy decayed exponentially to a baseline value β 0 , which was estimated according to the data gathered during a 40-day period starting on the 20th of March 2020. Notice that the Italian Government decreed a quite strict lockdown only some 10 days earlier. On the 3rd of November 2020 the Italian Council of Ministers approved a set of containment measures to stop the second surge of COVID-19 in Italy. The soaring values of β e (t) observed along October 2020 had already stabilized by the first days of November, and began to decrease steadily for about three weeks. The evolution of β e (t) (and µ * e (t)) since the onset of the pandemic until November 15th 2020 has already been described and modelled in [24] . We continue below that description starting from the 16th of November 2020. In this manuscript we have modelled the decay of the β e (t) s with an exponential (see Figure 4 ) On the other hand, unlike what happened at the onset of the pandemic, the µ * e (t)'s are relatively stable since the end of the first surge (about June 2020). Their values oscillate around their mean (see Figure 4) . The value chosen for µ * is exactly the average of the µ * e (t)'s over the 30 day calibration period, that is, µ * = 9.1416 × 10 −4 deaths per infected and day.

As for the calculation of the death rates µ * k (k = 1, 2, 3, 4) caused by COVID-19, we have employed the weekly reports elaborated by the Istituto Superiore di Sanità (see [30] ). These reports Table 3 . Infection Fatality Ratios of each of the four age groups considered calculated in a sufficiently long timespan (between June 2020 and mid February 2021).

include weekly information about deaths by COVID-19 in each age group since the beginning of June 2020. The IFR (Infection Fatality Ratio) of each of the considered four age groups can be calculated and the result is shown in Table 3 If IF R k is the Infection Fatality Ratio of the age group k, then the following relationships must hold

for k = 1, 2, 3, 4. Also, µ * should be a weighted average of the µ * k 's. As a matter of fact we can put

Observe that the weights I k I show a very stable behavior in the period prior to vaccination. Moreover, the approximation I k I ≈ N k N is legitimate (see Table 5 and Figure 5) , where N k is the number of people in the age group k and N the total population according to Table 1 . Since the quotients µ * k /µ * 1 can be estimated by (4.1) and µ * has already been fixed, we arrive at the following choice of the µ * k 's: Table 4 . Calibration of γ. We consider the average values of the γ e (t)'s in each period.

The value of γ will be fixed in accordance with the observed daily recuperation ratios γ e (t). On day t, γ e (t) is defined by

where r e (t) are the reported daily recuperations on day t. The average values of the γ e (t)'s are quite stable. We will assign to γ the average values of the γ e (t)'s in each of the periods shown in Table 4 .

Once the parameters of the model have been determined, in order to run a simulation we need to obtain the initial values of the variables. According to the downloadable data from Protezione Civile Italia [33, 34] , on the 15th of November 2020 there were roughly 712,000 active cases in Italy and about 421,000 cumulative recoveries. The latter two estimates are also provided by the real-time statistics site Worldometer [36] . Assuming that the population of the entire country was estimated in some 59,641,000 by the end of 2020 (see Table 1 ), it is easy to infer that the susceptibles on the 15th of November 2020 can be estimated in 59, 641, 000 − 712, 000 − 421, 000 = 58, 508, 000 people.

A further adjustment is needed to include in the category of the recovered the considerable amount of Italians who underwent the coronavirus during the first wave, specially in March and April 2020, without the knowledge of the authorities, and perhaps without even their own knowledge. The Health Ministry of Italy together with the statistics agency ISTAT carried out a seroprevalence study between the 25th of May and the 15th of July 2020 to estimate how many people had developed antibodies to the novel coronavirus, even in the absence of symptoms (see [31] ). This study considered a sample of nearly 65,000 people in 2,000 towns and cities, split by sex, occupation and six age brackets. One of the most shocking conclusions of the study was the fact that the number of people with antibodies to SARS-CoV-2 was estimated in 1,482,000 people. However, by mid July 2020, less than 250,000 cases had been reported since February 2020. In orther words, during the first months of the pandemic, only 1 out of 6 cases was detected. Assuming that this considerable gap in the detection and tracking of new cases is confined to the first months of the pandemic, we would only need to add about 1,500,000 additional people to the category of recovered. As a matter of fact, this is what we have done to establish the initial values of S, I and R: I(0) = 712, 000 people, R(0) = 421, 000 + 1, 500, 000 = 1, 921, 000 people, S(0) = 59, 641, 000 − 712, 000 − 1, 921, 000 = 57, 008, 000 people.

We also need to assign an initial value to the variables S k , I k , R k for k = 1, 2, 3, 4. The available reported data only include the distribution of daily infections within age ranges, neglecting the recoveries and active cases. In order to set a reasonable initial value to all the variables S k , I k , R k Table 5 . Relative weights of infections by age group compared to the relative weights of the age groups with respect to the entire population since the beginning of November 2020.

for k = 1, 2, 3, 4 we have analyzed the data of weekly new infections in each of our four age ranges. The data were obtained, once again, from the weekly reports published by the Istituto Superiore di Sanità [30] . Interestingly, the pattern of weekly infections since the beginning of November 2020 is quite stable, as the reader can see in Figure 5 . Moreover, the relative weight of new infections in each category is similar to the weights of the categories with respect to the whole population (see Table 5 ).

There does not appear to be a particularly strong protective effect towards the most vulnerable age categories. On the contrary, it seems that infections maintain a narrow connection with each category weight. We understand that this justifies the election

With the above initial values of the variables and the calibration performed in Section 4, the accuracy of the model is tested in Figure 6 .

The modeling of the second wave of COVID-19 in Italy was done in [24] satisfactorily using (1.1) by adapting the transmission function β(t). Actually it is only β(t) that needs to be adapted in order to describe the behavior of the pandemic in the event of a new wave. The existence of vaccines is also a remarkable fact that need to be treated appropriately. The fact that those Figure 6 . Assesing how the model's outcomes fit the data in four specific epidemiological parameters: the active cases (up-left), the daily cases (up-right), the cumulative deaths (down-left) and the cumulative incidence in 14 days per 100,000 inhabitants who are at higher risk of death are vaccinated first must have a substantial impact in the average mortality of the whole population that we need to estimate. 6.1. A model for β(t) during the third wave. The authors of [24] have not observed any significant alteration of µ * (t) nor γ (not to mention λ nor µ(t)) when studying the second waves of COVID-19 in Spain or Italy, or the second and third waves in the USA. This is consistent with the fact that mortality and recuperation are, in principle, determined by the reaction of the human body to the virus. On the contrary, β(t) depends heavily on sociological behaviour and is subjected to substantial changes due, for instance, to a higher mobility of the population or the relaxation of social distancing measures. In this manuscript we do not assess quantitatively the different factors that shape the function β(t). We refer the interested reader to [17] . The emergence of new and more contagious strains of the virus might also be directly linked to a substantial growth of β(t) (see [41] ). A good example about how to adapt β(t) to model a coronavirus wave can be found in [24] . In this manuscript we suggest a plausible transmission function β(t) to model the third surge of COVID-19 in Italy. Unfortunately, the daily values of β e (t) have been steadily increasing since mid February 2021 (see Figure 7) , and by the time this paper has been finished (beginning of March 2021), it is widely accepted that Italy is on its way to suffer her third surge of COVID-19. Figure 7 also shows the proposed model for β(t) during the third wave. The explicit form of β(t) depends on parameters that we have defined in previous sections, in particular, we recall that With all that in mind β(t) would be

if t ≤ 220.

Observe that t = 73, t = 114, t = 142, t = 152 and t = 220 correspond, respectively, to the 27th of January, the 9th of March, the 6th of April, the 16th of April and the 23rd of June 2021. The reader can compare β(t) with β e (t) in Figure 7 . The choice of β(t) may appear whimsical at first sight. As a matter of fact, the function β(t) is just an educated guess on how the third surge of COVID-19 will evolve in Italy. Until the 9th of March 2021, β(t) fits the observed values of β e (t) either by logarithmic regression (0 ≤ t ≤ 73), or by linear regression (73 ≤ t ≤ 114). For t ≥ 114 we just try to mimic the behaviour of β(t) during the second wave (see [24] ) both in length and in intensity.

6.2. Vaccination programmes. The third wave overlaps with the vaccination programmes. Here we describe two plausible vaccination programmes and assess their effect on infections and deaths separately. The two vaccination schemes are nowadays in use. One of the vaccination schemes is being employed in most of Europe and consists, basically, of administering two separate shots to every citizen (with, approximately, a three week time-gap between them). This plan will be labelled "A". At the other end of the scale, an alternative scheme is being developed in the United Kingdom, where second shots are postponed indefinitely until the entire population have received the first one. This plan will be labeled "B". In either case, we need to fix the daily availability of vaccines. According to the latest data from Protezione Civile (see [35] ), during the first days of March 2021, about 200,000 doses were being administered every day. We have assumed in our simulations that an average of 300,000 shots per day will be available.

There remain to determine four important issues regarding the vaccination process:

(1) The efficacy of vaccines, θ 1 and θ 2 , for the first and second shots respectively.

(2) The delay of the vaccine effect in providing immunity.

(3) The proportion of vaccine-target citizens who do not belong to the recovered, that is, the value for w k (t), where k = 1, 2, 3, 4 represents the age group. (4) The number of vaccines administered to each group every day, namely, v kj (t), k = 1, 2, 3, 4

where j = 1, 2 will stand for the first or second dose. Le us first make a reasonable choice of the efficacies θ 1 and θ 2 . Table 6 shows the efficacies of the main vaccines in use in Italy according to several studies [26, 27, 25, 37] . Table 6 . Efficacy of the 3 main vaccines in use in the Italian vaccination programme.

On the other hand, not all vaccines are available in the same proportion at any time. The estimated availability during each quarter of 2021 corresponding to the vaccines shown in Table 6 has been published by the Italian Ministry of Health in [32] . Considering this information we have calculated a weighted average for the efficacy of the first and second shots of all vaccines in Table 7 Table 7 . Calculation of a weighted average for the efficacies of all vaccines.

Secondly, we have considered that the average delay for the vaccines to be effective is 7 days, which conveys a time lapse of 7 days in the functions v k,j , in other words, v k,j (t) will be replaced by v k,j (t − 7) in the model's equations.

To determine the proportion of people who are susceptible to receive the vaccine and do not belong to the recovered group, we will divide R k (t), the recovered cases in group k and day t, into two subroups

will stand for the recovered cases in group k and day t who had been previously infected and R v k (t) will stand for those susceptible ones who were vaccinated and succeeded in gaining immunity. With this notation, if S k (t − 1) is the amount of susceptible people on day t − 1 and group k,

. Finally, we have assigned to v kj (t) real data about daily vaccinations from the end of 2020 until the 4th of March 2021. We describe below how the functions v kj (t) have been defined in each of the two vaccination schemes under study: 6.2.1. Vaccination plan A. We try to mimic the official vaccination plan in Italy. The first stage of the plan was already finished when this paper was submitted by the beginning of March 2021. The data corresponding to the first stage were downloaded and included in the model. In a second stage all the most vulnerable people will be vaccinated. This group consists of about 2,000,000 in the age groups 1, 2 and 3, and the whole of group 4. Assuming that there are up to 300,000 available shots every day, we have put until all recipients have received two shots. In the next stage we repeat the same procedure with the citizens in group 3 who have not been vaccinated yet, until all of them have received the corresponding two shots. Afterwards, the same is implemented with groups 2 and 1.

In this plan all available shots considered in plan A are administered as first shots until all the members in each of the age categories have received one shot. Second shots would only be administered once the whole population has received one shot, in which case, the process is repeated in order to administer a second shot to each citizen, starting with the most vulnerable people. Nevertheless, second shots at the considered vaccination rate would not start until mid August 2021.

6.3. Impact of the vaccination plans in µ * . If vaccines were administered randomly, the global death rate by COVID-19, µ * , would remain unaltered. However the vaccination plans are designed to minimize the impact of the virus in the most vulnerable sectors of the population. Those who are more exposed to the virus or under higher risk of death must be vaccinated first. Doing so, the percentage of vulnerable people among the new infections will be reduced. Therefore, µ * will also decrease. As vaccinations progress, it is reasonable to assume that µ * will change as follows:

Notice that I k (t) I(t) is the relative weight of group k among all active cases.

7. Quantitative assessment of the ongoing third wave of coronavirus in Italy and the vaccination process

In this section we analyse the outputs of our model in several situations. We are particularly interested in measuring the impact of vaccines in alleviating the consequences of a new surge of COVID-19 in Italy. For this reason we have considered a scenario where vaccines do not exist, comparing the results with a much more optimistic situation where all citizens older than 16 will eventually be vaccinated. The benefits of vaccines are reflected in Figure 8 , where a prediction of active cases is shown using any of the two vaccination schemes described earlier. If any of the vaccination programmes described above is carried ot at the pace of 300,000 doses per day, the impact of the third wave would be reduced to one half in the number of active cases by the end of May 2021. Interestingly, by that time approximately one forth of the population would have received at least one shot.

The results obtained for the two vaccination schemes are not markedly divergent, although model B seems to be slightly better. In particular, the active cases are about a 14.72% less in model B by the end of May.

The number of deaths due to COVID-19 is also expected to decline thanks to vaccines. Figure 9 shows a significant reduction of mortality due to COVID-19 in both vaccination schemes with respect to the model with no vacctination programme. This reduction is about 44.29% for the Figure 8 . Evolution of active cases with or without vaccination programmes vaccination scheme A, and 47.95% for vaccination scheme B. In absolute terms vaccines could save between 38,000 and 41,000 lives in March, April and May alone, depending on whether scheme A or B are considered respectively. The differences observed between the active cases and cumulative deaths predicted under the vaccination programmes A and B suggest that Scheme B is slightly better at least during the first stages of the vaccination process, but it is not significantly better in any case. The predicted number of cumulative deaths in all the considered models can be seen in Figure 9 . It must be warned that although the beneficial effect of vaccination includes tens of thousands of saved lives as well as a notable reduction in infections, it will not avoid the catastrophic impact of the third wave. Both vaccination programme A and B predict more than 44,000 deaths during the months of March, April and May. In addition to that, very high cumulative incidences are expected along the month of April even if the vaccination programes progress with normality. In particular, the cumulative incidence in 14 days per 100,000 inhabitants could reach values close to 1,000 (see Figure 10 ).

The model's outcomes predict a reduction of the population in Italy of about 129,000 people in the months of March, April and May in the scenario with no vaccines. This reduction would be between 91,000 and 94,000 people depending on the vaccination scheme.

It is important to notice that all the predictions made in this section are based on the particular choice of β(t) done in (6.1) and the vaccination rate of 300,000 shots per day considered above. If β(t) evolves differently or the vaccination capacity is improved, model (2.2) outcomes could be completely different.

Recall that the main objectives of this manuscript are the following:

(1) Assess the capacity of vaccines to hamper the spread of the virus and reduce the death toll. (2) Compare two different vaccination Schemes. As can be seen in Figures 8, 9 and 10, both deaths and infections are greatly reduced with the use of vaccines. It requires only one dose in a 25% of the population to achieve approximately a 50% reduction in active cases and deaths. However, infections are far from being stopped since a cumulative incidence of close to 1,000 cases in 100,000 inhabitants in 14 days can be easily attained. Deaths by COVID-19 can still be extremely high. We estimate that more than 44,000 additional deaths may happen in March, April and May 2021, even if the vaccination plans are accomplished as described in this paper. This significant increase of deaths may lead to a sensitive variation of the whole population. The model's outcomes predict a reduction of the whole population in Italy of 0.22% in the scenario with no vaccines, which is quite significant in just a three month period. This reduction would be about a 0.15% in any of the two sceranios decribed for vaccines. The only way to reduce this exceedingly high death toll is to act promptly to lower the rate of effective contacts β(t) combined with a higher rate of vaccination (in this paper we have used a capacity of 300,000 shots per day).

On the other hand, the comparison of the two vaccination plans evaluated in this manuscript does not provide compelling evidence about which programme is best. It seems that plan B prevents more infections than plan A by a narrow margin, however, both plans converge significantly in predicting deaths. 

Infectious diseases un humans: Dynamics and Control

Population biology of infectious diseases: Part I Nature

Case fatality models for epidemics in growing populations

One year of the covid-19 pandemic in galicia: A global view of age-group statistics during three waves

Age-Structured Modeling of COVID-19 Epidemic in the USA

Mathematical models in epidemiology

Model-informed COVID-19 vaccine priorization strategies by age and serostatus

Lockdown measures and their impact on single-and two-age-structured epidemic model for the COVID-19 outbreak in Mexico

Mathematical tools for understanding infectious disease dynamics

Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission

Early Transmission Dynamics in Wuhan, China, of Novel Coronavirus-Infected Pneumonia

Topics in mathematical biology

Exact analytical solutions of the Susceptible-Infected-Recovered (SIR) epidemic model and of the SIR model with equal death and birth rates

Stability and sensitivity analysis of the epidemiological model BE-CODIS predictiong the spread of human diseases between countries

Age-structured population dynamics in demography and epidemiology

Mathematical modeling of the spread of the coronavirus disease 2019 (COVID-19) taking into account the undetected infections. The case of China, Communications in Nonlinear Science and Numerical Simulation

A simple but complex enough θ-SIR type model to be used with COVID-19 real data. Application to the case of Italy

A Contribution to the Mathematical Theory of Epidemics

A conceptual model for the coronavirus disease 2019 (COVID-19) outbreak in Wuhan, China with individual reaction and governmental action

An introduction to mathematical epidemiology

A conceptual model for the coronavirus disease 2019 (COVID-19) outbreak in Wuhan, China with individual reaction and governmental action

Strong correlations between power-law growth of COVID-19 in four continents and the inefficiency of soft quarantine strategies

Simple mathematical models with very complicated dynamics

Seoane-Sepúlveda. A SIR-type model describing the successive waves of COVID-19

Organization: EMA assesment report on COVID-19 mNRA vaccine

Organization: FDA briefing Document Moderna COVID-19 Vaccine

Organization: FDA briefing Document Janssen Ad26.COV2.S Vaccine for the Prevention of COVID-19

Organization: Ministero della Salute and ISTAT on seroprevalence of SARS

Organization: Ministero della Salute: on anti-SARS-CoV-2/COVID-19 strategic vaccination plan

situaziones Italia

Protezione Civile Italia. Data on daily vaccinations

Safety and efficacy of the ChAdOx1 nCoV-19 vaccine (AZD1222) against SARS-CoV-2: an interim analysis of four randomised controlled trials in Brazil, South Africa, and the UK. The Lancet

Mathematical modeling of COVID-19 transmission dynamics with a case study of Wuhan

Fractional model of COVID-19 applied to Galicia, Spain and Portugal

A modified age-structured SIR model for COVID-19 type viruses

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

Power-series solution of compartmental epidemiological models

Mathematics in population biology

Epidemic and demographic interaction in the spread of potentially fatal diseases in growing populations

Effects of age and sex on recovery from COVID-19: Analysis of 5769 Israeli patients

A new coronavirus associated with human respiratory disease in China

for the China Medical Treatment Expert Group for Covid-19, Clinical Characteristics of Coronavirus Disease 2019 in China