Attainment of Herd Immunity: Mathematical Modelling of Survival Rate
AAttainment of Herd Immunity: Mathematical Modelling of Survival Rate
Sayantan Mondal, Saumyak Mukherjee, and Biman Bagchi Solid State and Structural Chemisry Unit, Indian Institute of Science, Bengaluru-560012, India
We study the influence of the rate of the attainment of herd immunity (HI), in the absence ofan approved vaccine, on the vulnerable population. We essentially ask the question: how hard theevolution towards the desired herd immunity could be on the life of the vulnerables ? We employmathematical modelling (chemical network theory) and cellular automata based computer simu-lations 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 degreeof quantitative understanding of the time dependence of the outcome. In this paper, our mainobjective 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 thesusceptible population into two categories- (i) vulnerables and (ii) resilients, and study dynamicalevolution of the disease progression.
We achieve such a classification by employing different rates ofrecovery of vulnerables vis-a-vis resilients . We obtain the relative fatality of these two sub-categoriesas a function of the percentages of the vulnerable and resilient population, and the complex depen-dence on the rate of attainment of herd immunity. Our results quantify the adverse effects on therecovery rates of vulnerables in the course of attaining the herd immunity.
We find the importantresult that a slower attainment of the HI is relatively less fatal . However, a slower progress towardsHI could be complicated by many intervening factors.
I. INTRODUCNTION
The present COVID-19 pandemic is a dynamic andvolatile process with often unpredictable ups and downsin the infected populations that make it difficult to pre-dict its future course. In the absence of any vaccineor definitive drug in the immediate future [1] the fightagainst COVID-19 is a hard and long drawn bitter bat-tle, with two strategies being put forward. The first isthe widely enforced lockdown, quarantine, and social dis-tancing where the spread of the disease is contained atits inception and only a limited fraction of populationis allowed to be infected.[2] This model appears to besuccessful in South Korea and China, and some otherAsian countries.[3] The other model is to allow the virusto have a relatively unconstrained transmission so thata large fraction of the people develops the immunity.[4]This is called the herd immunity (HI) that is favouredby Sweden, and was initially discussed by Germany andEngland, but largely discarded later. HI can be achievedby two ways- (i) by vaccination, and (ii) by infection.The HI approach is based on the understanding that onecan obtain the herd immunity in the society if 60-70% ofthe population gets immunized. Needless to say this herdimmunity is preferable through vaccination as happenedin small pox and measles. Implementation of both themodels has difficulties. Implementation of lockdown andsocial distancing requires enormous effort, backed up byresources. On the other hand, the HI model could haveadverse consequence on the vulnerable citizens, a subjectnot adequately discussed. In fact, experiences in Italyand Spain show that the demography can be altered insome regions if HI is given an unconstrained run.Herd immunity ensures an indirect protection fromCOVID-19 (or any other infectious disease) when a largefraction of the population becomes immunized to the disease.[5–7] Herd immunity drastically decreases theprobability of the presence of an uninfected individualin the vicinity of a presently infected individual. Thatis, the infected person is effectively quarantined by thesurrounding immunized population. Hence, the chain ofpropagation breaks. In Fig. 1 we pictorially explain thephenomenon of herd immunity.
FIG. 1. A pictorial representation of the herd immunity phe-nomenon. In the left we have a region with the susceptiblepopulation and one infected person. The total susceptiblesare further divided into vulnerables and resilients. The in-fection propagates in an unconstrained manner and after acertain period the region possesses a large fraction of immu-nized population (right). After this immunisation any furtherinfection cannot propagate and indirectly protects the suscep-tibles. In addition to that, multiple infected persons cannotdo further harm. The colour codes are maintained throughoutthis paper.
This can happen by providing the population with avaccine or by getting cured after being infected. In thecase of COVID-19 pandemic, as of now, we are unsure a r X i v : . [ q - b i o . P E ] M a y regarding the success of a vaccine and the latter is theonly option to attain HI. However, the herd immunitythreshold (HIT), that is the minimum fraction of popu-lation needs to get immunized in order to eradicate thedisease, is different for different infectious diseases.[8, 9]For example, HIT for measles is 92-95% and for that ofSARS-1 it is in the range of 50-80%.Researchers around the world are exploring mainlytwo aspects of this disease- (i) the microscopic and clin-ical aspects which would eventually lead to drug dis-covery and vaccine preparation,[1, 10] (ii) the demo-graphic aspects which lead to policy making and time-line prediction.[2, 3, 11, 12] The latter requires effectivemathematical modelling and crowd simulations. How-ever, these models fail to predict the real scenario be-cause of some inherent assumptions and limitations. Al-though a lot of interesting new studies are emerging inboth categories in the context of the recent coronaviruspandemic, the issue of herd immunity and its fatality arenot studied.There are several mathematical models which havebeen employed in the context of epidemic modelling, forexample, the famous Kermack-McKendrick (KM) modelwhich has been used extensively to study the spread ofinfectious diseases like measles, small pox etc.[13, 14] Atthe core of this model lies a system of three coupled dif-ferential equations for susceptible (S), infected (I) andremoved (R) (cured and dead) populations, that is, thefamous SIR model (Eq. 1).[15–17] At the onset of anepidemic S becomes I and I eventually becomes R, but Rcan never become S or I because of acquired immunity. dSdt = − k S → I SIdIdt = k S → I SI − k I → R IdRdt = k I → R I (1)Eq. 1 describes the three coupled non-linear differen-tial equations of the KM model where k S → I is the rate ofinfection and k I → R is the rate of removal (recovery anddeath). In the conventional SIR model k S → I and k I → R are written as α and β respectively. In principle the rateconstants should be time and space dependent, that is,non-local in nature. But it is difficult to predict the func-tional form of the rate constants with time- it could beperiodic, decaying or stochastic in nature. The appli-cability of this model is for a homogeneous populationdistribution and mass transmission at a large scale.[13]An important quantity is the basic reproduction num-ber ( R ) which is an estimate of the number of secondaryinfection from one primary infection.[18] The value of R is intimately connected with the herd immunity thresh-old ( H t ) discussed above.[9, 19] (Eq. 2) Hence a correctdetermination of the basic reproduction parameter, R ,is important. H t = (cid:18) − R (cid:19) × R in-creases the herd immunity threshold. For SARS-Cov2the value of R shows a large dispersion and as a conse-quence we cannot predict the value of H t . For COVID-19the average value of R is estimated to be in the rangeof ∼ R tobe in the range of 2.0-3.0 the value of H t would be inbetween 50%-66%.In the light of SIR model [Eq.(1)] R can be defined as R = k S → I k I → R S (3)Eq. 3 provides a different definition of R and canbe understood as follows. If we assume that (the S thefraction of susceptible population) is near 1.0 at the be-ginning (as there are very few infections compared to ahuge population), then R could be equal to unity if thetwo rate constants are equal. This means that the num-ber of infection and recovery are same at any time. Inthis situation the disease remains under control. R > R could be different depending on the intensityof region wise preventive and healthcare measures.In this work we ask the following questions- (i) whatare the relative magnitude of the fatality to the vulnera-ble and resilient populations if we attempt to achieve HIwithout a vaccine? (ii) What is the dependence of thefraction 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 questionsby employing a modified Susceptible-Infected-Removed(SIR) model and cellular automata (CA) simulations.The rest of the paper is organised as follows. In sec-tion II we describe the mathematical model and the CAsimulation protocols. Section III consists of the resultsfrom numerical solutions of the modified SIR model andsimulations, accompanied by detailed discussions. Thissection is further divided into several sub sections. Insection IV we summarize and conclude our study. II. THEORETICAL FORMALISMA. Mathematical Modelling
We modify the celebrated SIR (Susceptible-Infected-Removed) model by dividing the entire susceptible pop-ulation into two parts, namely vulnerable (Vul) and re-silient (Res). In the context of the corona virus dis-ease, the vulnerable category consists of persons whoare above 60 years of age or have pre-existing medicalconditions like diabetes, heart and kidney disease, andlung conditions.[22] The rest of the population is termedas resilient who have a greater chance of getting cured.We achieve such classification by employing different rateconstants associated with their recovery. This is basedon the available data on the coronavirus disease. Thescheme of this classification is described in Fig. 2.
FIG. 2. Schematic representation of the modified SIR net-work model. Here the susceptible (S) population is dividedinto S V and S R that represent elderly and younger people re-spectively. A part of the fraction S V gets infected and creates I O fraction of infected population. A part of the remainingfraction of the population, that is, S R gets infected and cre-ates I R fraction of the infected population. Both I V and I R get either cured (C) or dead (D). Naturally the rate of recov-ery for the younger fraction of the population is more thanthat of the older infected population. On the other hand, therate of death for the older population is more than that of theyounger invectives. We follow the scheme described in Fig. 2 and formulatea system of eight coupled non-linear differential equations[Eqs. 4 - 11]. dS V ul ( t ) dt = − k S V ul → I V ul ( t ) S V ul ( t ) I ( t ) (4) dS Res ( t ) dt = − k S Res → I Res ( t ) S Res ( t ) I ( t ) (5) dI V ul ( t ) dt = k S V ul → I V ul ( t ) S V ul ( t ) I ( t ) − ( k I V ul → C V ul ( t ) + k I V ul → D V ul ( t )) I V ul ( t ) (6) dI Res ( t ) dt = k S Res → I Res ( t ) S Res ( t ) I ( t ) − ( k I Res → C Res ( t ) + k I Res → D Res ( t )) I Res ( t ) (7) dC V ul ( t ) dt = k I V ul → C V ul ( t ) I V ul ( t ) (8) dC Res ( t ) dt = k I Res → C Res ( t ) I Res ( t ) (9) dD V ul ( t ) dt = k I V ul → D V ul ( t ) I V ul ( t ) (10) dD Res ( t ) dt = k I Res → D Res ( t ) I Res ( t ) (11)In the following, we explain the complex set of equa-tions. Here I(t) is the number of total infectives at anytime t, that is I ( t ) = I V ul ( t ) + I Res ( t ). This is the vari-able that couples the two population sub-categories. k ( t )are the rate constants associated with processes that aredescribed in the subscript with an arrow.We would like to point out that the rates in aboveequations of motion are all assumed to be time depen-dent. These rate constants contain all the basic informa-tion and also connected with R . In our earlier study,we employed a time dependent rate to produce certainfeatures observed in the time dependence of new casessuch a double peaked population structure.[12] The timedependence of rate can be employed to include certaindynamical features like crossover from local contact tocommunity transmission. It is worth stressing that themodelling of these time dependent rate constants plays apivotal role in the SIR scheme.We propagate these equations numerically to obtainthe respective temporal evolution of each kind of popula-tion fraction. From the temporal profiles we can extractseveral important quantities after a long time (that is,the end of the spread), for example, (i) the peak heightof the active infected cases, (ii) the fraction of cured pop-ulation, (iii) the fraction of dead population, (iv) the frac-tion of uninfected population, (v) time required to reachthe immunity threshold etc. We can regard these equa-tions together to form a system of reacting species, as ina system of chemical reactions.We solve these equations with two different sets of therate constant values and aim to understand the relativedamages to the vulnerable and resilient population. Thevalues 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 thesame probability of getting infected for both the sub-categories. However, the rate constants associated withrecovery and death differs in orders of magnitude betweenVul and Res.We now discuss the procedure we follow to assign dif-ferent rate constants to the vulnerables and resilients. Ina previous study we estimated the values of k S → I and k I → R dC/dt ) and dead ( dD/dt )population against the infected population to find theslope that gives the rate. This procedure provides us withrequired estimates of k I → C and k I → D . For India, till 27 th May, the estimated values are k I → C = 0 . day − and k I → D = 0 . day − . That is, k I → C is approximately 20times of k I → D . However, for countries like Italy, Spain,and USA k I → D was significantly higher. This compari-son 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 enor-mous spatial heterogeneity. If we see the state wise (ordistrict wise) statistics we find a large dispersion. Onthe other hand, the determination of k S → I is not thatstraight forward as the equations containing k S → I arenon-linear in nature in the SIR model. Hence one needsto obtain a good estimate of R and calculate k S → I fromEq. 3. As mentioned above, R also exhibits spatiotem-poral heterogeneity which makes the problem of estimat-ing the rate constants even more challenging. For exam-ple, in Italy R has been estimated to be ∼ ∼ k S → I ) isassumed to vary from 0.59 to 1 . day − .[24]However, the data required to extract the rate con-stants associated with the two individual sub-categories,namely, vulnerable and resilient, are not available sepa-rately. As the values of the rate constants are connectedto the basic reproduction number ( R ), we choose theinputs, by preserving the basic features, such that theaverage value of R yields an acceptable number, in lightof acquired information. Next we tune the parameterssuch that the maximum of the active cases falls in therange of ∼ TABLE I. The values of rate constants used to solve the sys-tem of coupled differential equations [Eq.4 - 11]. The unit ofthe rate constants is day − .Rate Const. Set-1 Set-2 k S V ul → I V ul k I V ul → C V ul k I V ul → D V ul k S Res → I Res k I Res → C Res k I Res → D Res
We invoke two different values of R for the two dif-ferent sub-categories. For set-1 R V ul = 3 .
33 and R Res =3 .
33. The larger value of R for vulnerables arise fromslower rate of recovery, Eq. 3.On the other hand, for set-2 R V ul = 5 .
20 and R Res =1 .
42. We obtain these values by considering each of thepopulation to be individually normalised (that is 100%).In such a situation the effective R can be calculated asfollows (Eq. 12). R eff = R V ul N V ul + R Res N Res N V ul + N Res (12)Here N V ul and N Res represent the number of people inthe vulnerable and resilient category respectively. In all our calculations we start with total infected fraction as0.001 and vary the percentage of vulnerable populationsfrom 5%-40%. By using Eq. 12 we calculate the effec-tive R values for different ratio of vulnerable to resilientpopulation. We find R varies from 1.03 to 1.88 for set-1and 1.61 to 2.93 for set-2. In a way, set-1 represents amore controlled situation compared to set-2 (Table II). TABLE II. The basic reproduction number ( R ) for the pa-rameters described in Table I (set-1 and set-2) for variousratios of vulnerable to resilient population.% % R vulnerble resilient Set-1 Set-25 95 1.031 1.60910 90 1.152 1.79815 85 1.273 1.98720 80 1.394 2.17625 75 1.515 2.36530 70 1.636 2.55435 65 1.757 2.74340 60 1.878 2.932 B. Stochastic Cellular Automata Simulation
Stochastic cellular automata (CA) simulations give amicroscopic and nonlocal picture of the problem at hand.Such simulations are often used to model several physicalphenomena.[25–31] Unlike the mathematical model, CAsimulations can directly establish a physical map of thedisease-spread. Moreover, we incorporate several regionspecific and disease specific parameters in our CA simu-lations that give a general outlook to our investigations.A detailed list of the parameters and associated symbolscan 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 of-ten become important in determining the fate of a pan-demic in a given geographical region. In such a case,the population parameters are space-dependent. More-over, the rate constants also have a spatial distribution.Hence, solutions of these equations are highly non-trivialand a large scale cellular automata simulation may cap-ture these inherent spatiotemporal heterogeneities.In this work, we neglect the effects of social distancingand quarantine, since we aim at establishing a relationbetween the percentage of mortality and immunizationby an unhindered transmission of the disease within thewhole population. Calculation of the rates of transmis-sion and recovery/death can often be difficult due to sev-eral reasons like unavailability of data, political or de-mographic complications etc. This becomes particularlynontrivial when we consider the process with respect to agiven population distribution of vulnerable and resilientindividuals. The probabilistic approach employed in oursimulations makes it easier to study the process, since ob-taining an average probability for each of the processesis much more practical.We use the Moore definition [32–34] to denote theneighbourhood of a given person. The salient featuresof our simulation are detailed in our previous work.[12]Here, we summarize our CA simulation methodology.We start we a land randomly occupied by susceptiblesand infectives. The population distribution is such that5% and 0.05% of the total available land is covered by sus-ceptibles and infectives respectively. We divide the pop-ulation into vulnerable and resilient individuals with re-spect to their probabilities of recovery ( P V ulR and P ResR ).Vulnerables primarily include people above the age of 60.This also includes people with serious health issues, whoare more prone to get deceased if infected.[35–37] The re-silients, on the other hand, are the young fraction of thesociety with no severe health conditions. When an in-fective comes in the neighbourhood of a susceptible, thelatter is converted to an infective with a given probabil-ity of transmission which is considered to be equal andtime independent (constant) for both vulnerables and re-silients. The time period of infection is determined byprobability of recovery and the probability of remaininginfected in a given simulation step. In this work, we con-sider the latter to be 0.99.[12] An individual, once curedfrom 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 well-defined for this simulations. To get an estimate of time,the results need to be compared with our theoreticalmodel. III. RESULTS AND DISCUSSIONA. Numerical solutions of the SIR model
Here we present the results from the numerical solu-tions of Eqs. 4 - 11 in Fig. 3. We choose two sets ofrate constants, set-1 (Fig. 3a) and set-2 (Fig. 3b) andobtain the changes in the population of vulnerables andresilients. With our choice of parameters (Table I) forset-1 we observe 40.8% increase in the immuned popula-tion. In order to achieve the 40.8% immunity a regionloses 4.7% of its resilient population and 34.3% vulnera-ble population. On the other hand, for set-2 a region loses7.9% of its resilient population and 57.1% of its vulnera-ble population in order to achieve ∼
68% immunity (thatcould be the HIT for COVID-19). Hence, it is clear thatfor both the two cases the vulnerables are significantlyaffected. We note that with an increased infection ratethe timescale of the saturation of the temporal profilesare drastically reduced. The graphs that are presented inFig. 3 are obtained for 20% initial vulnerable population.In Fig. 4a, we show the time evolution of the totalimmunity percentage. In order to study the effect of
FIG. 3. Population disease progression as obtained from thesolution of the system of eight coupled non-linear differentialequations presented in Eqs. 4 - 11 as function of time fortwo different situations described in Table I. Plots show theincrease in the total immunity (blue) with the decrease in thepopulations of vulnerable population (maroon) and resilientpopulation (green) for (a) Set-1 and (b) Set-2. In these twocalculations we start with V : R = 1 : 4. In both the twocases the percentage demise in the vulnerable population issignificantly higher. fast (early) vs slow (late) achievement of the immunitysaturation, we plot the percentage survival of the totalpopulation against the time required to attain the im-munity threshold ( t Im ) for different values of k S → I (Fig.4b). We find that the percentage of survival increaseslinearly with increasing t Im . This indicates that a quickachievement of immunity saturation could lead to fatalconsequences. If a society opts for herd immunity, it hasto be a slow process . FIG. 4. The effect of different rates of attaining herd immu-nity on the total population. (a) Plot of the time evolutionof the percentage of total immunized population for differ-ent values of susceptible to infected rate-constants. With in-creasing k S → I we see an increase in the percentage immunityand decrease the time required to reach saturation ( t Im ). (b)Percentage survival (uninfected and cured population) of thetotal population against t Im . The two quantities show lineardependence. That is, the percentage survival increases as wetake more time to reach immunity saturation. Note that boththe X and the Y axes are the outcome of the numerical so-lution and not provided as inputs. The calculations are doneusing a fixed Vul:Res=1:4 and the rate constants associatedwith recovery/death are also kept same as given in Table I. To make the immunity gaining process slow (whichleads 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 , dependsprimarily on the disease and partly on the presently avail-able healthcare facilities. k S → I can be controlled by em-ploying effective strategies like lockdown, quarantine, andsocial distancing. FIG. 5. The effect of the change in the initial percentage ofthe vulnerable population on the relative infection and recov-ery for the sub-categories, namely, vulnerables and resilients.Plots show the dependence of infection peak, percentage curedand dead population for vulnerable (maroon) and resilient(green) population with the initial fraction of vulnerable pop-ulation as obtained from the solution of the modified SIRmodel described in Eqs. 4 - 11. For figures (a)-(c) set-1 andfor (d)-(f) set-2. The quantities show a non-linear dependenceand enhanced fatality for the vurnerables.
Next we vary the % of initial vulnerable populationfrom 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 rangeis chosen in order to represent different regions/countries.For example, in India only ∼
8% of the entire populationis above 60 years whereas, in countries like Italy andGermany the number is over 20%.We obtain Fig. 5a - 5c for set-1 and Fig. 5d - 5f forset-2. In both the cases the variation of the infected peakmaxima with % vulnerable shows nearly linear increasewith 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 hugeto the vulnerable population when the % of vulnerablesincreases.We plot (Fig. 6) the percentage of deaths for both thesubcategories against the herd immunity threshold for agiven Vul:Res composition (1:4). This is to show the in-creasing damage with respect to Ht. We find that thetrend is linear for both the sets of parameters and therelative fatality is substantially higher for the vulnera-bles.
FIG. 6. Percentage outcome of different herd immunitythresholds ( H t ) on the vulberable and resilient population.Plot of percentage deaths against H t calculated from Eq. 2for (a) set-1, and (b) set-2. In both the two cases the depen-dence is linear with substantially more damage to the vulner-able population. The values on the Y axes are individuallynormalised. B. (Stochastic Cellular Automata Simulations
1. Dependence on the initial population distribution
Here, we keep the probability of transmission of diseasetime-independent and equal for both resilients and vul-nerables. We change the initial fraction of the vulnerablesection of the total population from 5% to 40%. In Fig. 7we plot the % of cured individuals (resilients and vulner-ables) against % of total immunization when the tem-poral progression of the population reaches saturation.As discussed earlier, herd immunity is obtained when amajor section of the population becomes immune, postinfection. However, apart from gaining immunity, thisprocess involves the death of many infected individualsaccording to their survival probability. The probabilityof recovery of the resilients is higher than that of the vul-nerables. Here, these two probabilities are taken as 0.95and 0.8 respectively.[36, 38]In Fig. 7 the abscissa is the percentage of the totalpopulation that becomes immune after recovering fromthe infection. The ordinate quantifies the percentage ofcured resilients and vulnerables with respect to the to-tal initial population.
With increase in the immunityattained in the society, a significant decrease in the per-centage of cured vulnerable individuals is observed . Thisimplies that higher the percentage of immunization inthe total population, greater is the probability of deathof the vulnerable section. Hence herd immunity resultsin the death of a major fraction of the vulnerable popu-lation. This stratum of the society includes mainly theold people (age greater than years) and people with se-rious health conditions or comorbidity.[22, 39] The geo-graphical regions with demographic distributions havinghigher fraction of the people of age above ∼
60 years areamong the worst affected. For example, Italy sufferedthe loss of many aged people as a result of the COVID-19 pandemic.[40, 41]
FIG. 7. Percentage of cured resilient and vulnerables in thepopulation on the course of attaining herd immunity. Thepercentage of cured individuals is shown as a function of thepercentage of total population immunized after getting in-fected. This is obtained by averaging over 100 CA simula-tions. Green shows the percentage of death for the resilientfraction of the society and maroon denotes the same for thevulnerable people.
In Fig. 8a , we show the time evolution of the fractionof vulnerables and resilients in the total population fordifferent % of initial number of vulnerables. The fractionsare 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 signif-icantly. This observation is clarified in Fig. 8b. Here weplot the absolute change in the fraction of resilients andvulnerables as functions of the initial % of vulnerables.Both show linear dependence. The gradient (slope) isnegative for resilients and positive for vulnerables. How-ever, we find that the absolute value of the slope for thelatter is ∼ FIG. 8. (a) Population dynamics represented as the temporalevolution of the fraction of resilient and vulnerable sections ofthe population are shown with varying initial distribution ofresilients and vulnerables. The colour bar on the right handside shows the initial % of vulnerables in the total population.(b) The absolute decrease in the resilient (green) and vulner-able (maroon) fractions of the total population as functionsof the initial percentage of vulnerables.
2. Dependence on the probability of recovery
Now, we keep the initial population distribution fixedat 20% vulnerable and 80% resilient individuals. Wechange the probability of recovery of these two categories( P V ulR and P ResR ) with the constraint P V ulR ≤ P ResR . Ac-cordingly, we change these two probabilities from 0.6 to0.8 and 0.8 to 0.95 respectively. We choose these valuesaccording to reported case fatality ratios for the SARS-CoV-2 pandemic.[36, 38]
FIG. 9. Interdependence of different fractions of the pop-ulation as the immunity evolves. Percentage of immunized(colour coded) represented as a function of the percentageof survival for vulnerables and resilients. The proportions arewith respect to the total initial population. The primary vari-ables are the probabilities of recovery of the vulnerables andthe resilients. The results are obtained after averaging over100 simulations.
For every pair of P V ulR and P ResR we get a value of per-centage of vulnerables and resilients who survive and afraction of the population that gets immunized. In Fig-ure 9 we plot the survival % of vulnerables and resilientsin the two perpendicular axes and represent the % immu-nized as colour codes according to the colour gradationbar on the right hand side. In this contour representa-tion, red denotes low immunity and blue denotes higherimmunity.The survival % of the vulnerables is lower than that ofthe resilients. The percentage of immunized populationis higher (blue) for maximum survival of the resilients ascompared to that of the vulnerables. This means that toattain higher immunity in the population, greater num-ber of old and vulnerable people suffer death as comparedto resilients. Hence, attainment of herd immunity comeswith the cost of a higher mortality of the vulnerable sec-tion of the society.
IV. SUMMARY AND CONCLUSION
Any epidemic is a dynamic process where time depen-dence plays a crucial role in the control of the spreadand the damage, that is, the outcome. COVID-19 is apandemic which is currently under intense scrutiny byall and sundry, and many aspects are yet to be under-stood. Every move by the government, and the popula-tion in general, is of crucial importance. Each pandemiccomes with unique characteristics that deserve specialtreatments, not just medical and clinical but also socio-logical. In each such epidemic, immunity plays a criticalrole. Spanish Flu mainly attacked the age group between20 and 30 years of age. This is the age group with max-imum immunity. In the case of COVID-19, again weface the sad reality that certain section of the society issubstantially more vulnerable than other sections. Thevulnerable section consists of age groups which are above60-65 years of age, and people with comorbidity. Thereis yet to further classification, although it is conceivablethat as we understand the disease better and more pre-cisely, better perception of danger would emerge.An epidemic often starts by a process of nucleationwhich is an important phenomenon often studied inphysics and chemistry. The process of nucleation is initi-ated by a sudden appearance of a group of infected indi-viduals in a region. This may be triggered a laboratoryaccident, 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 ge-ography and demography of the country or region. Theinitial period of the process is often slow. After the ini-tial nucleation, the disease spreads by a diffusion processinto the susceptible population. Hence, it is a percolation with a temporal evolution.In order to address the issue of vulnerability of thepopulation and the outcome with the progression of theepidemic, we carry out a theoretical analysis with theobjective to analyze the consequences of aiming for herdimmunity without vaccine, or a good drug, in the con-text of the present COVID-19 pandemic. We developand solve a modified SIR model numerically and by em-ploying cellular automata simulations. We particularlyprobed the following question: what is dependence ofmortality on the rate of herd immunity?One of the key results of the present study is the depen-dence of the percentage survival on the rate of attainmentof the immunity threshold. We find that a late attain-ment of the immunity saturation leads to relatively lesserfatality. We show that approximately 50-60% of the vul-nerables might lose their lives in order to attain 70%total immunized population. On the contrary the mor-tality of the resilient fraction of population is relativelylow, may be just about 10%. We find a non-linear trendin the dependence of the cured and dead population onthe initial population of the vulnerables. This is becauseas the number of vulnerables increases, the immunity byinfection from a larger fraction of population which can-not protect the vulnerables unless deliberate efforts aremade that requires intervention.While we discuss herd immunity by infection in thiswork, the other, more sustainable option is herd immu-nity by vaccination. For example, diseases like small pox,polio etc. have been wiped off the face of earth by vacci-nation. This is particularly crucial for diseases with highmortality rates. However, for any novel disease, prepa-ration of a vaccine can take years. In case of the presentCOVID-19 pandemic, for instance, extensive research isgoing on globally in search of a vaccine.[1] However, nopromising result has been obtained in almost five monthsand researchers believe it may take more than a year toprepare the vaccine.
ACKNOWLEDGMENTS
We thank Prof. Sarika Bhattacharyya (NCL, Pune)and Prof. Suman Chakrabarty (SNCNCBS, Kolkata) forseveral fruitful discussions and comments. The authorsthank the Department of Science and Technology (DST),India for financial support. BB thanks sir J. C. Bose fel-lowship for partial financial support. S. Mo. thanks Uni-versities Grants Commission (UGC), India for researchfellowship. S. Mu. thanks DST-INSPIRE programmefor providing research fellowship. [1] W.-H. Chen, U. Strych, P. J. Hotez, and M. E. Bottazzi,Current tropical medicine reports , 1 (2020).[2] K. Prem, Y. Liu, T. W. Russell, A. J. Kucharski, R. M.Eggo, N. Davies, S. Flasche, S. Clifford, C. A. Pearson, J. D. Munday, et al. , The Lancet Public Health (2020).[3] M. Shen, Z. Peng, Y. Guo, Y. Xiao, and L. Zhang,medRxiv (2020).[4] Y. Kamikubo and A. Takahashi, medRxiv (2020). [5] P. E. Fine, Epidemiologic reviews , 265 (1993).[6] R. M. Anderson and R. M. May, Nature , 323 (1985).[7] T. J. John and R. Samuel, European journal of epidemi-ology , 601 (2000).[8] N. T. Georgette, PloS one , e4168 (2009).[9] E. S. McBryde, Clinical infectious diseases , 685(2009).[10] D. Wrapp, N. Wang, K. S. Corbett, J. A. Goldsmith, C.-L. Hsieh, O. Abiona, B. S. Graham, and J. S. McLellan,Science , 1260 (2020).[11] R. Singh and R. Adhikari, arXiv preprintarXiv:2003.12055 (2020).[12] S. Mukherjee, S. Mondal, and B. Bagchi, arXiv preprintarXiv:2004.14787 (2020).[13] D. J. Daley and J. Gani, Epidemic modelling: an intro-duction , Vol. 15 (Cambridge University Press, 2001).[14] W. O. Kermack and A. G. McKendrick, Proceedings ofthe royal society of london. Series A, Containing papers ofa mathematical and physical character , 700 (1927).[15] A. Skvortsov, R. Connell, P. Dawson, and R. Gailis, in
MODSIM 2007 International Congress on Modelling andSimulation. Modelling and Simulation Society of Aus-tralia and New Zealand (Citeseer, 2007) pp. 657–662.[16] D. S. Jones, M. Plank, and B. D. Sleeman,
Differentialequations and mathematical biology (CRC press, 2009).[17] R. M. Anderson and R. M. May, Nature , 361 (1979).[18] K. Dietz, Statistical methods in medical research , 23(1993).[19] O. Diekmann, J. A. P. Heesterbeek, and J. A. Metz,Publications of the Newton Institute , 95 (1995).[20] S. Zhang, M. Diao, W. Yu, L. Pei, Z. Lin, andD. Chen, International Journal of Infectious Diseases ,201 (2020).[21] B. Tang, N. L. Bragazzi, Q. Li, S. Tang, Y. Xiao, andJ. Wu, Infectious disease modelling , 248 (2020).[22] J. Yang, Y. Zheng, X. Gou, K. Pu, Z. Chen, Q. Guo,R. Ji, H. Wang, Y. Wang, and Y. Zhou, InternationalJournal of Infectious Diseases (2020).[23] J. Wangping, H. Ke, S. Yang, C. Wenzhe, W. Sheng- shu, Y. Shanshan, W. Jianwei, K. Fuyin, T. Penggang,L. Jing, et al. , Frontiers in Medicine , 169 (2020).[24] Q. Lin, S. Zhao, D. Gao, Y. Lou, S. Yang, S. S. Musa,M. H. Wang, Y. Cai, W. Wang, L. Yang, et al. , Interna-tional journal of infectious diseases (2020).[25] C. A. Hollingsworth, P. G. Seybold, L. B. Kier, and C.-K. Cheng, International journal of chemical kinetics ,230 (2004).[26] P. G. Seybold, L. B. Kier, and C.-K. Cheng, The Journalof Physical Chemistry A , 886 (1998).[27] S. Wolfram, Reviews of modern physics , 601 (1983).[28] M. Bartolozzi and A. W. Thomas, Physical review E ,046112 (2004).[29] B. S. Soares-Filho, G. C. Cerqueira, and C. L. Pen-nachin, Ecological modelling , 217 (2002).[30] A. Goltsev, F. De Abreu, S. Dorogovtsev, and J. Mendes,Physical Review E , 061921 (2010).[31] R. M. Almeida and E. E. Macau, in Journal of Physics:Conference Series , Vol. 285 (IOP Publishing, 2011) p.012038.[32] S. Fu and G. Milne, in
Proc. of the Australian Conferenceon Artificial Life (2003).[33] S. H. White, A. M. Del Rey, and G. R. S´anchez, AppliedMathematics and Computation , 193 (2007).[34] G. C. Sirakoulis, I. Karafyllidis, and A. Thanailakis,Ecological Modelling , 209 (2000).[35] A. Remuzzi and G. Remuzzi, The Lancet (2020).[36] S. Ruan, The Lancet Infectious Diseases (2020).[37] J. T. Wu, K. Leung, M. Bushman, N. Kishore, R. Niehus,P. M. de Salazar, B. J. Cowling, M. Lipsitch, and G. M.Leung, Nature Medicine , 1 (2020).[38] R. Verity, L. C. Okell, I. Dorigatti, P. Winskill, C. Whit-taker, N. Imai, G. Cuomo-Dannenburg, H. Thompson,P. G. Walker, H. Fu, et al. , The Lancet Infectious Dis-eases (2020).[39] L. Fang, G. Karakiulakis, and M. Roth, The Lancet.Respiratory Medicine , e21 (2020).[40] E. Livingston and K. Bucher, Jama323