Sensitivity theorems of a model of multiple imperfect vaccines for COVID-19.

In response to the ongoing pandemic of COVID-19, several companies across the world have proposed a wide variety of vaccines of different mechanisms of action. As a consequence, a new scenario of multiple imperfect vaccines against the SARS-CoV-2 arose. Mathematical modeling needs to consider this complex situation with different vaccines, some of them with two required doses. Using compartmental models we can simplify, simulate and most importantly, answer questions related to the development of the outbreak and the vaccination campaign. We present a model that addresses the current situation of COVID-19 and vaccination. Two important questions were considered in this paper: are more vaccines useful to reduce the spread of the coronavirus? How can we know if the vaccination campaign is sufficient? Two sensitivity criteria are helpful to answer these questions. The first criterion is the Multiple Vaccination Theorem, which indicates whether a vaccine is giving a positive or negative impact on the reproduction number. The second result (Insufficiency Theorem) provides a condition to answer the second question. Finally, we fitted the parameters with data and discussed the empirical results of six countries: Israel, Germany, the Czech Republic, Portugal, Italy, and Lithuania.

Introduction

The introduction of vaccines in the current pandemic of COVID-19 has been highlighted by many scholars as to the best opportunity to curtail the spread of the Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2) [1], [2]. In response to the outbreak, more than 50 organizations have started the development of different vaccines, and some of them have been finally approved and implemented around the globe [3].

Many compartmental models have been proposed since the early development of the pandemic [4], [5], [6], [7], [8], [9], [10], [11], [12]. Recent epidemiological models has added a compartment related to the vaccinated population [1], [13], [14], [15]. Several papers in this topic are mainly focused on finding the optimal allocation of the vaccines in different partitions of the population, varying from age [2], [3], [16], work [17], [18], the presence of comorbidities [19] and the social contact network [20].

The presence of multiple types of vaccines arises the necessity of considering multiple compartments in the modeling. This problem is addressed in [21], which presented the -ij-SVEIHQRD model including multiple vaccines and new variants. In addition, some articles have included specific compartments for one and two doses [3], [22], [23], which is the case of some types of vaccines.

SEIARD model (which stands for Susceptible-Exposed-Infected-Asymptomatic-Recovered-Death) has been proposed in the context of the ongoing pandemic of COVID-19, due to the particular transmission dynamics of this disease [24], [25], [26]. It is worth mentioning that a similar model (SEIADR) was studied in papers such as [27], which includes vaccination and antiviral treatment, but the Death compartment is considered as infectious, an aspect that fits in diseases such as Ebola but not in the case of COVID-19. In [22], the SEIARD model was extended to incorporate one vaccine compartment.

In this article, we extended the multiple imperfect vaccines preliminary model of [28] to develop a SEIARD-based system (instead of the original SIR-based model), to particularize it to the specific case of COVID-19, as an attempt to correct some of the previous limitations. Similarly, some questions around the availability of multiple vaccines are to be addressed in this paper. Rather than considering how to distribute the vaccines, we are more interested in the complete campaign itself, and the presence of a variety of types of vaccine.

Is the vaccination campaign enough to eliminate or mitigate the global pandemic? This question is the title of the paper of [14], which was one of the first researches in studying the effects of vaccination as a key measure to curtail the outbreak. In the case of the United States of America, the mentioned article reported that a vaccine with an efficacy of 0.8 needs 82% of the population to be vaccinated to reach herd immunity. However, the situation is more complex with different vaccines with distinct “efficacies”, and current data is required to answer this question a posteriori. Moreover, analytical conditions would be useful to state if the vaccine campaign is enough or a combination of non-pharmaceutical interventions is also necessary.

Another big question around this phenomenon is the following: is the addition of new types of vaccines in the campaign always useful to reduce the effects of the pandemic? Although, intuition might drive us to answer this question with yes, in the preliminary preprint, the Multiple Vaccination Theorem states some conditions where this addition of vaccines can increase the reproduction number. Does that theorem hold for a SEIARD-based model?

To answer the main section using a more specific model for COVID-19, we developed the model directly from SEIARD but taking the same procedure of the SIR-based system, and including compartments from one and two doses and it is presented in Section 2. In Section 3, we include all mathematical results, which are the basic propositions about the model (computation of the reproductive number, stability of the Disease-Free Equilibrium), and in the remaining section, we provide the sensitivity analysis of the reproduction number. An important result of the sensitivity is the mentioned Multiple Vaccine Theorem, which answers the second main question, whereas the first question is discussed in the Insufficiency Theorem.

Mathematical models

Mathematical models with imperfect vaccines consider different ways of vaccine failure, which are, the primary vaccine failure (probability of not reaching protection just after vaccination), the “leakiness” (probability of infection after exposure with an infectious individual), and the waning rate (rate of immunity loss) [29], [30], [31], [32], [33]. Most of the cited articles present a SVIR model focused on childhood vaccination, and thus, this situation might be adapted to consider the case of COVID-19.

Our model adds one vaccination compartment with its respectively full-vaccination compartment for each -th vaccine. Let the total number of types of vaccines. Thus, . For each type, the primary vaccine failure is , the waning rate is , the leakiness for the uncompleted vaccinated individuals (for instance, with one dose) is , and is the leakiness for complete vaccination. Additionally, and , the vaccination rate of the first and second dose respectively, complete the vaccination parameters. The set are the intrinsic vaccination parameters and the rates are extrinsic vaccination parameters because the intrinsic parameters cannot be changed but the rates depend on the campaign.

The SEIARD base model is taking from [25], with two differences: all the compartments are directly normalized such that . Another difference is the supression of two recovered population and . In this manner, . Transmission of disease happens from , , and (Susceptible) to (Exposed) with a force of infection , with and the transmission rates of symptomatic and asymptomatic individuals. The final model is given by the following equations:

A full description of the parameters and the compartments are provided in tables 1 and 2 . A flowdiagram of the model is given in Fig. 1 .

Table 1

Description of the compartments of the model.

Compartment	Description
V_j	Vaccinated individuals with the label j-vaccine
	with one dose and whose protection has not wanned.
F_j	Vaccinated individuals with the label j-vaccine
	with two doses and whose protection has not wanned.
S	Susceptible individuals who have never been vaccinated
	or had experienced vaccine failure
E	Exposed individuals (in the latend stage)
I	Infected symptomatic individuals
A	Infected (and infectious) asymtomatic individuals
R	Recovered individuals
D	Dead individuals by COVID-19

Table 2

Description of the parameters of the model and their range.

Parameter	Range	Description
β_1	[0,∞)	Transmission rate from infected individuals
β_2	[0,∞)	Transmission rate from asymptomatic individuals
σ	[0,∞)	Incubation rate
γ	[0,∞)	Recovery rate
δ	[0,∞)	Death rate
w	[0,1]	Fraction of symptomatic cases
ε_j,A	[0,1)	Primary vaccine failure
ε_j,L1	[0,1)	Leakiness of the first dose of the j-th vaccine
ε_j,L2	[0,1)	Leakiness of the second dose of the j-th vaccine
α_j	(0,1]	Waning rate of the j-th vaccine
r_j	[0,1]	Vaccination rate of the first
		dose of the j-th vaccine.
λ_j	[0,1]	Vaccination rate of the second
		dose of the j-th vaccine.

Fig. 1

Schematic description of the model. represents the force of infection.

Basic results and stability

The Multiple Vaccination Model satisfies the following properties:

Let The system of equations is invariant in the set .

The Disease-Free Equilibrium is given by

The Basic Reproduction Number associated with the model with vaccination is given by

where , is the Basic Reproduction Number of a SEIARD model without vaccination.

The proof of 1.1 will be omitted, since it is standard. In the case of item 1.2, we can obtain (9) and (10) by setting and . The final equation is derived from item 1.1: .

The final item (1.3) uses the Next Generation Matrix technique [34]. Let . Thus,

Therefore, if is the DFE and ,

Thus,

On the other hand,

The spectral radious of is the reproduction number, which is where  □

An important quantity related with the reproduction number is the vaccine impact, represented with . In [33], this number is defined as “the relative reduction within a population per unit of vaccination”. We will address this definition by parts. Since, , then

Thus, the factor represents the reduction on the reproductive number. Dividing this quantity by the vaccination rate we get

If, the Disease-Free Equilibrium is locally stable. On the contrary, if, the DFE is unstable.

The linearization matrix of the system evaluated in the DFE is (without the final row) where .

The characteristic polynomial of the previous matrix is where is the following matrix

Using elemental operations where . Since is a triangular matrix, which yields the solutions of the equation : for , , for and . Those solutions are also zeros for . Finally, the second factor of can be expressed as where

By Routh-Hurwitz criterion, all the solutions of (25) have a negative real part provided that and . Since , , and are terms of , holds. is a direct consequence of , which implies and therefore . Note that if , the eigenvalues are positive. This procedure, yields the existence of solutions , , and with negative real-part. Lastly, . This shows that the DFE is stable if , and unstable if . □

Methodology for parameter estimation

Parameter estimation is the main topic of this section. Since some countries have provided data about their vaccination campaign, parameter optimization is possible for the model, although it relies on several parameters. Data was collected from the John Hopkins repository [35] and the source of Our World in Data [36], which adds a specific dataset for vaccination by the manufacturer. Data only provides the fraction of the vaccinated individuals with two and one doses in general, but the information about the number of individuals vaccinated by type (manufacturer) and the number of doses is not given and was estimated based on the available numbers. This situation adds some error in data, and it is shown in the non-strictly increasing amount of fully vaccinated individuals by type. Vaccine intrinsic parameters are provided by the extensive research on the vaccines (see Table 3 ). The incubation rate is set to following [37], and was set to 0.25.

Table 3

Parameters of the model. The waning rate remains unknown for all the considered cases. Nevertheless, since estimations are considering a window of fewer than six months, in this paper we assumed for all .

j	Vaccine	1−ε_j,L1	ε_j,A
	Manufacturer	1−ε_j,L2	λ_j
1	BNT162B2	0.52 (0.295,0.684)	3.77×10^−4[40]
	Pfizer	0.95 (0.903,0.976) [41]	1/21[42]
2	mRNA-1273	0.8 [43]	7.23×10^−4[43]
	Moderna	0.941 (0.893,0.968)	1/21[44]
3	AZD1222	0.463 (0.318,0.578) [45]	0.00268 [40]
	AstraZeneca	0.76 (0.68,0.82)	1/84[46]
4	Ad26.COV2.S	0.663 (0.599,0.718) [47]	0.0091 [47]
	Johnson&Johnson	(Single-dosed vaccine)	1 [48]

In general, this procedure will be done by using gradient-based methods. These algorithms are usually difficult to implement, since computing the exact gradient requires the analytical solution of the model, which is remarkably difficult in SIR-based systems [38]. Some authors, such as [39], fit the data with a polynomial and define a new loss with the derivatives of the polynomial and the model. In this case, we will use finite differences instead to approximate the gradient of the loss function Sum of Square Errors.

For some cases, we used a variation of Gradient Descend using finite differences and the backtracking algorithm. Instead of using a fixed learning , we used a vector , and hence, the Gradient Descend is given by where represents the pointwise product of the two vectors. This change was introduced to consider the range of parameters: for instance: is usually taken in a space whereas other parameters are restricted in [0,1]. Varying at the same rate as the rest of the parameters would only slow the process, and in practice, it could stack easily into a local minimum.

In the remaining cases, the optimization routine was done separately. For instance, (or its associate parameters, see the following subsection) can be optimized by using the “reported compartment” solely and the reported susceptible. However, at least the rates and need to be fitted with all the data.

Another problem related to this technique is that data is not fully provided in the sense of the model. A close look at the compartments shows that some of them need to be redefined in order to make a proper estimation. This is the case of compartments , , and , which will be redefined as . We can redefine the equations to consider to implement the parameter estimation

The changes are motivated for various reasons. Reported vaccinated, for example, is an increasing function and it does not consider failure. Nevertheless, we can define the dynamics of the reported sharing parameters with the original model, which is the rate in this case. The situation of is a bit complex: represents the fraction of reported vaccinated individuals who has not been fully vaccinated. In this sense, are individuals that have been vaccinated at least one, and is the compartment of vaccinated individuals whose vaccine has not failed. This establishes the difference between the two groups.

Finally, are reported individuals that have been infected. This definition presents a controversy: reported recovered might have been asymptomatic individuals or not. We need to compare one compartment with the data, and assumes that all the recovered that have been reported have been infected previously, whereas the original compartment takes into account the asymptomatic.

Time-dependent parameters

As mentioned in [28], one common problem for parameter optimization in the case of vaccination models was the fact that the vaccination rate is not constant. As a response, in this article, we used time-dependent parameters for the vaccination rates of the first dose. Time-dependent parameters have been used in different compartmental models of COVID-19. For instance, Caccavo [8] used exponential and logistic functions to model the general parameters of the pandemic in the case of Italy and China (transmission, death, and recovery rates). Not always the addition of time-dependent parameters is useful to model a particular situation (see [49]). The simplest format is the constant function, but in some cases a change on the vaccination rate is induced. A simple model for this discrete change is given by the piecewise function:

Another possibility is the exponential growth of the vaccination, given by the expression:

A combination between the exponential and piecewise constant function is given by the piecewise exponential function:

These exponential functions might not be the best options for inferences when , and the vaccination rates might be, instead, logistic functions:

Transmission rates and were originally considered fixed but in some cases was ideal to introduce some variations, following the structure of piecewise functions such as for . This function was used in the parameter estimation of Italy and Lithuania, whereas in the case of Germany and Czechia the function was set constant but we analyzed the possibility of introducing a modification.

Results

Description by country

Six countries were selected to perform parameter optimization on their publicly available data: Israel, Germany, Czech Republic, Portugal, Italy, and Lithuania. Two strategies are used by the elected countries: one single type of vaccine (Israel) or four vaccines (European countries), which are BNT162B2 (Pfizer), mRNA-1273 (Moderna), AZD1222 (AstraZeneca), and Ad26.COV2.S (Johnson & Johnson). The vaccines are indexed from 1 to 4 following the given order. The first three vaccines require two doses, but Ad26.COV2.S needs only one. To address this dissimilarity we set . In addition, the temporal windows were selected according to the last peak of infection. Wider windows are possible, but it needs to model with more breakpoints.

The majority of the parameters were optimized separately and only and were part of the global optimization. This technique separates the global loss function into smaller loss functions avoiding an unnecessary increment of the local minima of the global loss. Main results of the parameter optimization are provided in Table 4 , in the case of European countries, and in 7 in the case of Israel. Vaccine parameters of European countries are given in Table 5. Confidence intervals of and are provided (see also Table 8), and were constructed based on the confidence intervales of the vaccination parameters. Figures from 2 to 9, and 11 summarize the estimated dynamics of the model per country. To see simulations of the future behavior of , and possible scenarios, see Figs. 6, and the subsection of Variation of parameters.

Table 4

Parameters of European coutries. In the case of piecewise transmission rates, stands for and stands for . TW stands for temporal window.

Aspect	Germany	Czechia	Portugal	Italy	Lithuania
Final date	2021-07-15	2021-07-16	2021-07-23	2021-07-21	2021-07-18
TW	76	121	172	112	67
R_0	0.7962	0.5242	3.0586	0.8221	1.8953
CI upper	0.3809	0.2165	1.2018	0.3419	1.079
R_c	0.3453	0.1883	1.1061	0.2962	0.9905
CI lower	0.3379	0.1853	0.9491	0.2891	0.959
β_1	0.102524	0.046501	0.047020	0.001552	0.06
β_12	x	x	0.056424	0.007262	0.017218
β_13	x	x	0.14896	x	x
β_2	0.045884	0.032254	0.018513	0.030266	0.008703
β_22	x	x	0.067708	0.039335	0.115664
β_23	x	x	0.124133	x	x
γ	0.075026	0.068011	0.046517	0.061499	0.048019
δ	0.000924	0.000976	0.000725	0.000611	0.000483
w	0.25	0.25	0.25	0.25	0.25

Table 7

Parameters and analytical values of Israel .

β_11	0.062029
β_12	0.006911
β_21	1.3957
β_22	0.1555
γ	0.098765
δ	0.000652
w	0.25
r_0	0.011797
r_1	0.002063
θ_r	50
R_0	4.6905
R_c	3.6367
S^0	0.707

Table 5

Analytical values per country and per vaccine. In the case of Lithuania, stands for .

Parameter	Germany	Czechia	Portugal	Italy	Lithuania
	Pfizer
r_1	0.005912	x	x	x	x
r_10	x	80	120	80	40
r_11	x	0.014237	0.012465	0.015479	0.008476
r_12	x	0.027476	0.026379	0.024412	−0.004783
	Moderna
r_2	0.000814	x	x	x	x
r_20	x	100	100	30	30
r_21	x	0.001326	0.0016	0.001139	0.001131
r_22	x	0.012585	0.041129	0.055246	−0.018878
	AstraZeneca
r_3	0.001007	x	x	x	x
r_30	x	100	100	60	30
r_31	x	0.001307	0.00253	0.002096	0.002252
r_32	x	0.012893	0.02186	−0.002049	0.009902
	Johnson & Johnson
r_40	0	7.67×10^−5	0.000333	0.000412	50
r_41	0.001016	0.013872	0.007533	−7.27×10^−5	0.002237
θ_r_4	20	40	85	20	−0.005818

Table 8

Values of , , and per country and per vaccine.

Aspect	Israel	Germany	Czechia	Portugal	Italy	Lithuania
CI upper	0.7885	0.4785	0.4131	0.3929	0.4159	0.5691
η	0.7754	0.4337	0.3593	0.3616	0.3603	0.5226
CI lower	0.7677	0.4245	0.3533	0.3103	0.3516	0.506
	Pfizer
ϵ^(1)	0.0949	0.0949	0.0949	0.0949	0.0949	0.0949
η_1	x	0.6835	0.7596	0.5043	0.7620	0.734
	Moderna
ϵ^(2)	0.0737	0.0737	0.0737	0.0737	0.0737	0.0737
η_2	x	0.4560	0.3718	0.3426	0.388	0.5373
	AstraZeneca
ϵ^(3)	0.3345	0.3345	0.3345	0.3345	0.3345	0.3345
η_3	x	0.4441	0.3607	0.3231	0.3709	0.5544
	Johnson & Johnson
ϵ^(4)	0.3351	0.3351	0.3351	0.3351	0.3351	0.3351
η_4	x	0.4638	0.3667	0.3631	0.3608	0.5387

Fig. 2

Parameter estimation of Germany. The dashed lines stand for the estimated curves of reported compartments (), whereas the solid lines represents the original definition of the compartments () and the dots represent the data points. The green curve of recovered individuals is . (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 9

Parameter estimation of Lithuania. The dashed lines stand for the estimated curves of reported compartments (), whereas the solid lines represents the original definition of the compartments () and the dots represent the data points. The green curve of recovered individuals is . (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 11

Dynamics of the model in the case of Israel. The dashed lines stand for the estimated curves of reported compartments (), whereas the solid lines representes the original definition of the compartments . The green curve of recovered individuals is . (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 6

Estimated dynamics of the behaviour of in future days.

In the following lines, we will describe the situation of the European countries and the parameter optimization. The case of Israel will be described in the subsection of Reduction of vaccine efficacies, because further considerations are needed to estimate its parameters.

Germany

Germany presented a reduction in the cases from April to June, but this situation seems to change in July. The first three vaccines were set constant, but in the fourth vaccine (Johnson & Johnson), we selected an exponential piecewise function. The fact that , implies that solutions are approaching the Disease-Free Equilibrium, provided that the parameters will remain constant, including the vaccination rates. A reduction in the number of infected individuals is observed from the data. Nevertheless, this reduction is slow, and a projection shows that ( the number of inhabitants of the country) only for , implying that only 255 days after the final date of estimations, the model forecast the end of the pandemic on Germany if the conditions do not change. However, as we shall see, it is difficult to believe that parameters are not changed in the case of Germany. It is also important to mention that , meaning that as approaches infinity, the fraction of vaccinated individuals whose protection has not failed would be almost 65%.

In general, at the time of the parameter estimation, the number of cases declined, varying from 22,262 (April 26) to less than 1000 in the first half of July. An increment in recent days is observed, but it is difficult to decide whether this local tendency would change the global slope of the curve. Around 58.6% of the total population has been vaccinated at least with one dose, and it is estimated that around 5.11% present vaccination failure.

Czech republic

In the case of Czechia (Fig. 4), the dynamics of the vaccination rates of the first three vaccines were logistic, and in the case of the Johnson & Johnson vaccine, we select the exponential piecewise model. The reduction of new cases started on the first days of March. On March 3, 16,816 new cases were reported, but during the last 30 days, less than 400 new positives appear by day. If this tendency does not change, in around 143 days the pandemic would be over. Since , if the vaccination campaign does not change its rate, it is expected that around 72% of the population would be vaccinated (with protection), and the rest of the susceptible. 49.59% of the Czechs have been vaccinated, but 43.58% have not experienced failure. However, a small increment is observed in the last days, and therefore, changes in the parameters must be analyzed.

Fig. 4

Parameter estimation of Czech Republic. The dashed lines stand for the estimated curves of reported compartments (), whereas the solid lines represents the original definition of the compartments () and the dots represent the data points. The green curve of recovered individuals is . (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

In order to select the proper functions for the vaccination rates, comparison between the functions were performed, however, if a simpler function fitted correctly with the data, no improvement was included. For instance, in the case of Germany, constant functions were sufficient, but this was not true for Czechia. The estimation error using a constant function for was 0.18534, whereas the estimation error with the logistic function was 0.00126. In Fig. 3 this difference is shown visually.

Fig. 3

Comparison between the models of . In the presented case, parameters of (green line, which is also the parameters of ) are fitted with the data of vaccination in the Czech Republic (red dots). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Portugal

Portugal (Fig. 5) has particular dynamics on the development of its pandemic. In order to estimate correctly the behaviour of , it was needed to add three changes on the transmission parameters and , following the function for . In this case, , and , specifying the changes on the dynamics of the transmission rates. This behavior can be summarized as follows: in January 2021, Portugal reaches its peak of new cases with 16,432 new reported cases in one single day (January 28). In February, a rapid descend of the cases was presented, followed by a period of a relatively low number of cases (around 500 new daily cases), which was ended in the last two months (June and July), when a new increment is presented. In the next month, we forecast that the cases would continue growing until 37 days, reaching a peak of 71,770 active infected. On the other hand, the vaccine parameters of the first three vaccines were logistic, but in the case of the last vaccine, we used the piecewise exponential function.

Fig. 5

Parameter estimation of Portugal. The dashed lines stand for the estimated curves of reported compartments (), whereas the solid lines represents the original definition of the compartments () and the dots represent the data points. The green curve of recovered individuals is . (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

64.01% of the total population has been vaccinated, but only 47.66% are estimated to be in the vaccinated compartments. It was also found that , which means the DFE includes roughly 75% of the total population in the vaccination compartments.

Italy

Italy (Fig. 7) has experienced a reduction in the cases until the recent days when a change in the behavior of both active and new cases is observed. On March 13, Italy suffered a peak of new daily cases with a total of 26031. Then, the number of cases declined but nowadays the situation has been reverted. In this case, the first three vaccines rates were modeled using a logistic function, and the fourth vaccine was modeled by using the piecewise exponential function.

Fig. 7

Parameter estimation of Italy. The dashed lines stand for the estimated curves of reported compartments (), whereas the solid lines represents the original definition of the compartments () and the dots represent the data points. The green curve of recovered individuals is . (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

The value , meaning that more than 70% of the population would be in the vaccination compartments in the DFE. Around 60.82% of the population has been vaccinated currently, but only 53.8% have not experienced any type of vaccine failure.

Lithuania

In the case of Lithuania, all the vaccination rates were set as logistic functions, since the alternative options did not fit appropriately with the data. On April 29, this country suffered the last recent peak of new cases with a total of 1341. A reduction is observed in the following days reaching a minimum of around 15 new cases. However, in recent days, a new increment has surpassed 100 new cases, and this situation might be considered relevant. Therefore, we added a piecewise function in the transmission rates using a threshold of 45, which provides a better estimation of the function as seen in Fig. 8 . This change makes it possible to see details of the reproduction number in the recent days, which is estimated as , which implies that a recent modification of the dynamics of the transmission rate occurred. The confidence interval, however, shows that it is not possible to assert that . Moreover, even if the outbreak is controlled, more than 3000 days are needed to finish it.

Fig. 8

Comparison between two models of the behaviour of the transmission rates and their effect on . In the first model (A), and are constant whereas in the second model are piecewise functions.

At the moment of the parameter estimation, 49.09% of the population have received at least one dose, and 48.48% have not experienced vaccine failure. The value of implies that less than 60 & would remain in the vaccination compartments if the DFE is meet.

Sensitivity analysis of

In order to study the influence of each parameter on the reproduction number, we can use the normalized forward sensitivity index, introduced in [50]. This index, also called elasticity index [51], measures the relative change of the reproduction number when a specific parameter is changed.

If is differenciable with respecto to a parameter , the normalized forward sensitivity index of is given by

Using the above definition, we can calculate the forward sensitivity of the parameters

Note that are always positive and are always negative. The fact that the death rate has a negativity sensitivity should not be unexpected, since a very deadly infectious disease cannot spread easily. Just to provide one example, the MERS, an infection caused by another coronavirus, has a higher death rate, but a lower reproduction number [52].

Nevertheless, a natural assumption is to expect , which yields , indicating that if is larger than , the recovery rate contributes more to the reduction of the reproductive number. Another possibility is that the transmission rate from infected is lower than the transmision rate from asymptomatic , in such case, any increment on the percentage of asymptomatic would increase the reproduction number. In general, this situation happens if and only if .

In the case of vaccine parameters, let . Thus, and

Index of vaccine parameters are always positive, which means that they must be as lower as possible. A much complex situation arises when analyzing key parameters , and , since, against common sense, they do not show a clear tendency. This topic will be covered in the subsection of Multiple Vaccination Criteria, as part of the Multiple Vaccination Theorem.

Sensitivity analysis of the selected countries

The sensitivity analysis of the studied countries presents the values of elasticity of some relevant parameters in Table 6 . In the case of European countries, we also provide a bar chart of the elasticity values of (Fig. 10 ). We can see that the Pfizer vaccine has more negative value of elasticity for all the analyzed countries, whereas the other vaccines are closer to zero. In Portugal, Italy, and Czech Republic, the Moderna vaccina has second most negative elasticity, while in Germany this is the case of Johnson & Johnson, and in Lithuania is the case of AstraZeneca. Nevertheless, in terms of the sensitivity itselft (partial derivative), most negative values are observed in Johnson & Johnson vaccine. This dissimilarity might be explained by the higher vaccination rates of Pfizer vaccines.

Table 6

Sensitivity analysis per country.

Aspect	Israel	Germany	Czechia	Portugal	Italy	Lithuania
Γβ_1R_c	0.7482	0.4239	0.3214	0.2148	0.0571	0.0468
Γβ_2R_c	0.2518	0.5761	0.6785	0.7852	0.9429	0.9532
ΓγR_c	−0.9951	−0.9948	−0.9954	−0.9967	−0.9991	−0.9995
ΓδR_c	−0.0049	−0.0052	−0.0045	−0.0033	−0.0009	−0.0005
ΓwR_c	0.6643	0.2318	0.0953	−0.0469	−0.2572	−0.2709
	Pfizer
∂R_c∂r_1	−393.168	−19.36	−7.8152	−29.7082	−12.5859	−67.4425
Γr_1R_c	−0.2048	−0.3315	−0.4431	−0.2656	−0.4511	−0.27064
	Moderna
∂R_c∂r_2	x	−20.5641	−8.4389	−32.9065	−13.5288	−70.7599
Γr_2R_c	x	−0.0485	−0.0334	−0.0452	−0.0514	−0.027
	AstraZeneca
∂R_c∂r_3	x	−7.4276	−0.9592	8.6145	−2.4558	−38.8745
Γr_3R_c	x	−0.0217	−0.0037	0.0163	−0.0082	−0.052
	Johnson & Johnson
∂R_c∂r_3	x	−22.029	−9.4323	−39.0995	−14.949	−73.3008
Γr_4R_c	x	−0.0649	−0.0203	−0.0426	−0.0206	−1×10^−23

Fig. 10

Values of per country of European countries.

An interesting situation is given in the positive elasticity of the AstraZeneca vaccine of Portugal. This seems to be counterintuitive, but we shall see that this rare situation can be explained. In addition, it is worth mentioning that the recovery rate provides the most negative elasticity in all cases, with a value of almost -1.

Variation of parameters

In this subsection, we study variations of vaccine efficacies and transmission rates and , motivated by the fact that the countries are showing a change in the behavior of the outbreak and the emergence of new variants of the SARS-CoV-2, and their effect on the vaccine efficacies. Evidence of this phenomenon is presented in several papers such as [53] in Alpha and Beta variants and [54] for Epsilon variant. One of the first countries to present a clear increment on the daily cases in Israel. To model this particular situation, we reduced by a factor of 0.33 and by a factor of 0.88 when . It was also needed to vary the transmission rates by a factor of 22.5, although both changes affect increasing the number of infected, is difficult to decide which change explains better the increment of the cases.

In addition, we studied the situation of other countries (Germany, Czech Republic, Italy, Lithuania) by analyzing different panoramas of parameter variation. Portugal was not included in this analysis, because we performed a parameter estimation by introducing changes in the behavior of the transmission rates, and that was not the case for the remaining countries.

The variation of parameters can also be used to determine how can we increment the vaccination rates, to reach . In particular, it was observed that in the case of Portugal and Israel. In the first situation, if we increase the vaccination rates by a factor of 1.26, we can obtain . However, this change must be immediate, which is an unpractical measure. Another possibility is to increase the vaccination rate by a factor of 1.15 and in one month (30 days), we can reach the value of . In total, this involves getting a covering of 90.9% of the total population with at least one dose, and 79.75% with two doses. The case of Israel is a bit dramatic, because . This implies that the increments on the vaccination rate are unable to reach . In this case, reduction of the transmission rates seems to be the most effective option to reduce the impact of the pandemic.

The case of Israel

Israel is a remarkable case where only one type of vaccination was primarily used, which is the Pfizer (BNT162B2) vaccine [55]. This situation induces the presence of an authentic VSIR model. The final considered date was July 25, 2021, and the temporal window of the analysis consists of 194 days. The vaccination parameter was time-dependent, using a piecewise function. Moreover, another piecewise function was added to model the dynamics of the transmission rates, since in an increment on (set to ) and (set to ) was implemented to fit properly with data. All the parameters and some analytical values are presented in Table 7 and the estimated dynamics are given in Fig. 11.

In general, this change in the transmission rates has some consequences on the behavior of , which is increasing. We estimate that a new peak might occur in 50 days, taken 309,212 infected (see Fig. 12 ). 63.6% of the population have received at least one dose, whereas 58.4% have received two doses, but only 33.81% of the population belongs to the vaccine compartments.

Fig. 12

Estimated dynamics of in the following days.

Impact of the vaccination coverage on

The Eq. 20 allow us to restate the reproduction number in terms of the final values of , in the case , which is the situation of Israel. A complex panorama appears when , but for visuallization purposes, we will consider only this case. Let , . Therefore,

can be understood as the control reproduction number as a function of the final proportions of the vaccinated population with one and two doses respectively. These proportions must satisfy the constraint .

In Fig. 13 , three contour plots of are depicted, considering the confidence interval of the vaccine with the reduction of efficacies established at the beginning of this subsection of Israel. The isosceles triangle in the right upper corner is represented with the value 0 but is not included in the domain . Higher proportions of seems to be better to reach the value (more than 80%, with 20% of people vaccinated with one dose). Nevertheless, it is important to mention that in this specific case, is particularly high and increasing the vaccination rate solely is unable to find , as shown previously.

Fig. 13

Three different contour plots of . The values of the lower bound are , the values of the baseline case are , and the values of the upper bound are .

Other countries

Germany

An interesting situation is modeled when we change the values of vaccine efficacies, and double the values of the transmission rates and , in order to observe the reduction of efficacies with incoming SARS-CoV-2 variants. Four different panoramas are expressed in Fig. 14 :

Fig. 14

Four different scenarios of variation of parameters in the case of Germany, and its effect on .

Case 0: No change on the parameters.

Case 1: Duplication on the transmission rates.

Case 2: Duplication on the transmission rates and reduction on the vaccine efficacies by rate 0.75.

Case 3: Duplication on the transmission rates and reduction on the vaccine efficacies by rate 0.5.

Despite the parameter fitting, data seems to tend to the case 1, which yields and , whereas in case 2, , and in case 3 yields . If the transmission parameters are not duplicated, for instance, multiplied by a factor of 1.5, , and in the basic case. This means that the local tendency of Germany can surpass the global tendency in the following dates of the parameter estimation, and in case 3, the Disease-Free Equilibrium is unstable.

Czech Republic

The three considered cases were

Case 0: No change on the parameters.

Case 1: Duplication on the transmission rates. , .

Case 2: Duplication on the transmission rates and reduction on the efficacies by 0.5. .

Case 3: Triplication on the transmission rates and reduction on the efficacies by 0.5. , .

In contrast to Germany, Czech Republic required time-dependent vaccine parameters, as stated previously, specifically the logistic function. This particular situation motivates to study of the variation of the function across time and influences the behavior of . In the four cases it was found that . A consequence of this limit is observed in Fig. 15 , when all the four cases yield the same panorama.

Fig. 15

Four different scenarios of variation of parameters in the case of Czech Republic, and its effect on .

Italy

In order to simulate possible changes on the behaviour of the parameters, we considered the following panoramas:

Case 0: No chance on the parameters.

Case 1: Duplication on the transmission rates. , .

Case 2: Duplication on the transmission rates. , .

Case 3: Triplication on the transmission rates. , .

Case 4: Triplication on the transmission rates and reduction of the vaccine efficacies by rate 0.5. , .

Simulation of the referred panoramas is presented in Fig. 16 . Case 3 seems to be more related to the local trend of the data, which means that the parameters might have already changed.

Fig. 16

Four different scenarios of variation of parameters in the case of Italy, and its effect on .

Lithuania

The following cases were considered to simulate different panoramas:

Case 0: No change on the parameters.

Case 1: Reduction on the vaccine efficacies by 0.5. .

Case 2: Reduction on the vaccine efficacies by 0.75. .

Case 3: Increment on the transmission rates by a factor of 1.5. , .

Case 4: Increment on the transmission rates by a factor of 1.5, and reduction on the vaccine efficacies by a factor of 0.75. , .

All the dynamics of the five cases are presented in Fig. 17 . These cases seem consistent with the tendency a shows interesting behavior: a peak is reached in the following days and then it decays. This pattern could be relevant since the dynamics of the infected in most countries have the same tendency, and in this case, could be explained as an increment in the transmission rates, reduction in the vaccine efficacies, or a combination of both events.

Fig. 17

Five different scenarios of variation of parameters in the case of Czech Republic, and its effect on .

Multiple vaccination criteria

Multiple vaccination theorem

In this section, we will continue analyzing the sensitivity of key parameters for the vaccination, which are and . The importance of the mentioned parameters relies on the fact that they might be changed by governmental measures in terms of the addition and selection of vaccines. Other vaccine parameters are subject to the vaccine itself and not to the campaign. For that reason, they are fixed and cannot be used to design a better strategy, but are important to consider in the selection of the types of vaccine.

A generalization of the so-called Multiple Vaccination Theorem need the following remark: the Reproduction Number can be rewritten as

Multiple Vaccination Theorem

Letthe reproduction number of the model withvaccines,the reproduction number withvaccines.

Let

Let. The following properties are satisfied:

is constant over any chances of the vaccine parameters of the-th vaccine.

If,is a decreasing function and.

If,is constant and.

If,is an increasing function and.

Let , , . Note that and

Therefore, which does not include any vaccine parameter of the -th vaccine. In addition, the line given by equals provided , and if , , and if , . Finally, recall that where

Since , then which implies that if , , then ; if , ; and if , , indicating the behavior of . □

In the case of , a similar criterion holds for the function , but it is weaker since the threshold is not a constant. if and only if

Nevertheless, selection of better also depends on epidemiological criteria. Another similar condition holds for : if and only if -analysis

The Multiple Vaccination Theorem provides tools to analyze the effect of one particular vaccine on the reproduction number. This is particularly useful for the model with multiple vaccines. First of all, we need to compute . Second, we wonder what is the value of is the model lacks a particular vaccine? This question motivates the definition of : the value of computed without the parameters of the -th vaccine. More specific,

This definition is particularly useful because, as a corollary of the Multiple Vaccination Theorem, if and only if the addition of this vaccine is justified since it decreases .

Insufficiency theorem

An important consequence of the Multiple Vaccination Theorem is the Insufficiency Theorem, a criterion that establishes the condition of the existence of a set of values such that . In such a case, the DFE is unstable (Proposition 2).

Hence, the Insufficiency Theorem ensures the condition where the vaccination is unsufficient to end the pandemic. If this happens, the outbreak is likely to persist even if the vaccination rates are increased (and can also be theoretically infinite). Other public measures would be required to reduce the reproduction number, focusing on minimizing the transmission rates and . Firstly, we need to prove the following basic result:

Letand. Ifexists, thenexists, and the following relationship holds

Let and . Thus, exists , such that implies . If , then and , which yields . □

For the Insufficiency Theorem, the notion of the best vaccine is needed to use it to evaluate the possibility of reaching with the extrinsic vaccination parameters, in particular, . This notion is derived from the following definition.

A vaccine with label is said to be with the best better parameters than the -vaccine if for each vaccine parameter (excluding the rates), . Vaccines are said to be equivalent if for all , .

The sensitivity analysis gives prior conditions about when the parameters of a vaccine are better. For instance, less values of are better, but a more complex situation is observed with . The next proposition permits to study of a lower bound for inequivalent vaccines, and the limit when all the rates tend to infinite.

For equivalent vaccines,

This a straightful computation. Using Lemma 1, which implies  □

As a consequence, for equivalents vaccines, there exist the possibility of not reaching the DFE, depending on the parameters of the model. This includes , which means that for a greater reproductive number without vaccination, the combined leakiness must be set as lower as possible to prevent the undesired situation of unstability of the DFE. If a best vaccine exists, we can generalize the last result in a form of criteria.

Insufficiency Theorem

Let us consider the existence of the vaccine with the best parameters (labeled with). The inequality

holds if and only if for any values.

The converse is direct. If the remaining vaccines are equivalent,

The Multiple Vaccination Theorem ensures that any restriction is a monotonic function. Without loss of generality, let us consider the first vaccines (labeled with ) such that . Therefore, .

Since is the index of the vaccine with better parameters, any change of the parameters with respect of the model with equivalent vaccines, yields a higher reproduction number. which proves that . □

In other words, the Insufficiency Theorem states that if the best vaccine exists, the reproduction number is greater or equal than one if and only if the best vaccine with a theoretically infinite vaccination rate is at least one. If there is no a “best vaccine”, a weak result can be stated if we construct an idealistic best vaccine choosing the best parameters. In such a case, the converse of the theorem does not hold but provides us a condition when the desired cannot be reached. Due to its importance, the quantity will be named as insufficiency parameter, and this name will be refered in the empirical results.

Application of the multiple vaccine criteria

Both the Multiple Vaccination Theorem and Insufficiency Theorem can be used for specific purposes in this study. For instance, the -analysis is an alternative method to evaluate the effect on a particular vaccine. In this case, for European countries (where multiple vaccines were applied) we can see that the inequality holds for most cases, except for (AstraZeneca) in Italy, which is confirmed by the positive value in both sensitivity and elasticity, which means that the particular type is not helping to reduce the reproduction number.

A direct application of the Multiple Vaccination Theorem can be found in the case of Israel. In this country, since (Moderna), (AstraZeneca), and (Johnson & Johnson), whilst , the Multiple Vaccination Theorem do not advises agaisnt the application of the last two vaccines.

The Insufficiency Theorem is also useful in this context, in particular when . This is observed in the analysis of Portugal, where , but the Insufficiency parameter is 0.2255. In other countries, the Insufficiency Theorem can be used for the cases of variation of parameters. This is Case 3 of Germany, when the Insufficiency parameter is 0.2137 (using Moderna vaccine), which implies that even in the most dramatic and considered situation, it is possible to reach . This situation is also observed in Case 4 of Lithuania because the insufficiency parameter is lower than 1, (0.8679), which means that even in a dramatic situation vaccines could still be useful.

Discussion

The “Sensitivity Theorems” (Multiple Vaccination and Insufficiency Theorems) are a useful tool to analyze the behavior of a system of multiple (and different) types of vaccines. The importance of the results relies on the fact that answers two important questions about the multiple vaccination system: What is the effect of one particular type of vaccine on the reproductive number? Is the vaccination campaign enough to curtail an outbreak? The first question is the topic of the Multiple Vaccination Theorem whereas the second one is covered in the Insufficiency Theorem.

In practical terms, the Multiple Vaccination Theorem is a useful criterion to verify if a particular vaccine is decreasing (or increasing) the reproduction number. The empirical results showed that, indeed, vaccines with higher leakiness might not be appropriate for the vaccination campaigns of specific countries. The fact that in at least one country (Portugal), all vaccines satisfy the needed condition for reduction of , indicates that vaccines are not optimally distributed across the world to reduce contagion. Wider studies across the planet are needed to prepare better planning in the global distribution of vaccines. Nevertheless, this approach has only considered vaccination in a theoretical (mathematical) view, and other epidemiological-biological studies about the usage of multiple vaccines are also welcome to complement this discussion.

The Multiple Vaccination Theorem also indicates what vaccines should not be added to the campaign. This situation was analyzed only for Israel, but by a similar procedure, we can extend our analysis to other countries. The Insufficiency Theorem was implemented in the cases where was observed and in all the cases was possible to find that it is possible to adjust the campaign to get .

At the moment of the parameter estimation, most of the studied countries were experienced a new change in the dynamics of the outbreak, despite the massive vaccination campaign. Moreover, for the analyzed cases, the tendency of the active infected and new cases was declining and in fact, the projected curves indicated the end of the pandemic in these countries. Perhaps the problem of national-based models relies on the assumption that countries are not completely closed and the introduction of new cases is possible, explaining why the referred increment is observed almost simultaneously in the nearest countries.

Further research in this direction must consider the incoming variants of SARS-CoV-2 that can resist vaccination. Perhaps this is the reason for the new increments on the daily cases worldwide. Therefore, long-term predictions for the analyzed countries are difficult to claim. Complex scenarios might arise with the spread of the variants. New research of the reduction on the vaccine efficacy and consistent estimations of the variants are required to develop and fit a compartmental model.

Conclusions

In this paper, we presented a model with multiple types of vaccines adapting some of the structure of the VSIR model of Magpantay [31] for adult vaccination and including a general number of imperfect vaccines. This model was specifically for the case of the ongoing pandemic of COVID-19, which motivates the developed theory.

Rather than focusing on the distribution of vaccines across special compartments or the effects of the honeymoon period described by [30], we studied the design of the vaccination campaign. This special interest follows the current situation of the pandemic: while the first papers of vaccination tried to optimize the distribution by minimizing deaths or other metrics, nowadays is important to decide how can the vaccination campaign be improved by incrementing or reducing vaccination rates. In further research, it would be useful to see the effects of the referred honeymoon period and the potential risk associated with the new variants.

Two relevant results of the Sensitivity Analysis are the Multiple Vaccination Theorem and Insufficiency Theorem. The first one states that the reduction of the reproduction number only holds if one derived vaccination parameter () is less than a threshold called . Although it contradicts common sense, if this condition is not met, the reproduction number would increase, and more intriguing, any increment in the vaccination rate results in an increment in the reproduction number. This is odd since we could think that more vaccines always improves the situation, but theory shows evidence against our beliefs. It was also observed that such believes were considered in the vaccination campaigns of some European countries.

Insufficiency Theorem establishes another counterintuitive proposition: vaccines might not be sufficient to curtail an outbreak, especially if the is high enough. Luckily, in all the cases where this theorem applies, it was found that is possible if the rate of the vaccine with better parameters is increased. We expect that this research could bring some mathematical tools to study the vaccination campaigns as objects of research and improve the public decisions in this matter.

Funding Not applicable.

Availability of data and materials Data was collected from two different sources: Records from reported recovered and cumulative infected were taken from the John Hopkins repository [35] and for vaccination data, we used the repository of Our World in Data [36], including the datasets of vaccines per country and vaccines by the manufacturer.

Code availability All implementations of the simulations and symbolic computations were written in Python and are provided in https://github.com/Pherjev/MultipleVaccination.

CRediT authorship contribution statement

Fernando Javier Aguilar-Canto: Writing – original draft, Formal analysis, Software, Investigation. Ugo Avila-Ponce de León: Conceptualization, Methodology, Writing – review & editing, Supervision. Eric Avila-Vales: Conceptualization, Methodology, Writing – review & editing, Formal analysis, Supervision, Project administration.

Declaration of Competing Interest

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