Assessing the impact of vaccination in a COVID-19 compartmental model.

BACKGROUND: The aim of this research is to understand the role played by vaccination in the dynamics of a given COVID-19 compartmental model. Most of all, how vaccination modifies the stability, sensitivity, and the reproduction number of a susceptible population.
METHODS: The proposed COVID-19 compartmental model (SVEIRD) has seven compartments. Namely, susceptible (S), vaccinated (V), exposed (E, infected but not yet infectious), symptomatic infectious (I s ), asymptomatic infectious (I a ), recovered (R), and dead by Covid-19 disease (D).We have developed a computational code to mimic the first wave of the coronavirus pandemic in a state like New York (NYS).
FINDINGS: First a stability analysis was carried out. Next, a sensitivity analysis showed that the more relevant parameters are birth rate, transmission coefficient, and vaccine failure. We found an alternative procedure to easily calculate the vaccinated reproductive number of the proposed SVEIRD model. Our graphical results allow to make a comparison between unvaccinated (SEIRD) and vaccinated (SVEIRD) populations. In the peak of the first wave, we estimated 21% (2.5%) and 6% (0.8%) of the unvaccinated (vaccinated) susceptible population was symptomatic and asymptomatic, respectively. At 180 days of the NYS pandemic, the model forecast about 25786 deaths by coronavirus. A vaccine with 95% efficacy could reduce the number of deaths from 25786 to 3784.
CONCLUSION: The proposed compartmental model can be used to mimic different possible scenarios of the pandemic not only in NYS, but in any country or region. Further, for an unvaccinated reproductive number R > 1, we showed that the vaccine's efficacy must be greater than the herd immunity to stop the spread of the COVID-19 disease.

Introduction

Nowadays several vaccines intended to provide immunity against COVID-19 have been already authorized worldwide for emergency use. Clinical trials have shown that a given vaccine's efficacy varies from 60% to 95% [1]. Unfortunately, vaccines' production cannot keep with the world demand, bringing consequently a slow immunization process, i.e., a low coverage rate. Furthermore, in many countries the total number of coronavirus positive cases are increasing due to the emerging of several new coronavirus strains. Because of this fact, coronavirus' vaccines will need to be boosted, and perhaps the human population will have to be vaccinated periodically. These uncertainties suggest that the coronavirus disease will stay with us for many years to come. Therefore, a better understanding of the role played by a vaccine in the coronavirus epidemic is of great health interest. One step in that direction, it is to consider mathematical compartmental models to capture the key features of the coronavirus pandemic. Many authors already have applied different mathematical models to study the coronavirus pandemic [2], [3], [4], [5]] but only few have included vaccination [[6], [7], [8]].

In this paper, we have modified a SEIRD compartmental model by adding an extra vaccinated compartment to make a comparison between vaccinated and unvaccinated populations. In particular, we would like to have an answer for several queries related to the coronavirus dynamics. Namely: i) How much an unvaccinated infected population could be reduced by a vaccine? ii) How the reproduction number changes when the susceptible population is vaccinated? iii) For a given country or region what is the optimal vaccine (in terms of efficacy and coverage rate) to stop the spread of the coronavirus disease?

This paper is organized as follows. In section 2, we write seven first order coupled differential equations to describe mathematically our proposed SVEIRD compartmental model (Fig. 1 ). Disease free and endemic equilibrium points are derived in a closed form in sections 3, 4, respectively. The stability of these equilibrium points is discussed in Sec. 5. In section 6, we use an alternative procedure (other than the Next Generation Matrix [9,10]) to obtain an analytical expression of the vaccinated reproduction number (. Finally, in section 7, we applied the proposed SVEIRD compartmental model to study the coronavirus epidemic in a state like New York. Because there are still many uncertainties in the SVEIRD parameters values and in the understanding of the corona virus disease, we do not claim that our results for NYS should be considered definitive, Although we are aware that NYS had to endure two coronavirus waves, in this work we shall attempt only to simulate one wave of six months of duration. Moreover, our analytical, graphical, and numerical results could be useful as test-bed calculations of more complex and realistic corona-virus models. Here, we shall mention that Schneider et al. [11] have developed a COVID-19 pandemic simulation tool based on an extension of a SEIR model.

Fig. 1

The SVEIRD compartmental mode.

Notice that for Sections 2, 3 all derivations are obtained in a closed form and presented as analytical expressions. This is to benefit any interested reader who could use these results in a country or region besides NYS. Numerical results and related figures for NYS are discussed in Section 4.

Methods

In Fig. 1, the proposed biomathematical (SVEIRD) compartmental model assumes vaccination only to the susceptible population, no reinfection, and neither inherited immunity. The seven compartments are: susceptible (S), vaccinated (V), exposed (E, infected but not yet infectious), symptomatic infectious (I ), asymptomatic infectious (I ), recovered (R), and dead by coronavirus (D).

The meaning of the parameters connecting the different compartments is given bellow in Table 1 .

Table 1

The SVEIRD model's parameters of Fig. 1.

Parameters	Name
Λ	Recruitment
Β	COVID-19 transmission coefficient rate
γ1	Natural mortality rate
γ2	COVID-19 symptomatic mortality rate
γ_3	COVID-19 asymptomatic mortality rate
ω1	Vaccine coverage rate
θ1	Vaccine failure
Ε	Incubation rate
K	Proportion of exposed to the symptomatic population
k_1	Proportion of symptomatic to recovered population
k_2	Proportion of asymptomatic to symptomatic population
k_3	Proportion of asymptomatic to recovered population
k_4	Proportion of asymptomatic to death by COVID-19 population
Φ	Proportion of symptomatic to symptomatic
α1	Recovery rate of the symptomatic population
α2	Recovery rate of the asymptomatic population

Eqns. (1), (2), (3), (4), (5), (6), (7) mathematically describe the dynamics among the different compartments showed in Fig. 1. All compartments’ populations of Fig. 1, have been normalized by

the total population (N = S + V + E + I  + I  + R + D). The new normalized variables will be still denoted by the same letters as in Fig. 1. Namely, and where

Adding Eqns. (1), (2), (3), (4), (5), (6), (7) we obtain This implies when , and therefore solutions of Eqns. (1), (2), (3), (4), (5), (6), (7) must be in the following positively invariant subset of R7: {()/ }

Results

Equilibrium points

To obtain the equilibrium points we set Eqns. (1), (2), (3), (4), (5), (6), (7) to zero. For example, by setting to zero Eqn. (5) we solve it for and inserted into Eqn. (4), to obtain Eqn. (13). From there, it easily follows Eqn. (12). Using a similar algebraic procedure, we obtain Eqns. (14), (15), (16) where

Next, we set Eqn. (3) to zero and substitute Eqns. (12), (13), (14), (15) to obtain the following cubic equationwhere

Disease free equilibrium points

The trivial solution of Eqn. (20); i.e.,  = 0 defines the disease equilibrium point (DFE). Thus, from Eqns. (12), (13), (14), (15), (16) the only non-zero solutions are

Disease endemic equilibrium points

To obtain the disease endemic equilibrium points () we first must solve from.

Eqn. (20) the quadratic equation. Namelywhere and are given by Eqns. (21), (22), (23). Then I is to be substituted in Eqns. (12), (13), (14), (15), (16) to obtain the remaining endemic point values.

Stability

To study the stability of the disease free and endemic equilibrium points we built the Jacobians J and J , respectively. They are explicitly given bellow in Eqns. (27), (28).

As is well known, a local stable equilibrium is achieved when all the eigenvalues of matrices J and J are negative. Otherwise, there is an unstable equilibrium. By setting in Eqns. (27), (28) we shall obtain the respective Jacobians for the unvaccinated population. Here, we shall point out that the eigenvalues of matrices given by Eqns (27), (28) are numerically calculated for NYS in Section 4.

An alternative procedure to calculate the effective reproductive number

The Basic Reproductive Number ( is perhaps one of the most important variables in epidemiology. It is related to the average number of secondary infections and allows to estimate the spread of the disease. In fact, if R > 1 the disease will persist, and it will die out when R < 1. Although the majority of the scientists agree with this interpretation of R, there are some caveats that some authors still have [12,13]. Although the Next Generation Matrix method [9,10,14] is widely used to obtain R, there are other methods to make the same calculation [12]. They are: The Survival Function, the Jacobian, Constant Term of the Characteristic Polynomial, The Graph-Theoretic Method, and the Existence of the Endemic Equilibrium. The pro and cons of using these methods to obtain R, are discussed by Jing et al. [12].

We argue here, that in the case of an infectious disease with no more than three forces of infection the Existence of the Endemic Equilibrium may be the best method to calculate R. This is because, to obtain R all remaining methods face cumbersome calculations. Also, notice that the standard procedure to obtain R is the Next Generation Matrix method [8,9,12]. This procedure does not provide a clear biological interpretation, and occasionally it gives a wrong result [12].

In this section, we shall use a simpler although not trivial calculation within the framework of the Existence of the Endemic Equilibrium method to calculate . To fulfill this purpose, we shall solve Eqn. (20) to obtain the equilibrium point . Notice, >0 means an infectious population threshold which is maintained at a certain baseline level. As is well known, the quadratic equation solutions of Eqn. (20) can be written as showed in Eqn. (29), where , , and  = , are explicitly given by Eqns. (21), (22), (23), respectively. Thus,

Notice from Eqns. (21), (23) that and where

Now we can easily from Eqn. (29) identify the vaccinated reproduction number as

to satisfy the necessary condition Therefore, from Eqns. (29), (30), (31)

Of course, we have verified (Appendix 1) that the value of Eqn. (34) agrees with the one obtained using the standard method i.e., the Next Generation Matrix [9,10,14].

The unvaccinated reproductive number (R) is found by setting in Eqn. (34). Thus

Now, from Eqns. (34), (35) we obtain a relationship between and as follows

Eqn. (36) implies that for any , . This shows analytically that a vaccine diminishes the spread of the coronavirus disease. Also notice, when Eqn. (36) reduces to Eqn. (37)

Therefore, in Eqn. (37), the vaccine's efficacy () will stop the spread of the coronavirus provides the following inequality is satisfied. Namely

Eq. (38) can also be written as follows

Without a vaccine the coronavirus will die out when the population reach herd immunity. However, Eqn. (39) provides an additional interpretation. Namely, to stop the spread of COVID-19 the optimal vaccine's efficacy must be greater than the herd immunity. In other words, the optimal vaccine will move the endemic disease equilibrium to the disease-free equilibrium. A similar conclusion was discussed by Beckley et al. [15].

Thus, provide R > 1, Eqn. (39) allows to choose the optimal vaccine to stop the coronavirus spread in each country or region.

Sensitivity analysis

For the sensitivity analysis we have used the method: Normalized Forward Sensitivity Index of a Variable [[16], [17], [18]]. The vaccinated reproduction number () as given by Eqn. (34) depends on several parameters. To obtain the most relevant parameters of , a sensitivity analysis is carried out. The sensitivity of each parameter xn is given by

From Eqn, (36), the relationship between the unvaccinated ( and vaccinated sensitives between iswhere .

Therefore, for all parameters except and ,  = Also notice that does not depend of and For NYS model's parameters sensitives are given in Table 2 .

Table 2

This Table shows for NYS the positive (, and negative sensitivities (

Parameter	Sensitivity
Λ	1
Β	1
θ1	0.990774
κ_1	0.447102
ε	0.00016164
γ_3	−0.0002613
κ_4	−0.0002613
κ_2	−0.0076717
Φ	−0.0076717
ω_1	−0.0087607
γ_2	−0.0124228
κ	−0.0829919
κ_3	−0.196499
α2	−0.196499
α1	−0.782753
γ_1	−0.991793

The most sensitive parameters for are and θ To diminish (increase) the value of we must diminish (increase) the positive (negative) sensitivities.

Discussion

One of the aims of the proposed model is to understand how vaccination change the dynamics of a SEIRD compartmental models. In Section 1., we have provided analytical expressions to describe the disease free and endemic equilibrium points for SVEIRD and SEIRD compartmental models. These results are needed to study the stability of these solutions and to determine whether the epidemic will persist (i.e., I ≠ 0). If >1 there is only one biological meaningful solution of Eqn. (29). When <1, there are two real negative or two imaginary solutions. This means that there are not endemic points and therefore the disease died out.

We have also developed a novel procedure to derive which is easier to understand and calculate than the standard method, i.e., the Next Generation Matrix [9,10]. We have check that this new alternative procedure to derive analytically , it is valid for any compartmental model with three or less forces of infection.

To illustrate our theoretical procedure, we have chosen the first wave of the NYS corona-virus pandemic. In the first 65 days we expect a highly contagious environment and we have assumed a value of R = 14.6. Thus, using Eqn. (35) we obtained a value of . The next 115 days, social distance, masks, and the lockdown will diminish the value of R and to 2.39 and 0.09, respectively. Next, a MATHEMATICA [19] code is developed to plot the numerical solutions of Eqns. (1), (2), (3), (4), (5), (6), (7). In Table 3 we list the parameters’ values used in the simulations.

Table 3

The SVEIRD parameters values for NYS first wave coronavirus pandemic.

Parameters	Estimated Values	References
Λ	0.3640 10^−4(1/day)	www.osc.state.ny.us [20]
β	0.55/day (first period), 0.09/day (second period)	Estimated
γ1	0.2328 10^−4 (1/day)	https://webbi1.health.ny.gov [22]
γ2	0.0949 (1/day)	[23]
γ_3	0.0697 (1/day)	[23]
ω1	0.0500 (1/day)	Estimated
θ1	0.050	Estimated
ε	0.144	[23]
k	0.749	[23]
k_1	0.001	Estimated
k_2	0.100	Estimated
k_3	0.899	Estimated
Φ	0.200 (1/day)	Estimated
α1	0.0604	[23]
α2	0.0583	[23]

Fig. 2, Fig. 3, Fig. 4, Fig. 5, summarize for NYS first wave the dynamics of the susceptible, vaccinated, infectious and death populations during a six-month period.

Fig. 2

In red and blue, unvaccinated, and vaccinated susceptible populations, respectively. In green, the vaccinated population. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Fig. 3

The unvaccinated symptomatic (red) population. Vaccination reduces the symptomatic population (blue). (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Fig. 4

The unvaccinated asymptomatic (red) population. Vaccination reduces the asymptomatic population (blue). (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Fig. 5

The accumulated unvaccinated (Red) COVID-19 deaths. Vaccination reduces the number of COVID-19 deaths (blue). (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

We show in Fig. 2, the evolution on time of the vaccinated (blue) and unvaccinated (red) susceptible populations. After 180 days about 80% of the initial susceptible population will be vaccinated. As can been seen in Fig. 3, Fig. 4, vaccination will diminish the symptomatic (from 0.21 to 0.025) and asymptomatic (from 0.06 to 0.008) population's proportions) in the peak of the first wave. Accumulated covid-deaths are plotted in Fig. (4). The proposed model forecast for the first six months of the NYS pandemic 25785 deaths due to the coronavirus. Vaccination could have reduced the number of deaths from 25785 to 3784. The total NYS coronavirus deaths reported in the web, after six-months is 33139 [24].

For NYS, vaccination reduces the unvaccinated reproductive number from 14.6 to 0.76 (first 65 days), and from 2.39 to 0.12 for the remaining period (65 days–180 days). It is showed that for R = 2.39, only exist one real positive endemic equilibrium point ( Thus, without vaccination the Covid-19 disease will persist and converge to the endemic equilibrium point (. It is expected that all meaningful solutions of Eqns. (1), (2), (3), (4), (5), (6), (7) will converge to this unique equilibrium point and the disease will become endemic. Regarding the unvaccinated disease-free solutions all of them will converge to the value of Using the respective Jacobians we found that both of these equilibrium points are unstable. For  = 0.12 all solutions of Eqns. (1), (2), (3), (4), (5), (6), (7) will converge to the disease-free equilibrium points (and the Covid-19 disease will be died out.

In Eqn. (38), we showed that the optimal vaccine (1- ) must be greater than the herd immunity to stop the coronavirus spread in a given country or region. This information may be useful in choosing what Covid-19 vaccine's must be acquired. Notice that herd immunity values depend on R, and thus, on the mathematical models and parameters obtained for a specific country. For example, in NYS, at the beginning of the pandemic our model estimated an unvaccinated reproduction number of 14.60. Therefore, our SVEIRD proposed model advocates that the optimal vaccine for NYS should have at least 93% efficiency.

Fig. 6 shows in red (without vaccination) that must be greater than 0.05/day to become endemic. On the other hand, if the population was vaccinated should be less than 1.1/day to died out. Therefore, after the second period, a second wave was expected.

Fig. 6

The unvaccinated (red) and vaccinated (blue) reproductive numbers in function of the transmission coefficient rate (In green.. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

In Fig. 7, Fig. 8 , we plot three-dimensional figures of the intersecting a plane of value 1. These figures are useful to determine what are the values of the transmission coefficient, vaccine failure, or recruitment rate to stop the coronavirus disease.

Fig. 7

In blue, the vaccinated reproduction number (R). In green, a plane of value 1. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Fig. 8

In red, the unvaccinated reproduction number (R). In green, a plane of value 1. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Limitations

There are several limitations in this paper. First, as any compartmental model, it cannot make a perfect forecast. The proposed SVEIQRD compartmental model only provides an approximation that is accurate enough for a better understanding of the pandemic's first wave. Although there is NYS coronavirus data published in different websites [[24], [25], [26]] we cannot make it a direct comparison with our respective results (except for COVID-deaths). This is because the number of symptomatic and asymptomatic patients reported in the literature is based on the number of tests taken daily. There is some reliable data available: the number of COVID-19 cumulative deaths, and the average number for the reproductive number. The compartments' parameters describing the biology of the disease are perhaps the more reliable because they are based in COVID-19 clinical data. In summary, there is not yet high-quality data for the SVEIRD compartmental model parameters. Regarding the effective reproductive number, its value depends on the compartmental model chosen.

Conclusions

In this paper, we focus on how vaccination modifies the dynamics of a coronavirus SEIRD compartmental model. To this purpose a simple but not trivial SVEIRD compartmental model was built. In the first part, we derived in a closed form analytical expressions for the endemic and disease-free equilibrium points and its respective Jacobians. Next, we derive an alternative procedure to obtain the reproductive number. This approach is easier to calculate and understand than the standard method: The Next Generation Matrix. Moreover, we show that the optimal vaccine (in terms of efficacy and coverage rate) for a region or country must be greater than the herd immunity to stop the spread of the Covid-19.

As an application, we developed a MATHEMATICA code to mimic the first wave of the coronavirus pandemic in a state like New York. Two periods (65 days, 65 days-180 days) were considered. We assumed for the first two months a highly infectious (R = 14.6) period, which diminish to R = 2.39 in the subsequent four months.

As expected, Fig. 1, Fig. 2, Fig. 3, Fig. 4 shows that the vaccine reduces the prevalence of the disease. At the peak of the first wave (Fig. 3, Fig. 4), the program forecasts that a quarter of the susceptible population will be infectious (symptomatic and asymptomatic). Vaccination (95% efficacy, 1% daily coverage rate) will reduce the number of infectious patients for an approximate factor of 10. Accumulated COVID-19 deaths (25785) after six months of the pandemic have a 22% difference with the reported values (33139) [24]. Vaccination (95% efficacy, 1% daily coverage rate) will reduce the number of deaths from 25785 to 3784.

Finally, to contain the pandemic in NYS, the optimal vaccines should have at least 93% and 60% efficacies for the first and second period, respectively.

Ethics statement

The data and parameters values were collected via public domain, and thus, neither ethical approval nor individual consent is applicable.

Code availability

(MATHEMATICA code used in this research can be requested to Ernesto.Esteban@upr.edu)

Declaration of competing interest

The authors declare that they have no competing financial interest that has influence in the work reported in this paper.