key: cord-0473314-wb6qdm6t
authors: Mondal, Sayantan; Mukherjee, Saumyak; Bagchi, Biman
title: Attainment of Herd Immunity: Mathematical Modelling of Survival Rate
date: 2020-05-29
journal: nan
DOI: nan
sha: dbc8ce417f4961b02bcffc20d43a4c7674bffe06
doc_id: 473314
cord_uid: wb6qdm6t

We study the influence of the rate of the attainment of herd immunity (HI), in the absence of an approved vaccine, on the vulnerable population. We essentially ask the question: how hard the evolution towards the desired herd immunity could be on the life of the vulnerables? We employ mathematical modelling (chemical network theory) and cellular automata based computer simulations to study the human cost of an epidemic spread and an effective strategy to introduce HI. Implementation of different strategies to counter the spread of the disease requires a certain degree of quantitative understanding of the time dependence of the outcome. In this paper, our main objective is to gather understanding of the dependence of outcome on the rate of progress of HI. We generalize the celebrated SIR model (Susceptible-Infected-Removed) by compartmentalizing the susceptible population into two categories- (i) vulnerables and (ii) resilients, and study dynamical evolution of the disease progression. We achieve such a classification by employing different rates of recovery of vulnerables vis-a-vis resilients. We obtain the relative fatality of these two sub-categories as a function of the percentages of the vulnerable and resilient population, and the complex dependence on the rate of attainment of herd immunity. Our results quantify the adverse effects on the recovery rates of vulnerables in the course of attaining the herd immunity. We find the important result that a slower attainment of the HI is relatively less fatal. However, a slower progress towards HI could be complicated by many intervening factors.

We study the influence of the rate of the attainment of herd immunity (HI), in the absence of an approved vaccine, on the vulnerable population. We essentially ask the question: how hard the evolution towards the desired herd immunity could be on the life of the vulnerables? We employ mathematical modelling (chemical network theory) and cellular automata based computer simulations to study the human cost of an epidemic spread and an effective strategy to introduce HI. Implementation of different strategies to counter the spread of the disease requires a certain degree of quantitative understanding of the time dependence of the outcome. In this paper, our main objective is to gather understanding of the dependence of outcome on the rate of progress of HI. We generalize the celebrated SIR model (Susceptible-Infected-Removed) by compartmentalizing the susceptible population into two categories-(i) vulnerables and (ii) resilients, and study dynamical evolution of the disease progression. We achieve such a classification by employing different rates of recovery of vulnerables vis-a-vis resilients. We obtain the relative fatality of these two sub-categories as a function of the percentages of the vulnerable and resilient population, and the complex dependence on the rate of attainment of herd immunity. Our results quantify the adverse effects on the recovery rates of vulnerables in the course of attaining the herd immunity. We find the important result that a slower attainment of the HI is relatively less fatal. However, a slower progress towards HI could be complicated by many intervening factors.

The present COVID-19 pandemic is a dynamic and volatile process with often unpredictable ups and downs in the infected populations that make it difficult to predict its future course. In the absence of any vaccine or definitive drug in the immediate future [1] the fight against COVID-19 is a hard and long drawn bitter battle, with two strategies being put forward. The first is the widely enforced lockdown, quarantine, and social distancing where the spread of the disease is contained at its inception and only a limited fraction of population is allowed to be infected. [2] This model appears to be successful in South Korea and China, and some other Asian countries. [3] The other model is to allow the virus to have a relatively unconstrained transmission so that a large fraction of the people develops the immunity. [4] This is called the herd immunity (HI) that is favoured by Sweden, and was initially discussed by Germany and England, but largely discarded later. HI can be achieved by two ways-(i) by vaccination, and (ii) by infection. The HI approach is based on the understanding that one can obtain the herd immunity in the society if 60-70% of the population gets immunized. Needless to say this herd immunity is preferable through vaccination as happened in small pox and measles. Implementation of both the models has difficulties. Implementation of lockdown and social distancing requires enormous effort, backed up by resources. On the other hand, the HI model could have adverse consequence on the vulnerable citizens, a subject not adequately discussed. In fact, experiences in Italy and Spain show that the demography can be altered in some regions if HI is given an unconstrained run.

Herd immunity ensures an indirect protection from COVID-19 (or any other infectious disease) when a large fraction of the population becomes immunized to the disease. [5] [6] [7] Herd immunity drastically decreases the probability of the presence of an uninfected individual in the vicinity of a presently infected individual. That is, the infected person is effectively quarantined by the surrounding immunized population. Hence, the chain of propagation breaks. In Fig. 1 we pictorially explain the phenomenon of herd immunity.

FIG. 1. A pictorial representation of the herd immunity phenomenon. In the left we have a region with the susceptible population and one infected person. The total susceptibles are further divided into vulnerables and resilients. The infection propagates in an unconstrained manner and after a certain period the region possesses a large fraction of immunized population (right). After this immunisation any further infection cannot propagate and indirectly protects the susceptibles. In addition to that, multiple infected persons cannot do further harm. The colour codes are maintained throughout this paper.

This can happen by providing the population with a vaccine or by getting cured after being infected. In the case of COVID-19 pandemic, as of now, we are unsure regarding the success of a vaccine and the latter is the only option to attain HI. However, the herd immunity threshold (HIT), that is the minimum fraction of population needs to get immunized in order to eradicate the disease, is different for different infectious diseases. [8, 9] For example, HIT for measles is 92-95% and for that of SARS-1 it is in the range of 50-80%.

Researchers around the world are exploring mainly two aspects of this disease-(i) the microscopic and clinical aspects which would eventually lead to drug discovery and vaccine preparation, [1, 10] (ii) the demographic aspects which lead to policy making and timeline prediction. [2, 3, 11, 12] The latter requires effective mathematical modelling and crowd simulations. However, these models fail to predict the real scenario because of some inherent assumptions and limitations. Although a lot of interesting new studies are emerging in both categories in the context of the recent coronavirus pandemic, the issue of herd immunity and its fatality are not studied.

There are several mathematical models which have been employed in the context of epidemic modelling, for example, the famous Kermack-McKendrick (KM) model which has been used extensively to study the spread of infectious diseases like measles, small pox etc. [13, 14] At the core of this model lies a system of three coupled differential equations for susceptible (S), infected (I) and removed (R) (cured and dead) populations, that is, the famous SIR model (Eq. 1). [15] [16] [17] At the onset of an epidemic S becomes I and I eventually becomes R, but R can never become S or I because of acquired immunity.

Eq. 1 describes the three coupled non-linear differential equations of the KM model where k S→I is the rate of infection and k I→R is the rate of removal (recovery and death). In the conventional SIR model k S→I and k I→R are written as α and β respectively. In principle the rate constants should be time and space dependent, that is, non-local in nature. But it is difficult to predict the functional form of the rate constants with time-it could be periodic, decaying or stochastic in nature. The applicability of this model is for a homogeneous population distribution and mass transmission at a large scale. [13] An important quantity is the basic reproduction number (R 0 ) which is an estimate of the number of secondary infection from one primary infection. [18] The value of R 0 is intimately connected with the herd immunity threshold (H t ) discussed above. [9, 19] (Eq. 2) Hence a correct determination of the basic reproduction parameter, R 0 , is important.

It is clear from Eq. 2 that a higher value of R 0 increases the herd immunity threshold. For SARS-Cov2 the value of R 0 shows a large dispersion and as a consequence we cannot predict the value of H t . For COVID-19 the average value of R 0 is estimated to be in the range of ∼2.0-3.0 but it can possess spatial heterogeneity and time dependence in reality. [20, 21] If one considers R 0 to be in the range of 2.0-3.0 the value of H t would be in between 50%-66%.

In the light of SIR model [Eq.(1)] R 0 can be defined as

Eq. 3 provides a different definition of R 0 and can be understood as follows. If we assume that (the S the fraction of susceptible population) is near 1.0 at the beginning (as there are very few infections compared to a huge population), then R 0 could be equal to unity if the two rate constants are equal. This means that the number of infection and recovery are same at any time. In this situation the disease remains under control. R 0 > 1 causes an epidemic as it challenges the capacity of the healthcare facilities. However, for different region the value of R 0 could be different depending on the intensity of region wise preventive and healthcare measures.

In this work we ask the following questions-(i) what are the relative magnitude of the fatality to the vulnerable and resilient populations if we attempt to achieve HI without a vaccine? (ii) What is the dependence of the fraction of survival on the rate of the attainment of HI? These two issues are widely discussed all over the world. Here we seek answers to these two important questions by employing a modified Susceptible-Infected-Removed (SIR) model and cellular automata (CA) simulations.

The rest of the paper is organised as follows. In section II we describe the mathematical model and the CA simulation protocols. Section III consists of the results from numerical solutions of the modified SIR model and simulations, accompanied by detailed discussions. This section is further divided into several sub sections. In section IV we summarize and conclude our study.

We modify the celebrated SIR (Susceptible-Infected-Removed) model by dividing the entire susceptible population into two parts, namely vulnerable (Vul) and resilient (Res). In the context of the corona virus disease, the vulnerable category consists of persons who are above 60 years of age or have pre-existing medical conditions like diabetes, heart and kidney disease, and lung conditions. [22] The rest of the population is termed as resilient who have a greater chance of getting cured. We achieve such classification by employing different rate constants associated with their recovery. This is based on the available data on the coronavirus disease. The scheme of this classification is described in Fig. 2 .

FIG. 2. Schematic representation of the modified SIR network model. Here the susceptible (S) population is divided into SV and SR that represent elderly and younger people respectively. A part of the fraction SV gets infected and creates IO fraction of infected population. A part of the remaining fraction of the population, that is, SR gets infected and creates IR fraction of the infected population. Both IV and IR get either cured (C) or dead (D). Naturally the rate of recovery for the younger fraction of the population is more than that of the older infected population. On the other hand, the rate of death for the older population is more than that of the younger invectives.

We follow the scheme described in Fig. 2 and formulate a system of eight coupled non-linear differential equations [Eqs. 4 -11] .

In the following, we explain the complex set of equations. Here I(t) is the number of total infectives at any time t, that is I(t) = I V ul (t) + I Res (t). This is the variable that couples the two population sub-categories. k(t) are the rate constants associated with processes that are described in the subscript with an arrow.

We would like to point out that the rates in above equations of motion are all assumed to be time dependent. These rate constants contain all the basic information and also connected with R 0 . In our earlier study, we employed a time dependent rate to produce certain features observed in the time dependence of new cases such a double peaked population structure. [12] The time dependence of rate can be employed to include certain dynamical features like crossover from local contact to community transmission. It is worth stressing that the modelling of these time dependent rate constants plays a pivotal role in the SIR scheme.

We propagate these equations numerically to obtain the respective temporal evolution of each kind of population fraction. From the temporal profiles we can extract several important quantities after a long time (that is, the end of the spread), for example, (i) the peak height of the active infected cases, (ii) the fraction of cured population, (iii) the fraction of dead population, (iv) the fraction of uninfected population, (v) time required to reach the immunity threshold etc. We can regard these equations together to form a system of reacting species, as in a system of chemical reactions.

We solve these equations with two different sets of the rate constant values and aim to understand the relative damages to the vulnerable and resilient population. The values of rate constants are provided in Table I . We keep k S V ul →I V ul and k S Res →I Res the same which depicts the same probability of getting infected for both the subcategories. However, the rate constants associated with recovery and death differs in orders of magnitude between Vul and Res.

We now discuss the procedure we follow to assign different rate constants to the vulnerables and resilients. In a previous study we estimated the values of k S→I and k I→R by fitting the infected/cured/death vs. time data for India (source: www.covid19india.org). [12] We plot the rate of change of the cured (dC/dt) and dead (dD/dt) population against the infected population to find the slope that gives the rate. This procedure provides us with required estimates of k I→C and k I→D . For India, till 27 th May, the estimated values are k I→C = 0.026 day −1 and k I→D = 0.0013day −1 . That is, k I→C is approximately 20 times of k I→D . However, for countries like Italy, Spain, and USA k I→D was significantly higher. This comparison however takes no cognition of the relative time scales, and therefore should be taken with care.

These values are mean field in nature and contain enormous spatial heterogeneity. If we see the state wise (or district wise) statistics we find a large dispersion. On the other hand, the determination of k S→I is not that straight forward as the equations containing k S→I are non-linear in nature in the SIR model. Hence one needs to obtain a good estimate of R 0 and calculate k S→I from Eq. 3. As mentioned above, R 0 also exhibits spatiotemporal heterogeneity which makes the problem of estimating the rate constants even more challenging. For example, in Italy R 0 has been estimated to be ∼3.0-6.0 and in the Hunan province of China it is ∼1.73-5.25. [23] In a recent study on Wuhan, the transmission rate (k S→I ) is assumed to vary from 0.59 to 1.68 day −1 . [24] However, the data required to extract the rate constants associated with the two individual sub-categories, namely, vulnerable and resilient, are not available separately. As the values of the rate constants are connected to the basic reproduction number (R 0 ), we choose the inputs, by preserving the basic features, such that the average value of R 0 yields an acceptable number, in light of acquired information. Next we tune the parameters such that the maximum of the active cases falls in the range of ∼60-90 days, as observed for most countries. We note that we consider these values only to study the trends and do not strictly correspond to any particular region in reality. On the other hand, for set-2 R V ul 0 = 5.20 and R Res 0 = 1.42. We obtain these values by considering each of the population to be individually normalised (that is 100%). In such a situation the effective R 0 can be calculated as follows (Eq. 12).

Here N V ul and N Res represent the number of people in the vulnerable and resilient category respectively. In all our calculations we start with total infected fraction as 0.001 and vary the percentage of vulnerable populations from 5%-40%. By using Eq. 12 we calculate the effective R 0 values for different ratio of vulnerable to resilient population. We find R 0 varies from 1.03 to 1.88 for set-1 and 1.61 to 2.93 for set-2. In a way, set-1 represents a more controlled situation compared to set-2 (Table II) . 

Stochastic cellular automata (CA) simulations give a microscopic and nonlocal picture of the problem at hand. Such simulations are often used to model several physical phenomena. [25] [26] [27] [28] [29] [30] [31] Unlike the mathematical model, CA simulations can directly establish a physical map of the disease-spread. Moreover, we incorporate several region specific and disease specific parameters in our CA simulations that give a general outlook to our investigations. A detailed list of the parameters and associated symbols can be found in our previous work. [12] The spread of COVID-19 is strongly inhomogeneous. So, a homogeneous model fails to capture many aspects. In a real-world scenario, the non-local description may often become important in determining the fate of a pandemic in a given geographical region. In such a case, the population parameters are space-dependent. Moreover, the rate constants also have a spatial distribution. Hence, solutions of these equations are highly non-trivial and a large scale cellular automata simulation may capture these inherent spatiotemporal heterogeneities.

In this work, we neglect the effects of social distancing and quarantine, since we aim at establishing a relation between the percentage of mortality and immunization by an unhindered transmission of the disease within the whole population. Calculation of the rates of transmission and recovery/death can often be difficult due to several reasons like unavailability of data, political or demographic complications etc. This becomes particularly nontrivial when we consider the process with respect to a given population distribution of vulnerable and resilient individuals. The probabilistic approach employed in our simulations makes it easier to study the process, since obtaining an average probability for each of the processes is much more practical.

We use the Moore definition [32] [33] [34] to denote the neighbourhood of a given person. The salient features of our simulation are detailed in our previous work. [12] Here, we summarize our CA simulation methodology.

We start we a land randomly occupied by susceptibles and infectives. The population distribution is such that 5% and 0.05% of the total available land is covered by susceptibles and infectives respectively. We divide the population into vulnerable and resilient individuals with respect to their probabilities of recovery (P V ul R and P Res R ). Vulnerables primarily include people above the age of 60. This also includes people with serious health issues, who are more prone to get deceased if infected. [35] [36] [37] The resilients, on the other hand, are the young fraction of the society with no severe health conditions. When an infective comes in the neighbourhood of a susceptible, the latter is converted to an infective with a given probability of transmission which is considered to be equal and time independent (constant) for both vulnerables and resilients. The time period of infection is determined by probability of recovery and the probability of remaining infected in a given simulation step. In this work, we consider the latter to be 0.99. [12] An individual, once cured from infection, becomes immune to the disease.

We run our simulations for a given number of steps (N ). It should be noted that the time unit is not welldefined for this simulations. To get an estimate of time, the results need to be compared with our theoretical model.

Here we present the results from the numerical solutions of Eqs. 4 -11 in Fig. 3 . We choose two sets of rate constants, set-1 (Fig. 3a) and set-2 (Fig. 3b ) and obtain the changes in the population of vulnerables and resilients. With our choice of parameters (Table I) for set-1 we observe 40.8% increase in the immuned population. In order to achieve the 40.8% immunity a region loses 4.7% of its resilient population and 34.3% vulnerable population. On the other hand, for set-2 a region loses 7.9% of its resilient population and 57.1% of its vulnerable population in order to achieve ∼68% immunity (that could be the HIT for COVID -19) . Hence, it is clear that for both the two cases the vulnerables are significantly affected. We note that with an increased infection rate the timescale of the saturation of the temporal profiles are drastically reduced. The graphs that are presented in Fig. 3 are obtained for 20% initial vulnerable population.

In Fig. 4a , we show the time evolution of the total immunity percentage. In order to study the effect of Table I . Plots show the increase in the total immunity (blue) with the decrease in the populations of vulnerable population (maroon) and resilient population (green) for (a) Set-1 and (b) Set-2. In these two calculations we start with V : R = 1 : 4. In both the two cases the percentage demise in the vulnerable population is significantly higher. fast (early) vs slow (late) achievement of the immunity saturation, we plot the percentage survival of the total population against the time required to attain the immunity threshold (t Im ) for different values of k S→I (Fig.  4b) . We find that the percentage of survival increases linearly with increasing t Im . This indicates that a quick achievement of immunity saturation could lead to fatal consequences. If a society opts for herd immunity, it has to be a slow process. Percentage survival (uninfected and cured population) of the total population against tIm. The two quantities show linear dependence. That is, the percentage survival increases as we take more time to reach immunity saturation. Note that both the X and the Y axes are the outcome of the numerical solution and not provided as inputs. The calculations are done using a fixed Vul:Res=1:4 and the rate constants associated with recovery/death are also kept same as given in Table I. To make the immunity gaining process slow (which leads to relatively less casualty), the rate of infection (k S→I ) needs to be brought down. On the other hand, the rate of removal (recovery and death), k I→R , depends primarily on the disease and partly on the presently available healthcare facilities. k S→I can be controlled by employing effective strategies like lockdown, quarantine, and social distancing. Next we vary the % of initial vulnerable population from 5% to 40% and obtain the % of highest active cases (that is the maxima in the temporal variation of I V (t) or I R (t)), % of cured population and % of death. The range is chosen in order to represent different regions/countries. For example, in India only ∼ 8% of the entire population is above 60 years whereas, in countries like Italy and Germany the number is over 20%.

We obtain Fig. 5a -5c for set-1 and Fig. 5d -5f for set-2. In both the cases the variation of the infected peak maxima with % vulnerable shows nearly linear increase with a higher slope for the vulnerables (Fig. 5a and 5d) . Interestingly the % cured ( Fig. 5b and 5e ) and % dead ( Fig. 5c and 5f) shows a nonlinear dependence on % vulnerable. It clearly shows that the damage is huge to the vulnerable population when the % of vulnerables increases.

We plot (Fig. 6 ) the percentage of deaths for both the subcategories against the herd immunity threshold for a given Vul:Res composition (1:4) . This is to show the increasing damage with respect to Ht. We find that the trend is linear for both the sets of parameters and the relative fatality is substantially higher for the vulnerables.

FIG. 6. Percentage outcome of different herd immunity thresholds (Ht) on the vulberable and resilient population. Plot of percentage deaths against Ht calculated from Eq. 2 for (a) set-1, and (b) set-2. In both the two cases the dependence is linear with substantially more damage to the vulnerable population. The values on the Y axes are individually normalised.

Here, we keep the probability of transmission of disease time-independent and equal for both resilients and vulnerables. We change the initial fraction of the vulnerable section of the total population from 5% to 40%. In Fig. 7 we plot the % of cured individuals (resilients and vulnerables) against % of total immunization when the temporal progression of the population reaches saturation. As discussed earlier, herd immunity is obtained when a major section of the population becomes immune, post infection. However, apart from gaining immunity, this process involves the death of many infected individuals according to their survival probability. The probability of recovery of the resilients is higher than that of the vulnerables. Here, these two probabilities are taken as 0.95 and 0.8 respectively. [36, 38] In Fig. 7 the abscissa is the percentage of the total population that becomes immune after recovering from the infection. The ordinate quantifies the percentage of cured resilients and vulnerables with respect to the total initial population. With increase in the immunity attained in the society, a significant decrease in the percentage of cured vulnerable individuals is observed. This implies that higher the percentage of immunization in the total population, greater is the probability of death of the vulnerable section. Hence herd immunity results in the death of a major fraction of the vulnerable population. This stratum of the society includes mainly the old people (age greater than years) and people with serious health conditions or comorbidity. [22, 39] The geographical regions with demographic distributions having higher fraction of the people of age above ∼60 years are among the worst affected. For example, Italy suffered the loss of many aged people as a result of the COVID-19 pandemic. [40, 41] FIG. 7. Percentage of cured resilient and vulnerables in the population on the course of attaining herd immunity. The percentage of cured individuals is shown as a function of the percentage of total population immunized after getting infected. This is obtained by averaging over 100 CA simulations. Green shows the percentage of death for the resilient fraction of the society and maroon denotes the same for the vulnerable people.

In Fig. 8a , we show the time evolution of the fraction of vulnerables and resilients in the total population for different % of initial number of vulnerables. The fractions are calculated with respect to the total initial population. We see that with increase in the initial % of vulnerables, the number of resilients dying show a slight decrease, whereas the number of dead vulnerables increases significantly. This observation is clarified in Fig. 8b . Here we plot the absolute change in the fraction of resilients and vulnerables as functions of the initial % of vulnerables. Both show linear dependence. The gradient (slope) is negative for resilients and positive for vulnerables. However, we find that the absolute value of the slope for the latter is ∼5 times higher than that of the former. This denotes that countries with higher population of elderly and vulnerable people in the society incur a greater loss in the number of vulnerable individuals. 

Now, we keep the initial population distribution fixed at 20% vulnerable and 80% resilient individuals. We change the probability of recovery of these two categories (P V ul R and P Res R ) with the constraint P V ul R ≤ P Res R . Accordingly, we change these two probabilities from 0.6 to 0.8 and 0.8 to 0.95 respectively. We choose these values according to reported case fatality ratios for the SARS-CoV-2 pandemic. [36, 38] For every pair of P V ul R and P Res R we get a value of percentage of vulnerables and resilients who survive and a fraction of the population that gets immunized. In Figure 9 we plot the survival % of vulnerables and resilients in the two perpendicular axes and represent the % immunized as colour codes according to the colour gradation bar on the right hand side. In this contour representation, red denotes low immunity and blue denotes higher immunity.

The survival % of the vulnerables is lower than that of the resilients. The percentage of immunized population is higher (blue) for maximum survival of the resilients as compared to that of the vulnerables. This means that to attain higher immunity in the population, greater number of old and vulnerable people suffer death as compared to resilients. Hence, attainment of herd immunity comes with the cost of a higher mortality of the vulnerable section of the society.

Any epidemic is a dynamic process where time dependence plays a crucial role in the control of the spread and the damage, that is, the outcome. COVID-19 is a pandemic which is currently under intense scrutiny by all and sundry, and many aspects are yet to be understood. Every move by the government, and the population in general, is of crucial importance. Each pandemic comes with unique characteristics that deserve special treatments, not just medical and clinical but also sociological. In each such epidemic, immunity plays a critical role. Spanish Flu mainly attacked the age group between 20 and 30 years of age. This is the age group with maximum immunity. In the case of COVID-19, again we face the sad reality that certain section of the society is substantially more vulnerable than other sections. The vulnerable section consists of age groups which are above 60-65 years of age, and people with comorbidity. There is yet to further classification, although it is conceivable that as we understand the disease better and more precisely, better perception of danger would emerge.

An epidemic often starts by a process of nucleation which is an important phenomenon often studied in physics and chemistry. The process of nucleation is initiated by a sudden appearance of a group of infected individuals in a region. This may be triggered a laboratory accident, or infection from eating wild animals like bats, pangolin etc. or by arrival of infected tourists and so on. The process may be dependent on the nature of the geography and demography of the country or region. The initial period of the process is often slow. After the initial nucleation, the disease spreads by a diffusion process into the susceptible population. Hence, it is a percolation with a temporal evolution.

In order to address the issue of vulnerability of the population and the outcome with the progression of the epidemic, we carry out a theoretical analysis with the objective to analyze the consequences of aiming for herd immunity without vaccine, or a good drug, in the context of the present COVID-19 pandemic. We develop and solve a modified SIR model numerically and by employing cellular automata simulations. We particularly probed the following question: what is dependence of mortality on the rate of herd immunity?

One of the key results of the present study is the dependence of the percentage survival on the rate of attainment of the immunity threshold. We find that a late attainment of the immunity saturation leads to relatively lesser fatality. We show that approximately 50-60% of the vulnerables might lose their lives in order to attain 70% total immunized population. On the contrary the mortality of the resilient fraction of population is relatively low, may be just about 10%. We find a non-linear trend in the dependence of the cured and dead population on the initial population of the vulnerables. This is because as the number of vulnerables increases, the immunity by infection from a larger fraction of population which cannot protect the vulnerables unless deliberate efforts are made that requires intervention.

While we discuss herd immunity by infection in this work, the other, more sustainable option is herd immunity by vaccination. For example, diseases like small pox, polio etc. have been wiped off the face of earth by vaccination. This is particularly crucial for diseases with high mortality rates. However, for any novel disease, preparation of a vaccine can take years. In case of the present COVID-19 pandemic, for instance, extensive research is going on globally in search of a vaccine. [1] However, no promising result has been obtained in almost five months and researchers believe it may take more than a year to prepare the vaccine.

Epidemic modelling: an introduction

MODSIM 2007 International Congress on Modelling and Simulation. Modelling and Simulation Society of Australia and New Zealand

Differential equations and mathematical biology

Statistical methods in medical research

Infectious disease modelling

Reviews of modern physics

Journal of Physics: Conference Series

Proc. of the Australian Conference on Artificial Life

The Lancet

The Lancet Infectious Diseases

The Lancet Infectious Diseases

The Lancet

We thank Prof. Sarika Bhattacharyya (NCL, Pune) and Prof. Suman Chakrabarty (SNCNCBS, Kolkata) for several fruitful discussions and comments. The authors thank the Department of Science and Technology (DST), India for financial support. BB thanks sir J. C. Bose fellowship for partial financial support. S. Mo. thanks Universities Grants Commission (UGC), India for research fellowship. S. Mu. thanks DST-INSPIRE programme for providing research fellowship.