A population structure-sensitive mathematical model assessing the effects of vaccination during the third surge of COVID-19 in Italy.

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.

Introduction

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 [3], [10], [15], [17], [19], [20], [21], [25], [26], [40], [45], [46]. 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], [5], [6], [9], [11], [12], [13], [14], [18], [22], [23], [41], [42].

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 where the system variables, S, I and R stand, respectively, for the susceptible, the infected, and the recovered. The system parameters λ, μ, β, γ, μ and are not assumed to be necessarily constant. Their meaning is described as follows: β is the effective contact ratio, λ is the birth rate, is the death rate due to COVID-19, μ is the death rate due to other causes and γ is the recuperation or recovering rate. Model (1.1) describes the evolution of infections, deaths and other epidemiological parameters, like cumulative incidence, in the three countries under study in [24]. Essentially only two assumptions account for the success of model (1.1) in describing relevant epidemiological parameters: When , β, γ, μ and λ are appropriately chosen, the system (1.1) is effective in modelling some relevant epidemiological data (see [24]). Fig. 1 shows two additional examples about Italy that do not appear in [24]. Here and in the rest of this manuscript the calibration of the models is done on a daily basis.

Fig. 1

Daily number of cases (left) and the cumulative incidence in 14 days per 100,000 inhabitants (right) can be accurately described by the model, when properly calibrated.

Under total or partial lockdown, β and decay exponentially to baseline values and according to reported data of infections, deaths and recoveries. This is what happened since mid March until the end of June 2020.

Once β and reach their baseline values and , the observed values of remain stable around , 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.

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 order to improve the response of model (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 for Germany constructed out of data from the Robert Koch Institute between March and April 2020 can be seen in Fig. 3 (right). The same exponential decay can be observed in Fig. 3 (left), where a model for in the USA is obtained out of data from Johns Hopkins University. In all cases, the exponential 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. Fig. 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.

Fig. 3

In the left it is shown how the β(t)'s decay exponentially to β0 = 5.02742 × 10−11 (day×person)−1 in a seventy day period in the USA. On the right, a similar exponential decay is observed in Germany to β0 = 3.36449 × 10−11 (day×person)−1.

Fig. 2

Daily mortality ratio in Germany between the end of March and mid June 2020.

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 favouring 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.

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 ageing 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) allow 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 on infections and deaths of two different vaccination schemes. Section 8 presents the main conclusions.

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 [4], [8], [38]). In the model we are about to formulate the whole population at a given time t will be divided into several disjoint age segments , , …, , , where the ages are to be determined according to available data. If we set and (no living human being is known older than 125), the k-th age interval is given by for all . Let us denote by 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 susceptible in the k-th age segment contains people at a given time t. We define

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 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 members at a given time t. We define

To assess the impact of vaccines we define and 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 (deaths per inhabitant and day).

The rate of deaths due to COVID-19 within the age group k is given by (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 and 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 and range in the interval , where 0 means that vaccines never provide immunity, and 1 that the vaccines provide immunity in all cases. The most typical values attributed to range between 0.5 and 0.8, whereas 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 and 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 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 . Obviously, the probability of getting immunity after the first shot is . The probabilities and will be considered in the model.

For short or medium term predictions, transitions between age groups by ageing 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 , and for , where represents the group of younger age and the group of the elderly. Also, , , and , represents the entire population.

Table 1

Age distribution of the Italian population by the end of 2020 according to ISTAT (Istituto Nazionale di Statistica).

Age group	Number of people
From 0 to 39	23,537,624
From 40 to 59	18,351,424
From 60 to 79	13,332,737
80 and older	4,419,703
Total	59,641,488

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 [44] that there is a correlation between age and recuperation rates which is statistically significant, but small. Actually, the study carried out in [44] 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 . Observe that, for convenience, we have put and .

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 , , , and 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 . 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, , , , and , 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 λ, , , γ 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.

Qualitative assessment of the disease free equilibrium points

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 [16] which, in its turn, is fundamented on [43].

Procedure to compute the basic reproduction number

Let us write a given epidemiological model in the form where is the vector of all the variables representing infectious states, contains the variables of the non-infected states, the functions and are such that () is the rate of creation of new infections in compartment i and () is the rate at which individuals are transfered into compartment j (for negative values) or leave compartment j (for positive values). Finally, 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 with . If certain conditions are accomplished (see a bit further down), then the basic reproduction number corresponding to is the spectral radius of the so-called next generation matrix. We define the next generation matrix as , where Of course, the matrix has to be non-singular. Hence The hypotheses that allow us to apply the previous formula are listed below:

The functions , and are in .

A disease free population remains disease free, i.e., .

for every if for all and .

if and whenever .

for all X, Y with () and ().

The disease-free system has a unique equilibrium which is globally asymptotically stable in .

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, , from which , where is the initial population. It could also happen that if 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 . Therefore the immigration rate has to be 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 be the ratio of immigration with a sufficient degree of immunity gained either by contact with COVID-19 or by effective vaccination. Then adding and 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 . Also, as the reader may easily check (we spare the details), the system (3.3) has a unique disease free equilibrium point , where In order to apply formula (3.2) we need to check whether conditions (1) through (6) are fulfilled. Before doing that, it is essential to show that (3.3) is positively invariant in that is, if , then for all .

The system(3.3)is positively invariant in Ω.

It suffices to show that for all whenever . Indeed, if positivity is proved, then by adding the equations in (3.3) we have , from which for all . If , then we would necessarily have that and hence . To prove the positivity, notice that if for some , then for for some by continuity of the functions . Assume now that for some with either , or . If , then either or . Using the second equation in (3.3) we arrive at Since the coefficients and are strictly positive, it follows that . This shows that for all . A similar argument proves that for all . Finally, if then the uniqueness of local solutions of (3.3) shows that must be the free of disease equilibrium point, and therefore for all . □

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

We just need to follow the steps introduced in Section 3.1 and check whether condition (1) through (6), shown on page 8, is satisfied. The system (3.3) can be written as in (3.1) by setting and On the other hand, the disease-free equillibrium point is where and .

We spare the details of the easy proofs of conditions (1) through (5) on page 8. As for condition (6), the disease-free system is given by By means of elementary integration one can check easily that all the solutions of the linear system in have the form for some . It follows straightforwardly that is the unique equilibrium point in 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 [43, Theorem 2] we have:

Consider the system(3.3), its unique disease free equilibrium point, whereandare as in(3.4), and its basic reproduction numbergiven by(3.5). We have:

Ifthenis locally asymptotically stable.

Ifthenis unstable.

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 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 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, , and with 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 ageing of the population. Hence an ageing 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 . 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 . Reasoning as we did in Section 3.2, the loss of population can be compensated with immigration. Let 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 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 , using elementary linear algebra we obtain the unique disease free equilibrium point of (3.8), namely , where

The local existence and uniqueness of the solutions of (3.8) is an easy consequence of the fact that all the functions 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:

The basic reproduction number for the disease free equilibrium point of (3.8) is given by where and are given by (3.9) and (3.10) respectively.

Following the steps pointed out in Section 3.1, we identify and where K is given by and Therefore, and which allows us to conclude The characteristic polynomial of the matrix 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 Observe that is In (3.14), adding the first and third equations on the one hand, and the second and the fourth on the other yields where and . Obviously, . 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 where and The explicit solutions of (3.16) can be expressed in terms of , and as for some . Since the eigenvalues of M, namely −δ and , are both negative, it follows that To finish we focus on the second and fourth lines in (3.14) and the solutions for and which, according to (3.17), are of the form where f and g are linear combinations of , and . Now, plugging (3.19) and (3.20) into (3.18), by elementary integration where F and G are again linear combinations of , and . Hence which concludes the proof. □

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

Consider the system(3.8), its unique disease free equilibrium point, defined in(3.9)through(3.12), and its basic reproduction numbergiven by(3.13). We have:

Ifthenis locally asymptotically stable.

Ifthenis unstable.

Calibration of the parameters

The parameters λ and () are easily obtained with accuracy from the ISTAT databases, and are shown in Table 2 . The constant μ is a weighted average of the '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 's and the '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 [30]).

Table 2

Average daily birth and death rates according to ISTAT (Istituto Nazionale di Statistica) in the period 2015-2019.

Parameter	Value
λ	1.9 × 10−6 births per person and day
λ1	5.05 × 10−5 births per person and day
μ1	9.89103 × 10−7 deaths per day and inhabitant
μ2	4.01955 × 10−6 deaths per day and inhabitant
μ3	3.27992 × 10−5 deaths per day and inhabitant
μ4	2.73280 × 10−4 deaths per day and inhabitant
μ	2.92106 × 10−5 deaths per day and inhabitant

As for β and we have analysed the data of daily and cumulative infections, deaths and recoveries provided by Protezione Civile Italia [34], and downloaded from [35]. The period used to calibrate β and consisted of the 30 days between the 16th of November 2020 () and the 15th of December 2020. Time will be measured in days. Consider the definitions where, on day t, are the new reported infections, are the active cases, are the susceptible and are the new reported deaths. During the calibration period 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 we subtract from N the cumulative cases between the onset of the pandemic and day t.

Observe that and , as they have been defined, are natural choices for and respectively. This justifies that and . It was shown in [24] that the observed values of during the first wave of COVID-19 in Italy decayed exponentially to a baseline value , 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 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 (and ) 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 with an exponential (see Fig. 4 ) where the coefficients a and b are determined by linear regression applied to the points with . Recall that is the 16th of November 2020. The choice maximizes the coefficient of determination of the linear regression and in this case On the other hand, unlike what happened at the onset of the pandemic, the 's are relatively stable since the end of the first surge (about June 2020). Their values oscillate around their mean (see Fig. 4). The value chosen for is exactly the average of the 's over the 30 day calibration period, that is,

Fig. 4

In the left we show the exponential decay of β(t) to β0 = 4.072 × 10−10 in Italy during the last 45 days of 2020. On the right, we can see how μ⁎(t) oscillates around their mean during the same period of observation.

As for the calculation of the death rates () caused by COVID-19, we have employed the weekly reports elaborated by the Istituto Superiore di Sanità (see [31]). These reports 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 is the Infection Fatality Ratio of the age group k, then the following relationships must hold for . Also, should be a weighted average of the 's. As a matter of fact we can put Observe that the weights show a very stable behaviour in the period prior to vaccination. Moreover, the approximation is legitimate (see Table 5 and Fig. 5 ), where is the number of people in the age group k and N the total population according to Table 1. Since the quotients can be estimated by (4.1) and has already been fixed, we arrive at the following choice of the 's:

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).

Age group	IFT in deaths by infection
0–39	0.000167463
40-59	0.002781909
60-79	0.041082514
80+	0.189311149

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.

Age group	% infections	% population
From 0 to 39	38.01	39.47
From 40 to 59	34.18	30.77
From 60 to 79	19.28	22.35
80 and older	8.53	7.41

Fig. 5

Relative weights of weekly new infections by age group.

The value of γ will be fixed in accordance with the observed daily recuperation ratios . On day t, is defined by where are the reported daily recuperations on day t. The average values of the 's are quite stable. We will assign to γ the average values of the 's in each of the periods shown in Table 4 .

Table 4

Calibration of γ. We consider the average values of the γ(t)'s in each period.

Period	Considered value for γ (recuperations per infected and day)
16/11/2020–15/12/2020	0.03192
16/12/2020–15/1/2021	0.03099
16/01/2021–15/02/2021	0.03594
15/02/2021–08/03/2021	0.03125

Validation of the model

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 [34], [35], 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 [37]. 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 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 [32]). 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: We also need to assign an initial value to the variables , , for . 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 , , for 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à [31]. Interestingly, the pattern of weekly infections since the beginning of November 2020 is quite stable, as the reader can see in Fig. 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 Fig. 6 .

Fig. 6

Assessing 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.

A model for β, and vaccines

The modeling of the second wave of COVID-19 in Italy was done in [24] satisfactorily using (1.1) by adapting the transmission function . Actually it is only that needs to be adapted in order to describe the behaviour 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 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.

A model for during the third wave

The authors of [24] have not observed any significant alteration of nor γ (not to mention λ nor ) 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, 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 . 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 (see [39]). A good example about how to adapt to model a coronavirus wave can be found in [24]. In this manuscript we suggest a plausible transmission function to model the third surge of COVID-19 in Italy. Unfortunately, the daily values of have been steadily increasing since mid February 2021 (see Fig. 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 its third surge of COVID-19. Fig. 7 also shows the proposed model for during the third wave. The explicit form of depends on parameters that we have defined in previous sections, in particular, we recall that Also, let m and r be, respectively, the slope of the linear regression applied to the points where t ranges over the 42 integers between 73 and 114, and . We have With all that in mind would be Observe that , , , and 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 with in Fig. 7.

Fig. 7

Proposed model for β(t) compared to β(t) until the first days of March 2021.

The choice of may appear whimsical at first sight. As a matter of fact, the function is just an educated guess on how the third surge of COVID-19 will evolve in Italy. Until the 9th of March 2021, fits the observed values of either by logarithmic regression (), or by linear regression (). For we just try to mimic the behaviour of during the second wave (see [24]) both in length and in intensity.

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 [36]), 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:

The efficacy of vaccines, and , for the first and second shots respectively.

The delay of the vaccine effect in providing immunity.

The proportion of vaccine-target citizens who do not belong to the recovered, that is, the value for , where represents the age group.

The number of vaccines administered to each group every day, namely, where will stand for the first or second dose.

Le us first make a reasonable choice of the efficacies and . Table 6 shows the efficacies of the main vaccines in use in Italy according to several studies [7], [27], [28], [29].

Table 6

Efficacy of the 3 main vaccines in use in the Italian vaccination programme.

Brand	θ_1	θ_2
Moderna	0.802	0.945
Pfizer-BioNTech	0.524	0.905
Oxford-AstraZeneca	0.641	0.704

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 [33]. 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.

Quarter	θ_1	θ_2
Q1	0.615	0.933
Q2	0.629	0.775
Q3	0.618	0.812
Average	0.622	0.824

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 , in other words, will be replaced by 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 , the recovered cases in group k and day t, into two subgroups where will stand for the recovered cases in group k and day t who had been previously infected and will stand for those susceptible ones who were vaccinated and succeeded in gaining immunity. With this notation, if is the amount of susceptible people on day and group k,

Finally, we have assigned to real data about daily vaccinations from the end of 2020 until the 4th of March 2021. We describe below how the functions have been defined in each of the two vaccination schemes under study.

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.

Vaccination plan B

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.

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 is the relative weight of group k among all active cases.

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 Fig. 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.

Fig. 8

Evolution of active cases with or without vaccination programmes.

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. Fig. 9 shows a significant reduction of mortality due to COVID-19 in both vaccination schemes with respect to the model with no vaccination programme. This reduction is about 44.29% for the 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 Fig. 9.

Fig. 9

Evolution of cumulative deaths with and without vaccination programmes.

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 programmes progress with normality. In particular, the cumulative incidence in 14 days per 100,000 inhabitants could reach values close to 1,000 (see Fig. 10 ).

Fig. 10

Evolution of cumulative incidence with and without vaccination programmes.

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 done in (6.1) and the vaccination rate of 300,000 shots per day considered above. If evolves differently or the vaccination capacity is improved, model (2.2) outcomes could be completely different.

Conclusions

Recall that the main objectives of this manuscript are the following: As can be seen in Fig. 8, Fig. 9, Fig. 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 scenarios described for vaccines. The only way to reduce this exceedingly high death toll is to act promptly to lower the rate of effective contacts combined with a higher rate of vaccination (in this paper we have used a capacity of 300,000 shots per day).

Assess the capacity of vaccines to hamper the spread of the virus and reduce the death toll.

Compare two different vaccination Schemes.

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.