Mathematical epidemiological model of the propagation and prevention of infectious diseases considering COVID-19 vaccination

The duration of the necessity of current restrictions and preventive countermeasures against the COVID-19 pandemic is of great concern. While several previous epidemiological studies have discussed controlling the course of the epidemic with regard to social distancing, vaccination, care capacities, and future scenarios, we adopt an alternative approach to provide insights into the time scale of epidemic propagation in human populations. We show that the reproduction number modified by preventive countermeasures   implies a threshold   to reach the herd immunity level at time   . While   varies moderately for large   , it increases sharply around , and   is infinite below . The transition region for this increase is minute, demonstrating that prevention of infectious diseases must consider the relatively asymptotic nature of their propagation, which varies unpredictably between steady transmission and explosive outbreaks. These results suggest the continuation of preventive countermeasures to suppress the transmission of COVID-19 for many years; if they were discontinued or reduced such that , the disease would be transmitted throughout the considered community, and the required herd immunity level would be reached within 250 days if However, the implementation of vaccination programs could drastically alter this dynamic.


Introduction
At present, the containment of the severe acute respiratory syndrome coronavirus 2 (SARS-  infection that causes coronavirus disease 2019 (COVID-19) is the most urgent issue worldwide [1].
To date, the statistics concerning the global pandemic remain unprecedentedly grim as the largest global public health crisis in over a century continues. As of July 7, 2021, the US had suffered the largest reported damage; the cumulative number of confirmed cases was 33,451,965 (a ratio to the total population of 10.2%). In contrast, the size of the epidemic outbreak in Japan has been relatively small, with 813,976 confirmed cases (0.64%) [1].
Expectations concerning the course and potential end of the current global pandemic remains unclear, and the question of how long the restrictions on private and economic activities will continue remains as well. Socially, the question arises of whether to give priority to our activities or to the elimination of the pandemic. Four scenarios for the end of local epidemics in specified communities may be considered [2], including (1) the eradication of the epidemics by the successful development of effective vaccines and vaccination programs that completely end the crisis, (2) the attenuation of the infection force of the viruses by interventions including vaccination, leading to remarkable reduction in mortality and morbidity, (3) acceptable successes in the development of vaccines and drugs alongside preventative measures to reduce the severity of the epidemics such that they remain as a relatively small problem, and (4) the acquisition of herd immunity. We may be progressing to scenario (1); however, this is an ideal goal, and it is not yet known whether vaccines will be well-distributed globally or whether the developed vaccines will retain their effectiveness when facing virus mutations.
Regarding mutations, the realization of scenario (2) also involves uncertainties. The realization of scenario (4) implies that many individuals must acquire immunity, and society must be ready to accept many fatalities unless sufficient vaccination programs are implemented. These three scenarios are so uncertain that it may be necessary to proceed with scenario (3). In this case, the most immediate concern regarding the COVID-19 pandemic is the duration of the necessity of continued preventive measures such as social distancing, remote work, shutdown or curtailment of opening hours of public facilities, restrictions on social events, and restrictions on travel.
Regarding this question, the present work investigated the times required to achieve herd immunity with tools for mathematical epidemiology and obtained intriguing results. The herd immunity level ℎ is defined as a level of population immunity such that the spread of the disease will decline and cease even after all preventive measures have been relaxed [3]; stated differently, ℎ provides a threshold on which an epidemic starts to decline under the hypothetical condition of no preventive measures being applied.
According to mathematical epidemiology, the herd immunity level is expressed as ℎ = 1 − 1/ , where denotes the basic reproduction number [3,4,5,6]. This number is an important variable for the epidemiology of infectious diseases [7,8,9], defined as the average number of secondary infections produced when one infected individual is introduced into a host population where all members are susceptible. As an epidemic spreads, the number of susceptible individuals decreases, and the number of infectious individuals increases. When one infectious individual passes pathogens on to individuals, and − individuals among them are already immune, the effective number of infected individuals is . The fraction of immunized individuals is 1 − / . If is more (less) than 1, an epidemic grows (declines). Hence, an epidemic threshold can be defined at = 1. This provides the herd immunity level as ℎ = 1 − 1/ .
Suppose that the ratio of to equals the ratio of the total number of susceptible individuals to the total population : / = / . Then, / = 1/ on the threshold because = 1. According to the conventional compartmental model (see Method), / = −1 + / , where and denote contact and recovery rates, respectively. This equation indicates that takes on a maximum value at = / such that increases when is greater than / and decreases when is less than this value. Therefore, we obtain = / , which is consistent with the calculation of the herd immunity level and leads to = on the threshold. If all members of the host population are susceptible, = , and thus = 1 on the threshold; this is in accordance with the definition of the basic reproduction number .
Apart from the basic reproduction number , the present work considers two other reproduction numbers: the effective reproduction number defined as = ( ), and the modified reproduction number . The former considers the decrease in the number of susceptible individuals , which occurs as an epidemic develops. When an epidemic grows perceptibly and enters into community transmission in the host population, the effective reproduction number is given by analysis of case reports. The latter includes a modification in the contact rate representing the impact of nonpharmaceutical interventions and preventive measures, such as social distancing and isolation.
Therefore, we can write = / using the modified contact rate .
Given that changes in the number of susceptible individuals is expected to be negligible during a short time period, the compartmental model provides for the exponential growth of an epidemic; ( ) = (0) , where = / − represents the growth rate. In the early stage of an epidemic, ≅ , and = ( − 1) . In contrast, in the later stages = ( − 1) = ( − 1) . If the contact rate is modified, is replaced with in the growth rate.
While it may be difficult to incorporate the wide variety of complex, interacting factors that affect actual epidemics, the techniques of mathematical epidemiology provide us with a simplified perspective from which to analyse the propagation of infectious diseases in various populations [3,4,5,6,13,14,15,16,17,18]. To this end, it is important to observe the essential features of epidemics.
The capability of mathematical epidemiology to produce clear and precise projection based on various conjectures is among its key benefits [13][14][15]. A simple perspective has been provided to model the spread of infectious diseases, e.g. age-dependence of vaccination effects on rubella epidemics [3,19], biennial epidemics of measles [6,20,21], or influence of heterogeneity in sexual activity on HIV epidemics [3]. The present work presents the results of an investigation of mathematical models for epidemics focusing on the three reproduction numbers , , and , and especially discusses the time scale for disease propagation, as well as the results of a comparison between the actual data on the COVID-19 epidemics in Tokyo [10,11] and NYC [12] and the computational results derived mathematical models.
Since the beginning of 2020, numerous studies have been published on COVID-19 and SARS-CoV-2 in fields ranging from microscopic to macroscopic, including molecular biology, medical and clinical studies, epidemiology, and political decision making, with each field applying its characteristic methodologies. As examples of reports with similar focus and methodology to this study, research on COVID-19 modelling has considered countries or regions such as France [22], Italy [23], the Basque country [24], Spain [25], Portugal [26], Hong Kong [27], China [28,29,30], and the UK [31]. These studies discussed recurrence and rebound of the epidemics, effects of protection measures and lockdown, widespread testing, contact tracing, care capacities, or scenarios of pandemic management.
In this context, as an alternative approach to model the progression of the COVID-19 pandemic, the present work has evaluated the time required to reach the herd immunity level in a given community, using the epidemiological compartment model [3,4,5,6]. This approach enabled us to discuss the time scale of epidemic propagation through populations and the effects of vaccination. We provide insights into factors supporting the spread of infectious diseases and illustrate the implications of vaccination mathematically. The implementation of vaccination not only increases the number of immunized individuals but also accelerates the achievement of the herd immunity level. This latter effect may be rather important for the elimination of epidemics in general, and specifically for the current global pandemic. The -curve in Fig. 1(a) was calculated with = / = 2.5 (a value in the range plausible for COVID-19 [24,30,32]). When model parameters representing the implementation of preventive measures and changes in the contact rate of individuals were modified, the representative point traced a different trajectory. Figure 1(b) shows trajectories calculated with the modified reproduction number = / due to the reduction in the contact rate from to . When was set to 2.0, the curve could intersect the herd immunity level (green solid curve). When = 1.6, the curve could still pass across the herd immunity level (green broken curve); however, its final value point was very near the herd immunity level. When = 1.5 or less, the curve could not reach the herd immunity level (red curves).

Results
Next, suppose that preventive countermeasures are cancelled in the middle of an epidemic, and the value of the contact rate is restored from to . Then, the curves that cannot reach the herd immunity level under the preventive countermeasures come to reach that level. Figure 1(c) depicts this situation. For example, the blue curve almost converges to the + = 0 line around =0.65, i.e. the epidemic was almost eliminated, and the relaxation of preventative measures caused a further increase followed by an eventual decrease after the herd immunity level was reached. When the timing of countermeasures being rescinded was earlier, the trajectories were closer to the inherent case. However, the peak position of the trajectories remained at = 1/ ; this implies that the herd immunity level is unchanged by the implementation of vaccination.
With the derivation of the trajectories, we have calculated the time over which a trajectory reached the herd immunity level. Figure 2(a) shows as a function of (these calculations do not consider the suspension of preventive countermeasures). We found that increased gradually with decreasing ; was less than 250 days when was greater than 1.6 under the condition of E(0)/N ≥ 10 (see Discussion). However, slightly above = 1.5, this increase became sharp, and diverged for < 1.5. Figure 2(b) shows the effects of vaccination on . First, it is found that the implementation of vaccination tends to transform the curves into a sigmoid form, distorting and flattening them in the region of < 1.5. This effect removes the singularity of at = . Second, the height of the curves, especially in the flattened part, becomes smaller with increasing vaccination rate. When was larger than 0.02, the curve was not even sigmoid, but rather almost linear.
To investigate the actual situation of the COVID-19 epidemics, we estimated the effective reproduction numbers using public data on Tokyo [10,11] and NYC [12]. Figure   In NYC (Fig 3b), at the beginning of March 2020, the reproduction number was as high as 11.5 (95% The frequency distribution of the numbers estimated for the period until January 2021 are shown for Tokyo (Fig. 4(c)) and NYC (Fig.4(d)). This estimate did not consider the period after February 2021 to compare with the mathematical model excluding the mixing effects of reduction in the contact rates and vaccination programs. The mean value was 1.07 and 1.00, and the standard deviation was 0.47 and 0.25 for Tokyo and for NYC, respectively, for the period from April 2020 to January 2021.

Discussion
These results suggest that if the basic reproduction number was modified by preventive Because the proportion of the cumulative total number of individuals infected with COVID-19 was less than 1% and 2% as of January 2021 in Tokyo and NYC, respectively [10,11], the effective reproduction number can be approximated by . Except in the early stages of the epidemics, fluctuated between 0.0 and 3.0 in Tokyo and between 1.0 and 2.0 in NYC (Fig. 3), i.e. has continued to increase and decrease across the value of 1.5, which is the threshold for the divergence of . This fluctuation reflects the dilemma that periods of rapid epidemic propagation and of associated social restrictions have repeatedly alternated since the onset of the pandemic. The distributions of (Fig.4) also illustrate the challenge of determining whether to prioritize the prevention of disease or the continuance of economic activities, demonstrating the circumstances under which both cities have found the COVID-19 pandemic unmanageable. It may be considered likely that many countries have been experiencing similar situations. The abovementioned scenario (3) would oblige us to continue the general living conditions of the current emergency for many years, and it may not be possible to eliminate the pandemic within a few years. This may be a matter of concern in countries where vaccination programs are not progressing.
However, it is probable that the fluctuation of around 1.5 in Fig.3 was not intended but rather a fortuity because individuals generally do not know that ≅ 1.5 is a threshold for . This suggests that the threshold for may be of benefit as a reference value for evaluating the future course of the COVID-19 pandemic. Regarding this point, Fig. 3 illustrates that the propagation speed at the early stage of an epidemic is small, not depending on . Thus, unless we do not know the existence of a threshold for , we cannot find any singularity with asymptomatic increase in the epidemic growth at the early stage, which indicates the threshold. Figure 3 also shows that the spreading velocity remained small until half of , and then a notable increase occurred. Therefore, the trajectories in If the contact rate is modified, is replaced with , shifting the peak position and landing point of the representative point ( Fig.1(b)). This shift of thetrajectory leads to the sharp increase in (Fig.2) with the following mathematical explanation. Two factors affect this sharp increase.
First, when the landing point is shifted to the region of larger than the herd immunity level ( = 1/ ), the representative point cannot reach the herd immunity level and becomes infinitely large.
Second, the compartmental model shows that the rate of change of the number of infectious individuals is proportional to that number itself (see Methods). This is a typical rate equation. Therefore, the rate of change of was small when is small, implying that the motion of the representative point The trajectories displayed in Fig.1(d) have the same value of at the peaks irrespective of the vaccination rates, demonstrating that the implementation of vaccination did not change the herd immunity level [3]. This is because vaccination does not reduce the contact rate (hence, the basic denotes a vaccination rate. It is noted that, in this equation, the vaccination rate is not included in the first term of the right-hand side but rather in the second term. The essential point concerning the effects of vaccination is given below. We found in Fig. 3(b) that the vaccination removes the singularity in at = . Therefore, a reduction in did not prevent the representative point from reaching the herd immunity level; specifically, preventive countermeasures did not prevent the epidemic status from reaching the herd immunity level. Here, it is noteworthy that this achievement of the herd immunity level proceeds without considerable epidemic spreading and expansion, demonstrating the special role of vaccination. This stands in drastic contradistinction to the case without vaccination.
The mathematical interpretation is that the second term in the equation yielding plays a role to update the state of infectious individuals to the state indicating recovered or immunized individuals (see Method). Apart from the numerical calculations shown in the present study, the removal of the singularity by this vaccination term may be proved by calculus for the equations in the compartmental model.
Recent studies showed that heterogeneity of populations reduces the herd immunity level [32]. By introducing age and activity heterogeneities into the population models for SARS-CoV-2, herd immunity was achieved in the models at a population-wide infection rate of ∼ 40%, which is substantially less than the classical herd immunity level of 60% obtained through homogeneous immunization. Although this result is intriguing, such quantitative details of the estimate of the herd immunity level may not affect the qualitative feature obtained in the present work because the sharp increase in the time required to reach the immunization level occurs in the interrelationship and superposition of the fixed point and the landing point of the representative point, as described above; accordingly, the detailed location of the herd immunity level does not affect the nature of the sharp increase in . It is interesting, however, to examine effects of heterogeneities by using more complicated mathematical models.
The experimental results suggest that the classical herd immunity level of 60% would be obtained within a finite period of time if homogeneous immunization was implemented: ~500 days with the rate = 0.001, or ~100 days with = 0.02 (Fig.2b).
To compare the model calculations with the actual course of the COVID-19 pandemic, this study has estimated the effective reproduction numbers using reported cases for Tokyo [10,11]  This study has considered the time required to reach the herd immunity level as an index for evaluating the time scale of epidemic propagation and has provided several insights to illustrate the processes of such propagation. Further studies may yield a clearer understanding regarding this problem by using different indices.
In addition, the results derived herein are applicable to infectious diseases, which may emerge or reemerge in the future. Although each infectious disease possesses a distinctive parameter set for the compartmental model that determines the herd immunity level, the characteristic that increases sharply in a narrow range of arises irrespective of the parameter set. Therefore, precise values of , as a function of , can be calculated only by adjusting the parameter set to the target infectious disease and solving the compartmental model.
The results of the present work suggest that , which is connected with , may be appropriately used as a reference for the ascertainment of epidemic status and the timing of related decision-making.
When vaccination campaigns are conducted, it is important to estimate both and the vaccination rate and to inspect the status of the epidemic carefully. In addition, the results demonstrated that the implementation of vaccination not only increases the number of immunized individuals but accelerates the achievement of the herd immunity level. This latter effect may be of considerable importance for the elimination of epidemic diseases including the COVID-19 pandemic.
The SEIR compartmental model has proven useful for the understanding and prediction of infectious diseases, e.g. this model was able to provide epidemic curves for measles or rubella, which fit to case reports [3,20]. Thereafter, the SEIR model has been applied to other subjects through suitable modifications, e.g. incorporating seasonality into the model gave explanation for chaotic epidemic patterns of measles [6,21,22].

Dataset
The present work used the daily reports of laboratory-confirmed cases.

Back calculation
We used the back-calculation method [34]  The assumption that infections occur at random during an infectious period of an infectious individual leads to the conclusion that the average serial interval is the sum of the average latent period and half the average infectious period. Accordingly, it seems to be relevant to consider that the serial interval is comparative to the period from infection to PCRtesting. For SARS-CoV-2 infections, this distribution is a Gamma distribution with the mean and standard deviations of 7.5 and 3.4 days, respectively [33]. The reported cases were smoothed by 7day-average, and the above equation was solved by deconvolution.

Effective reproduction numbers
We evaluated effective reproduction numbers from the case reports adjusted by the back-calculation and the SEIR model. The SEIR model provides an effective reproduction number Re as = ( + )( + )/ [5]. We set 1/ =5.2 days and 1/ =4.6 days, following previous studies for Wuhan [33].    NYC (d). The mean value was 1.07 and 1.00, and the standard deviation was 0.47 and 0.25 for Tokyo and for NYC, respectively, for the period from April 2020 to January 2021.