Herd immunity under individual variation and reinfection
aa r X i v : . [ q - b i o . P E ] J u l A SEIR MODEL WITH VARIABLE SUSCEPTIBILITY OR EXPOSURETO INFECTION
ANTONIO MONTALB ´AN AND M. GABRIELA M. GOMES
Abstract.
We study a SEIR model considered by Gomes et. al. [3] and Aguas et. al. [1]where different individuals are assumed to have different levels of susceptibility or exposureto infection. Under this heterogeneity assumption, epidemic growth is effectively suppressedwhen the percentage of population having acquired immunity surpasses a critical level - theherd immunity threshold - that is much lower than in homogeneous populations. We findexplicit formulas to calculate herd immunity thresholds and stable configuration, and explorevariations of the model. Introduction
We analyse a SEIR (Susceptible-Exposed-Infectious-Recovered) model considered in [3]where each of the compartments S , E , I , and R are split into continuum many compart-ments S ( x ), E ( x ), I ( x ), and R ( x ) for x ∈ R + . We are modeling a situation where eachindividual has a level of susceptibility x . This individual will start in compartment S ( x ) andstay within the compartments S ( x ), E ( x ), I ( x ), and R ( x ) the whole time. They may infector be infected by individuals in other compartments. We will consider two types of models:In the variable susceptibility case the susceptibility of an individual at level x is propor-tional to x , or, in other words, if you compare an individual at level x and an individual at level y , the one at level x is xy times more likely to get infected than the one with susceptibility y .One may interpret this as variation in biological susceptibility which may be due to genetics,age, general health, etc.In the variable connectivity case the propensities for an individual at level x to acquireinfection and transmit to others are both proportional to x , or, in other words, if you comparean individual at level x and an individual at level y , the one at level x is xy times more likelyto get infected than the one in level y and also xy times more likely to infect someone else.One may interpret this as variation due to the connectivity, i.e., individuals that have manycontacts are both more likely to get infected and to infect others.For each x , we have a system of the form:S(x) xλ / / E(x) δ / / I(x) γ / / R(x)where λ is the force of infection which is formulated differently in the variable susceptibility or the variable connectivity cases:Variable susceptibility: λ = β Z I ( x ) dx, Variable connectivity: λ = β Z x I ( x ) dx. Saved: June 28, 2020Compiled: August 5, 2020
The system is given by the equations:˙ S ( x ) = − xλS ( x ) , (1) ˙ E ( x ) = xλS ( x ) − δE ( x ) , (2) ˙ I ( x ) = δE ( x ) − γI ( x ) , (3) ˙ R ( x ) = γI ( x ) . (4)We assume that the system has been scaled such that the total population is 1. The initialcondition satisfies S ( x ) = (1 − ǫ ) q ( x ), E ( x ) = ǫq ( x ) and I ( x ) = R ( x ) = 0, where q ( x ) is adistribution with mean 1 and coefficient of variation CV , and 0 < ǫ ≪ t to denote the time at which we are considering these compartments: S t ( x ), E t ( x ), I t ( x ), R t ( x ), and R eff ( x ). We use S t to denote the conjunction of all the com-partments S t ( x ) for x ∈ R + . We thus have S t = R + ∞ S t ( x ) dx . Same with E t , I t , R t , and R eff .We will use the first three momenta of S t ( x ), that we denote S t , ¯ S t and S t : S t = Z S t ( x ) dx, ¯ S t = Z xS t ( x ) dx, and S t = Z x S t ( x ) dx. When infection is absent ( ǫ = 0), we have S = 1, ¯ S = 1 and S = 1 + CV . But notethat S t ( x ) is not a probability density function for ǫ > S t becomes less than 1. Thequotient S t ( x ) /S t will be a probability distribution for ǫ > t with first and secondmomenta ¯ S t /S t and S t /S t which decrease as time passes. We will see that in the case wherethe initial configuration q ( x ) is a gamma distribution, all the distributions S t ( x ) /S t will bealso be gamma distributions with the same coefficient of variation CV but with lower mean.Similarly, we define the moments R eff , ¯ R t and R t for the recovered compartment, and thesame with E and I . Notice for instance that λ t is equal to βI t and to β ¯ I t in the variablesusceptibility and variable connectivity cases respectively.We obtain precise formulas for the effective reproductive number and the herd immunitythreshold in terms of the momenta of the susceptible population. We will see that in the casewhen q ( x ) is a gamma distribution, we can calculate this momenta in terms of S t and get anexact formula for the herd immunity threshold (HIT) in terms only of the basic reproductivenumber R and CV .In the variable susceptibility case we will get that R = βγ and R eff = ¯ S t βγ . This implies that herd immunity is achieved when ¯ S t < / R . If we assume that q ( x ) is agamma distribution the proportion of individuals that have been infected by the time herdimmunity is reached is HIT = 1 − R − CV . In the variable connectivity case we will get that R = (1 + CV ) βγ and R eff = S t βγ . This implies that herd immunity is achieved when S t < (1 + CV ) / R . If we assume that q ( x ) is a gamma distribution the proportion of individuals that have been infected by the time SEIR MODEL WITH VARIABLE SUSCEPTIBILITY OR EXPOSURE TO INFECTION 3 herd immunity is reached is
HIT = 1 − R − CV . We will call this model, the basic model . We will see now a few variations.1.1.
The model with reinfections.
In [1], the authors consider an extension of the modelwhere immunity after recovery is not fully protective, but only partially. A factor σ is added torepresent the quotient of the probability of getting reinfected after recovery over the probabilityof getting infected while fully susceptible.The model now looks like this:S(x) xλ / / E(x) δ / / I(x) γ / / R(x) σxλ f f with λ as above. The extended system is given by the equations:˙ S ( x ) = − xλS ( x ) , (5) ˙ E ( x ) = xλ ( S ( x ) + σR ( x )) − δE ( x ) , (6) ˙ I ( x ) = δE ( x ) − γI ( x ) , (7) ˙ R ( x ) = γI ( x ) − σxλR ( x ) . (8)The basic reproductive number is calculated exactly as in the no-reinfections case, but theeffective reproductive now depends not only on the distribution of S t ( x ), but also on thedistribution of R eff ( x ). We will see that • R eff = ( ¯ S t + σ ¯ R t ) · β/γ in the variable susceptibility case, and • R eff = ( S t + σR t ) · β/γ in the variable connectivity case.The system exhibits newer dynamics in comparison with the basic case. Depending onwhether σ is below or above 1 / R (known as the reinfection threshold [5, 4]) we get thateither the disease dies out after a while and a certain proportion of the population nevergets infected, or continues endemically and every individual is eventually infected and thenreinfected over and over again.1.2. A variation with a carrier stage.
In the original Gomes et.al. model [3], the exposedcompartments are not simply a latent stage but a carrier stage where individuals are infectiousbut to a lesser degree than individual in the infectious compartments. What changes is theforce of infection λ :Variable susceptibility: λ t = β Z ρE t ( x ) + I t ( x ) dx = β ( ρE t + I t ) , Variable connectivity: λ t = β Z x ( ρE t ( x ) + I t ( x )) dx = β ( ρ ¯ E t + ¯ I t ) . We will then get that the basic and effective reproductive numbers are as follows:Variable susceptibility: R = β (cid:18) ρδ + 1 γ (cid:19) and R eff = ( ¯ S t + σ ¯ R t ) β (cid:18) ρδ + 1 γ (cid:19) , Variable connectivity: R = (1 + CV ) β (cid:18) ρδ + 1 γ (cid:19) and R eff = ( S t + σR t ) β (cid:18) ρδ + 1 γ (cid:19) . ANTONIO MONTALB ´AN AND M. GABRIELA M. GOMES
In analogy with the previous models, we still get R eff = ( ¯ S t + σ ¯ R t ) · R and R eff = ( S t + σR t ) · R / (1 + CV ) respectively. The formulas for the herd immunity threshold in terms of R remain the same.1.3. A variation with a death rate.
Aguas et.al. [1] have one more feature: a death rate φ meaning that a proportion φ of the individuals that come out of I ( x ) go to a new compartment D . The model now looks like this:S(x) xλ / / E(x) δ / / I(x) (1 − φ ) γ / / φγ " " ❊❊❊❊❊❊❊❊❊ R(x) σxλ x x DWe will not analyse this model in detail. Let us just say that if the reinfection parameter σ is below the reinfection threshold 1 / R , the infection goes extinct slightly faster than before,and if the reinfection parameter is above the reinfection threshold, everybody eventually diesdue to the infection.2. Effective reproductive number in the basic model
Let us start by studying the basic model with no reinfections, no carrier stage, and no deathrate.
Definition 2.1.
The effective reproductive number R eff at time t is defined to be the averagenumber of secondary infections caused by an infected individual over their entire infectiousperiod.We make an approximation by assuming that while the individual is contagious, the sus-ceptible population does not change. That is, we disregard the fact that since the susceptiblepopulation declines, this individual infects more people at the beginning than at the end oftheir infection. The decline in the susceptible population is usually slow enough compared tothe length of the infectious period so that this does not make a big difference. Formally, wewill use R eff to analyze stable configurations, where this assumption holds, so our results willbe precise.Here we provide the derivations of the formulas presented in the introduction.Let’s look first at the variable susceptibility case: Consider an individual who getsinfected at time t (i.e., moving from S to E at time t ). They will eventually move to I wherethey will then spend, on average 1 /γ days. While in I , they will infect β R yS t ( y ) dy otherindividuals each day. We thus get R eff = 1 γ β (cid:18)Z y ( S t ( y )) dy (cid:19) = ¯ S t · βγ . In particular, we get that R = βγ and that R eff = ¯ S t · R .Let’s look now at the variable connectivity case: Consider an individual who gets infectedat time t (i.e., that enters E at time t ). It matters now what kind of individual they are, i.e.,what connectivity level they have, because individuals at different levels will infect differentnumbers of people. SEIR MODEL WITH VARIABLE SUSCEPTIBILITY OR EXPOSURE TO INFECTION 5
Let p ( x ) be the probability distribution function measuring the probability that this indi-vidual has connectivity level x . Their probability of becoming infected (i.e., of entering the E ( x ) compartment) is xλ . Thus, the value of p ( x ) is proportional to xS t ( x ). We get p ( x ) = xS t ( x ) R yS t ( y ) dy = x S t ( x )¯ S t . As above, an individual who enters E will eventually move to I where they are then going tospend, on average, 1 /γ days, and where they will infect β R yS t ( y ) dy other individuals eachday. We thus get R eff = Z (cid:16) x p ( x ) (cid:17)(cid:16) Z yS t ( y ) dy (cid:17) β/γ dx = Z (cid:16) x S t ( x )¯ S t (cid:17) ¯ S t β/γ dx. Moving things around we get R eff = S t · β/γ. In particular we get that R = q · β/γ = (1 + CV ) · β/γ and that R eff = S t · R / (1 + CV ).3. Herd immunity
Suppose we have a population with no infected individuals, so that all individuals are eithersusceptible or recovered. We say that this population has herd immunity if an introductionof the disease (i.e., a tiny increase in E ) does not trigger an epidemic. Mathematically, thismeans that any small enough deviation from the configuration with no infected individualswill quickly converge back to a configuration with no infected individuals. More formally, aconfiguration with E ( x ) = I ( x ) = 0 has herd immunity if, for every ǫ >
0, there is a δ > δ from the originalconfiguration, the system then converges to a configuration that is at a distance less than ǫ . If we visualize the dynamical system as modeling individuals who move between the com-partments S ( x ), E ( x ), I ( x ) and R ( x ), it is not hard to see that a configuration with E ( x ) = I ( x ) = 0 has herd immunity if and only if R eff ≤ . We say that a configuration with no infected individuals is at the herd immunity threshold if R eff = 1. Usually, in SEIR models with no variability, configurations are determined by thevalue of 1 − S = R , which is the number of people in the recovered compartment, and the herdimmunity threshold is defined as the value of 1 − S , the unique configuration that is at theherd immunity threshold, a value that is well-know to be equal to 1 − / R [2]. But, for theheterogenous models, there are many configurations which are at the herd immunity threshold.One such configuration is given by S ( x ) = q ( x ) / R for all x . This would be obtained, forinstance, if one vaccinates a proportion 1 − / R of the total population randomly withouttaking into account susceptibility levels. When immunity is acquired naturally, however,individuals with higher susceptibility get infected earlier and a configuration which has herdimmunity is reached before a proportion 1 − / R of the total population is infected.For now, let us observe that in the variable susceptibility case, herd immunity is achievedwhen ¯ S t = 1 / R , and that in the variable connectivity case, herd immunity is reached when S t = (1 + CV ) / R For a distance in the space of configurations, one may use the L distance, i.e., the sum of the integrals ofthe differences of the compartments R | S ( x ) − S ′ ( x ) | + | E ( x ) − E ′ ( x ) | + | I ( x ) − I ′ ( x ) | + | R ( x ) − R ′ ( x ) | dx . ANTONIO MONTALB ´AN AND M. GABRIELA M. GOMES
Next, we will see how we can calculate ¯ S t and S t under the assumption that the originaldistribution is a gamma distribution.4. Starting with the gamma distribution
Let us start this section studying how the distribution of susceptible compartments evolves,and then see how this fits nicely in the case where individual variation in susceptibility orconnectivity is gamma distributed.4.1.
The evolution of the susceptible compartments.
We claim that at every t , thereis a number k t ∈ R that only depends on t , such that(9) S t ( x ) = q ( x ) · e − x · k t . This holds in all models: in both the variable susceptibility case and the variable connectivitycase (with different values for k t ), in the case with reinfection, with carrier stage, and with adeath rate.Let’s start by proving equation (9). From the SEIR equation for ∂S t ( x ) ∂t we get1 S t ( x ) ∂S t ( x ) ∂t = x · λ t . Integrate with respect to t : Z t S t ( x ) ∂S t ( x ) ∂t dt = Z t x · λ t dt. Then by substitution Z t S t ( x ) dS t ( x ) = x · Z t λ t dt. Evaluating the integrals and letting k t = − R t λ t dt :ln( S t ( x )) − ln( q ( x )) = − xk t from which (9) follows.4.2. The gamma distribution.
We use the following notation for the gamma probabilitydensity function : Gamma a,b ( x ) = b a Γ( a ) x a − e − bx . The gamma distribution has mean a/b and coefficient of variation CV = 1 / √ a .Since our initial distribution q ( x ) has mean 1, we are using a gamma distribution with a = b , i.e., q ( x ) = Gamma a,a ( x ).Substituting q ( x ) for Gamma a,a ( x ) in (9) we get S t ( x ) = a a Γ( a ) x a − e − x ( a + k t ) = (cid:18) aa + k t (cid:19) a ( a + k t ) a Γ( a ) x a − e − x ( a + k t ) S t ( x ) = (cid:18) aa + k t (cid:19) a · Gamma a,a + k t ( x ) SEIR MODEL WITH VARIABLE SUSCEPTIBILITY OR EXPOSURE TO INFECTION 7
Using that Z Gamma a,a + k t ( x ) dx = 1 Z x Gamma a,a + k t ( x ) dx = aa + k t Z x Gamma a,a + k t ( x ) dx = a ( a + 1)( a + k t ) , we calculate S t , ¯ S t and S t : S t = (cid:18) aa + k t (cid:19) a · Z Gamma a,a + k t ( x ) dx = (cid:18) aa + k t (cid:19) a ¯ S t = (cid:18) aa + k t (cid:19) a · Z x Gamma a,a + k t ( x ) dx = (cid:18) aa + k t (cid:19) a · aa + k t = (cid:18) aa + k t (cid:19) a +1 S t = (cid:18) aa + k t (cid:19) a · Z x Gamma a,a + k t ( x ) dx = (cid:18) aa + k t (cid:19) a · a ( a + 1)( a + k t ) = (cid:18) aa + k t (cid:19) a +2 a + 1 a From the above we get ¯ S t = ( S t ) a +1 a and S t = ( S t ) a +2 a · a + 1 a . Using that a = 1 /CV , we get ( a + 1) /a = 1 + CV , and we can rewrite these as:¯ S t = ( S t ) CV and S t = ( S t ) CV · (1 + CV ) . Putting all together.
Recall that the herd immunity threshold (HIT) is 1 − S t at thetime t when R eff = 1 in the case where there are almost no people infected.In the variable susceptibility case we had that herd immunity is achieved when R eff =¯ S t · R . Thus, when R eff = 1, we have R − = ¯ S t = ( S t ) CV . It follows that
HIT = 1 − R − CV . In the variable connectivity case we had that herd immunity is achieved when R eff = S t R / (1 + CV ). Thus, when R eff = 1, we get R − = S t / (1 + CV ) = ( S t ) CV . It followsthat HIT = 1 − R − CV . The model with reinfection
Effective reproductive number.
Let’s consider now the model with a reinfection fac-tor σ as described in the introduction. We will see that R eff depends not only on S t ( x ) butalso on R t ( x ). When we consider configurations with no infected individuals, we will havethat R t ( x ) = q ( x ) − S t ( x ) and we will be able to express R eff in terms of S t ( x ) only.The formulas for the effective reproductive number R eff at time t are • R eff = ( ¯ S t + σ ¯ R t ) · β/γ in the variable susceptibility case, and • R eff = ( S t + σR t ) · β/γ in the variable connectivity case. ANTONIO MONTALB ´AN AND M. GABRIELA M. GOMES
In particular, we get R = β/γ and R = (1 + CV ) · β/γ as in the case with no reinfection.The derivation of these formulas is essentially the same as the derivations in Section 2. Thereare two differences: One, that each level- x individual infects β · ( ¯ S t + σ ¯ R t ) or xβ · ( ¯ S t + σ ¯ R t )other people in each day spent in I , in the respective cases, instead of β · ( ¯ S t ) or xβ · ( ¯ S t )Two, that when we consider an individual that gets infected in the variable connectivity case,the probability that they are a level- x individual is proportional to xS ( x ) + σxR ( x ) insteadof xS ( x ).5.2. Herd immunity.
Recall that a configuration with no infected individuals has herd im-munity if and only if R eff ≤ . Assuming that no one is infected, that is R ( x ) = q ( x ) − S ( x ), weget ¯ R = 1 − ¯ S and R = q − S . We can then understand the configurations at herd immunityin terms of ¯ S and S .In the variable susceptibility case a configuration with no infected individuals is at theherd immunity threshold if and only if 1 / R = ¯ S + σ (1 − ¯ S ), and hence, if and only if¯ S = R − − σ − σ . In the variable connectivity case a configuration with no infected individuals is at theherd immunity threshold if and only if 1 / R = ( S + σ ( q − S )) / (1 + CV ), and hence, if andonly if S = (1 + CV ) R − − σ − σ . A point worth mentioning is that the behavior of herd immunity differ from the basic model.When one runs the basic model, herd immunity is achieved when the number of infectionsstart to decline. In the case with reinfections, the number of infections may decline earlier.The reason is that when E t and I t are non-zero, the individuals in those compartments arenot susceptible to get infected at that moment, but will become partially susceptible oncethey recover. In other words, we may have R eff < R t , we could create a situation where R eff >
1. In practical terms, imagine strongisolation measures are imposed at that moment and all those infected people recover withoutinfecting anyone else. Then, when the measures are lifted, since R t becomes larger, R eff islarger and any reappearance of the disease my trigger a small epidemic. Thus, we need thecondition R eff ≤ R t ( x ) = q ( x ) − S t ( x ) to ensure that new introduction willnot trigger epidemics.5.3. The reinfection threshold.
For the formulas above to make sense, it is necessary thatwe have σ < R − . That is, the reinfection factor σ has to be below R − , a critical value known as the reinfectionthreshold [5, 4]. If this is verified, then all configuration with no infected individuals andsatisfying the conditions above (either ¯ S = ( R − − σ ) / (1 − σ ) or S = (1+ CV )( R − − σ ) / (1 − σ )depending on the case) are configurations which have herd immunity in the sense that anyincrease of the infected compartments quickly extinguishes as R eff won’t go above 1 again.If the reinfection factor is above the reinfection threshold, R eff will be greater than 1 in anysuch configurations, so there won’t be any configuration with no infected individuals and withherd immunity. This implies that there will always be a portion of the population infected, andhence that the population of susceptible individuals will eventually be completely depleted.The equilibrium configuration will now have S = 0. SEIR MODEL WITH VARIABLE SUSCEPTIBILITY OR EXPOSURE TO INFECTION 9
Suppose we have an equilibrium configuration when σ > R − . We will then have S = 0, andhence E ( x ) + I ( x ) + R ( x ) = q ( x ) for all x . Setting ˙ I ( x ) to zero, we get that E ( x ) = I ( x ) · γ/δ and from this we derive that I ( x ) = ( q ( x ) − R ( x )) · δ/ ( δ + γ ) and E ( x ) = ( q ( x ) − R ( x )) · γ/ ( δ + γ ).Using that R eff = 1 at any equilibrium configuration, we get that either ¯ R = σ R or R = σ R .5.4. Starting with a gamma distribution.
Recall from Section 3 that if we start with agamma distribution for q ( x ), we get that S t ( x ) /S t remains a gamma distribution for all t , andthat ¯ S t = S CV t and S t = S CV t (1 + CV ). We can then obtain the values of S t at themoment when herd immunity is achieved, and then calculate HIT = 1 − S t .In the variable susceptibility case we get: HIT = 1 − (cid:18) R − − σ − σ (cid:19) CV . In the variable connectivity case we get:
HIT = 1 − (cid:18) R − − σ − σ (cid:19) CV . The model with a carrier stage
Recall that, in this model, individuals in the E compartments are infectious, but by afactor ρ with respect to the individuals in I . This model is sometimes used when one wantsto differentiate pre-symptomatic form symptomatic stages, which is usually closer to what thereal word data represents, and, at the same time, allow for pre-symptomatic transmission.The calculation of the effective reproductive number R eff is slightly different. The differenceis that now one has to add the time an individual is in E to the infectious period, multipliedby the factor ρ . The average time an individual spends in E is 1 /δ . We then get • R eff = ( ¯ S t + σ ¯ R t ) · β ( ρ/δ + 1 /γ ) in the variable susceptibility case, and • R eff = ( S t + σR t ) · β ( ρ/δ + 1 /γ ) in the variable connectivity case,where ¯ S t and S t are the momenta of S t ( x ) defined above, and the same with R t .In particular, we get R = β ( ρ/δ + 1 /γ ) and R = (1 + CV ) · β ( ρ/δ + 1 /γ ) respectivelyas in the previous case. Also as in the previous cases we get R eff = ( ¯ S t + σ ¯ R t ) · R and R eff = ( S t + σR t ) · R / (1 + CV ). We then get the same formulas for the herd immunitythreshold in terms of R and CV as before. References [1]
R. Aguas , R. M. Corder , J. G. King , G. Gonc¸alves , M. U. Ferreira , M. G. M. Gomes ,Herd immunity thresholds for SARS-CoV-2 estimated from unfolding epidemics. medRxiv , medRxiv R. M. Anderson , R. M. May , Infectious Diseases of Humans: Dynamics and Control . Oxford UniversityPress Inc., New York (1991).[3]
M. G. M. Gomes , R. M. Corder , J. G. King , K. E. Langwig , C. Souto-Maior , J. Carneiro , G.Gonc¸alves , C. Penha-Gonc¸alves , M. U. Ferreira , R. Aguas , Individual variation in susceptibilityor exposure to SARS-CoV-2 lowers the herd immunity threshold. medRxiv
M. G. M. Gomes , E. Gjini , J. S. Lopes , C. Souto-Maior , C. Rebelo , A theoretical framework toidentify invariant thresholds in infectious disease epidemiology.
Journal of Theoretical Biology
M. G. M. Gomes , L. J. White , G. F. Medley , Infection, reinfection, and vaccination under suboptimalimmune protection: Epidemiological perspectives.
Journal of Theoretical Biology
Department of Mathematics, University of California, Berkeley, USA
E-mail address : [email protected] URL : ∼ antonio Department of Mathematics and Statistics, University of Strathclyde, Glasgow, United King-dom
E-mail address : [email protected]
URL ::