SEIR model

General equations

Clear[ages,stages,classes,b,freq];

We start by building code that allows an arbitrary number of age classes (j) and stage classes (jj), where stage refers to the activity level in the Britton et al. model.

Relative transmission rates (β) are assumed to depend on the ages of the two individuals involved (j,k), while the activity class depends on the individuals involved in a manner that potentially depends on its age, j, and distancing class, jj, captured by b[j,jj].  We will assume b is measured relative to the different classes, so that the average activity level is (OverBar[b])=1 across age and stage classes.

[In this notebook, we focus on only one age class.]

Social mitigation measures are incorporated by a constant term α reducing transmission rates ((Sqrt[α]) would describe the reduction in activity of each individual partner).

Disease transmission then depends on the efficacy of social mitigation measures (α), the activity of each partner (b[j,jj] b[k,kk]), the transmission rate between those partners (βI[j,k]), and the chance that one of them is susceptible and one is infected:

dsdt[j_,jj_]:=-Sum[α*b[j,jj] b[k,kk] (βI[j,k]*i[k,kk]),{k,1,ages},{kk,1,stages}]*s[j,jj];
dedt[j_,jj_]:=Sum[α*b[j,jj] b[k,kk] (βI[j,k]*i[k,kk]),{k,1,ages},{kk,1,stages}]*s[j,jj]-κE[j,jj]*e[j,jj];
didt[j_,jj_]:=κE[j,jj]*e[j,jj] - κI[j,jj] i[j,jj];

κE and κI are the transition rates from E to I and from I to outcome (recovery or death).

eqns:=Flatten[Table[{dsdt[j,jj],dedt[j,jj],didt[j,jj]},{j,1,ages},{jj,1,stages}]];

vars:=Flatten[Table[{s[j,jj],e[j,jj],i[j,jj]},{jj,1,stages},{j,1,ages}]];

start:=Flatten[Table[{s[j,jj]->s0[j,jj],e[j,jj]->0,i[j,jj]->0},{j,1,ages},{jj,1,stages}]];

We will want to relate the heterogeneity in movement rates to the mean (OverBar[b]) and variance (Subscript[V, b]) in b.  For the two activity classes, these are:

subb=Simplify[Solve[{(OverBar[b])==b1 freq[1]+b2 freq[2],(Subscript[V, b])==(((b1-(OverBar[b])))^2) freq[1]+(((b2-(OverBar[b])))^2) freq[2]},{b1,b2}],{ (freq[1]+freq[2])==1}]

We’ll focus on the second of these solutions, where type 1 is more active (b1 > b2):

subb=subb[[2]];

Note that activity is a relative measure with (OverBar[b])=1, which constrains 0≤(Subscript[V, b])≤1 when freq[1]=freq[2] because the two activities must be positive:
	{b1,b2}/.subb/.(OverBar[b])->1/.freq[1]->1/2/.freq[2]->1/2
	{2 ((1/2)+((Sqrt[(Subscript[V, b])])/2)),1-(Sqrt[(Subscript[V, b])])}

One class

ages=1; (* Age class *)
stages=1;(* Distancing class *)
classes=3; (*SEI*)

WIth only one distancing class, we can set all the b[j,jj] terms to 1:

freq[1]=1; (*Frequency of distancing class one*)
b[j_,jj_]=1; (*Relative activity is one with only one activity class*)

We will measure the spread of the disease in an initial population composed of a fraction, p, of susceptibles in the population (i.e., not vaccinated or previously infected and resistant).

The stability matrix will determine the spread of the disease:

stabmat=Table[D[eqns[[j]],vars[[k]]],{j,1,Length[eqns]},{k,1,Length[eqns]}]/.start/.s0[j_,jj_]->p*S*freq[jj]/.S->1/.κE[j_,jj_]->κE/.κI[j_,jj_]->κI/.βI[j_,jj_]->βI;

Its characteristic polynomial will determine the eigenvalues:

charpoly1=Det[λ IdentityMatrix[Length[eqns]]-stabmat]//Factor

Holding r0 constant: In the absense of immunity and no social mitigation measures (all individuals initially susceptible, p=1 and α=1), we could estimate transmissibility βI from the observed exponential growth rate of r0 and hold this rate constant:

subβr=Flatten[Solve[((charpoly1/λ)/.λ->r0/.p->1/.α->1)==0,βI]]

We then consider the factor p by which the susceptible population would have to be reduced to cause the growth rate to equal zero:

Factor[((charpoly1/λ)/.λ->0/.α->1/.subβr)];

charpolythreshold1=p/.Flatten[Solve[%==0,p]]//Simplify

This equals 1/R0 as found below, which provides a clearer way of thinking about the level of resistance needed to prevent the spread of the disease:

1/R0/.R0->((βI)/(κI))/.subβr//Factor

Holding R0 constant: We find R0 using the next generation method (see Hurford et al. 2010 J. R. Soc. Interface 7:561–571 and references therein), which is useful in more complicated cases with abritrary age and stages. 

First, we must identify the F matrix (reproduction) and V matrix (movement), such that stabmat = F-V.  Because we assume a large population of susceptibles relative to disease classes, we can ignore the first row and column of the stability matrix and focus on the elements describing the disease types:

Fmatrix=Transpose[Transpose[(stabmat/.κE->0/.κI->0)[[2;;3]]][[2;;3]]]

Vmatrix=-(Transpose[Transpose[(stabmat)[[2;;3]]][[2;;3]]]-Fmatrix)//Factor

which satisfies the conditions required for the next generation method (see Hurford et al.). 

R0 is given by the largest real eigenvalue of the simpler matrix: F (V^(-1))

Fmatrix.Inverse[Vmatrix]//Simplify//MatrixForm

repronum1=%[[1,1]]

This fits the expectation that if everybody were susceptible (p=1) and moving around normally (α=1), then R0=((βI)/(κI)) is the transmissibility βI times the average amount of time spent in the infectious class (1/κI).

In order to have R0 fall to one without any mitigation measures (α=1), we would need p to equal the inverse of the R0 in a fully susceptible population:

subβR=Flatten[Solve[(repronum1/.p->1/.α->1)==R0,βI]]

Factor[repronum1/.α->1/.%]

R0threshold1=p/.Flatten[Solve[%==1,p]]//Simplify

With only one infectious class, we can use the equation (1) of Britton et al. to determine the social distancing (α<1) needed such that the fraction of the population, sinf, that would remain susceptible after the spread of a disease with social distancing falls to the threshold p found above (assuming that enough time has passed with social mitigation that the disease has become rare again):

λ in equation (1) is a typo and should be α???

eqn1={sinf==Exp[-α βI (1-sinf)/κI]};

α/.Solve[eqn1/.sinf->1/((βI)/(κI)),α][[1]]/.subβR//Simplify

%/.R0->2.5

This is consistent with Figure 2, where α set to 0.6 was not quite enough to control the disease (yellow curves rise after the first phase is over).

Two classes - Heterogeneity in fraction of susceptibles

Entering equations

Clear[b,freq]

ages=1; (* Age class *)
stages=2;(* Distancing class *)
classes=3; (*SEI*)

Now let us change the differential equations such that the first stage corresponds to lower movement and the second to higher movement individuals (reduced from the three stage model of Britton et al.).

Although we could allow movement from non-distancing to distancing groups at rate ud and from distancing to non-distancing at rate ur (constant for all disease classes and ages), we will assume these are zero, as in Britton et al.:

tab1=Table[({dsdt[j,jj],dedt[j,jj],didt[j,jj]}/.jj->1)-(ud*{s[j,jj],e[j,jj],i[j,jj]}/.jj->1)+(ur*{s[j,jj],e[j,jj],i[j,jj]}/.jj->2),{j,1,ages}]/.ur->0/.ud->0;

tab2=Table[({dsdt[j,jj],dedt[j,jj],didt[j,jj]}/.jj->2)+(ud*{s[j,jj],e[j,jj],i[j,jj]}/.jj->1)-(ur*{s[j,jj],e[j,jj],i[j,jj]}/.jj->2),{j,1,ages}]/.ur->0/.ud->0;

eqns=Flatten[Join[tab1,tab2]];

WIth two distancing classes, we set the b[j,1] terms to b1 (not social distancing) and the b[j,2] terms to b2 (b1>b2), reflecting the current extent of social distancing.

b[j_,1]=b1;
b[j_,2]=b2;

The stability matrix will determine the spread of the disease:

stabmat=Table[D[eqns[[j]],vars[[k]]],{j,1,Length[eqns]},{k,1,Length[eqns]}]/.start/.s0[j_,1]->p1*S*freq[1]/.s0[j_,2]->p2*S*freq[2]/.S->1/.κE[j_,jj_]->κE/.κI[j_,jj_]->κI/.βI[j_,jj_]->βI/.α[j_,jj_]->α;

Its characteristic polynomial will determine the eigenvalues:

charpoly2=Det[λ IdentityMatrix[Length[eqns]]-stabmat]//Factor

In the absense of immunity and no social mitigation measures (all individuals initially susceptible, p1=p2=1 and α=1), we could estimate transmissibility βI from the observed exponential growth rate of r0 and hold this rate constant (rather than R0, as Britton et al. held constant):

subβr=Simplify[Flatten[Solve[((charpoly2/(λ^2))/.λ->r0/.p1->1/.p2->1/.α->1)==0,βI]],{freq[1]+freq[2]==1}]

Note, however, that holding r0,κE,κI constant gives the same relationship between βI and the activity classes as estimated from R0 (see below; {βI->((R0 κI)/((b1^2) freq[1]+(b2^2) freq[2]))}), given that all activity classes have equal transition rates.  This would not be true in an age-dependent model, however, if more active individuals are also more likely to proceed through the asymptomatic route of infection.  In that case, one would have to be careful about what one holds constant (r0 or R0).

Determining level of “infleuncer immunity”

Britton et al. estimate herd immunity levels with and without heterogeneous contact rates, assuming that the number of infected individuals per newly infected individual, R0, is held constant.  If all individuals were immunized/resistant at the same rate (leaving p susceptible), then p*R0 would have to be less than one to prevent the disease from spreading. If, however, the fraction of susceptible individuals of type i is pi, then the disease can stop spreading at a lower fraction of resistant individuals (higher average fraction of susceptibles (OverBar[p])).  Here we solve for (OverBar[p]). To do so, Britton et al. identify the maximum amount of social activity (α) that will cause just enough resistance to build up in the first wave that a seccond wave does not occur once social mitigation measures are relaxed (returning α=1)

Rather than calling this “herd immunity”, it is easier to think of this as “influencer immunity”, where influencers are those individuals most likely to spread the disease (here more active individuals with more contacts). The “influencer immunity” is given by the fraction of individuals that would have to be resistant, 1-(OverBar[p]), to prevent a second wave and is less than the fraction needed if all individuals were equally resistant (= “herd immunity” level).

Holding R0 constant: We find R0 using the next generation method (see Hurford et al. 2010 J. R. Soc. Interface 7:561–571 and references therein), which is useful in more complicated cases with abritrary age and stages. 

First, we must identify the F matrix (reproduction) and V matrix (movement), such that stabmat = F-V.  Because we assume a large population of susceptibles relative to disease classes, we can ignore the first row and column of the stability matrix and focus on the elements describing the disease types.

We can use Drop to drop the first and fifth rows and columns (corresponding to the S class whose eigenvalues are 0 when the disease is rare):

Fmatrix=Transpose[Drop[Transpose[Drop[(stabmat/.κE->0/.κI->0),{1,4,3}]],{1,4,3}]]

Vmatrix=-(Transpose[Drop[Transpose[Drop[(stabmat),{1,4,3}]],{1,4,3}]]-Fmatrix)//Factor

Correctly gives the same eigenvalue (stabmat = Fmatrix-Vmatrix):

Det[λ IdentityMatrix[Length[Vmatrix]]-(Fmatrix-Vmatrix)];
Factor[charpoly2/%]

which satisfies the conditions required for the next generation method (see Hurford et al.). 

R0 is given by the largest real eigenvalue of the simpler matrix: F (V^(-1))

Fmatrix.Inverse[Vmatrix]//Simplify//MatrixForm

Eigenvalues[%]

repronum2=Last[%]

Simplify[%/.subb,{freq[1]+freq[2]==1}]

Here, we hold R0 constant within a population that is fully susceptible (p1=p2=1) and not socially distancing (α=1), allowing us to relate βI to R0:

subβR=Simplify[Flatten[Solve[(repronum2/.α->1/.p1->1/.p2->1)==R0,βI]],{freq[1]+freq[2]==1}]

We then consider a social mitigation measure (α<1) that allows the disease to run its course and then relax the mitigation (returning to α=1).  At this point, we ask if the fraction of susceptibles is sufficiently low to prevent a second wave of the disease.

In order to avoid a second wave without any mitigation measures (α=1), we need p1 and p2 to balance the reproductive number R0 in a fully susceptible population such that the following falls below one when mitigation measures are relaxed:

Simplify[repronum2/.α->1/.subβR,{freq[1]+freq[2]==1}]

That is, we must choose α such that the above falls just to one (no further growth), i.e., where (((b1^2) p1 freq[1]+(b2^2) p2 freq[2])/((b1^2) freq[1]+(b2^2) freq[2]))==(1/R0) (call this the “first wave target”).

The simplest scenario is when there is extreme heterogeneity, where b2=0, in this case, we just need the frequency of susceptibles in the high movement class p1 to equal (1/R0), while non-movers remain susceptible (p2=1).

(((b1^2) p1 freq[1]+(b2^2) p2 freq[2])/((b1^2) freq[1]+(b2^2) freq[2]))/.b2->0

-> The best case for a first wave eliminating the disease is when there is extreme variation in movement, with one class moving at rate b1 and the other class not moving. In this case, “influencer immunity” is reached when p1=(1/R0), such that (OverBar[p])=p1 freq[1]+p2 freq[2]=((freq[1])/R0)+freq[2].

One minus this fraction, (1-(1/R0))freq[1], is the fraction of influencers that must be resistant to prevent the disease from spreading during the second wave. If both classes are equally frequent (freq[1]=1/2), this is half the number of immune individuals compared to the homogenous case where all individuals move about equally.  Even less resistance is needed if the active class is rarer.  Recall, though, that we hold R0 constant, so if the active class is rare, infected individuals in this class must lead to many new infections.

This “best case” conclusion holds no matter how many classes, because the best case always occurs when there is the maximal possible variability, which occurs when some individuals move a lot and others not at all.

Numerical exploration suggests that roughly the same α value will lead to the best case first wave, (((b1^2) p1 freq[1]+(b2^2) p2 freq[2])/((b1^2) freq[1]+(b2^2) freq[2]))==(1/R0), regardless of how different the activity rates are betweent the two groups.  The fact that α remains similar, regardless of heterogeneity in activity, isn’t obvious and we return to this point below.

With two infectious classes, we can use equation (1) of Britton et al. to determine the social distancing (α<1) needed such that the fraction of the two sub-populations, p1=((s1[∞])/(s1[0])) and p2=((s2[∞])/(s2[0])), that remain susceptible after the spread of a disease with social distancing falls to this threshold (assuming that enough time has passed with social mitigation that the disease has become rare again).

This requires solving equation 1 of Britton et al. for the susceptible fractions, p1=((s1[∞])/(s1[0])) and p2=((s2[∞])/(s2[0])), giving two equations that must equal zero:

eqn1={p1-Exp[-α freq[1] βI (b1^2) (1-p1)/κI-α freq[2] βI b1 b2 (1-p2)/κI],
p2-Exp[-α freq[1] βI b1 b2 (1-p1)/κI-α freq[2] βI (b2^2) (1-p2)/κI]}//Simplify

Without any heterogeneity (b1=b2=1), we regain the result of the one class model that α->((Log[(1/R0)])/(1-R0)):

Simplify[eqn1[[1]]/.subβR/.b1->b/.b2->b/.p1->(OverBar[p])/.p2->(OverBar[p]),{freq[1]+freq[2]==1}]//Flatten

Simplify[Solve[%==0,α]]//Flatten

Simplify[%/.(OverBar[p])->1/R0]

%/.R0->2.5

This is consistent with Figure 2 in Britton et al., where α set to 0.6 was not quite enough to control the disease (yellow curves rise after the first phase is over).

With heterogeneity in activity, we can use the relationship between p1, p2 and R0 from the first wave target to replace for p2:

Simplify[Solve[((((b1^2) p1 freq[1]+(b2^2) p2 freq[2])/((b1^2) freq[1]+(b2^2) freq[2])))==(1/R0),p2],{freq[1]+freq[2]==1}]//Flatten

Simplify[eqn1/.subβR/.%,{freq[1]+freq[2]==1}]//Flatten

Solving the first one for α and plugging into the second:

Simplify[Solve[%[[1]]==0,α],{freq[1]+freq[2]==1}]//Flatten

Simplify[%%[[2]]/.%,{freq[1]+freq[2]==1}]

This cannot be solved explicitly for p1, but it can be rearranged more simply as:

Simplify[Solve[%==0,R0]]//Flatten

This provides a useful insight, as it indicates that the solution must have the form (((b1^2) p1 freq[1]+ (b2^2) (p1^(b2/b1)) freq[2])/((b1^2) freq[1]+(b2^2) freq[2]))==(1/R0), which is identical to the “first wave target”  that causes the disease to neither grow nor decline after mitigation, except with p2=(p1^(b2/b1)).

-> Thus, the fraction of susceptibles that must be resistant for mitigation to just stop a second wave once relaxed satisfies  (p1^(1/b1))=(p2^(1/b2)).

At this point, we have an equation that we can solve for p1 and then get p2, and we do not need to find α to do so:

Clear[getme]
getme[R0_,Var_,fr_]:=getme[R0,Var,fr]=FindRoot[((((b1^2) p1 freq[1]+ (b2^2) (p1^(b2/b1)) freq[2])/((b1^2) freq[1]+(b2^2) freq[2]))==(1/R0)/.subb/.(OverBar[b])->1/.(Subscript[V, b])->Var/.freq[1]->fr/.freq[2]->1-fr),{p1,0.2,0,1}]

For example, when there is no heterogeneity ((Subscript[V, b])=0), we regain the one class result that p1=p2=2/5 to balance an R0 of 5/2:

getme[2.5,0,1/2]

But with heterogeneity, the disease can stop spreading after mitigation is relaxed (α returning to one) with a higher total number of susceptible individuals, as argued by Britton et al.

Table of the final fraction of susceptibles with increasing variance in activity (from 0 to 1):

Table[(OverBar[p])/.(OverBar[p])->freq[1]p1+freq[2]p2/.p2->(p1^(b2/b1))/.subb/.(OverBar[b])->1/.(Subscript[V, b])->Var/.getme[2.5,Var,1/2]/.R0->5/2/.freq[1]->1/2/.freq[2]->1/2,{Var,0,1,1/10}]

That is, “influencer immunity” requires fewer resistant individuals than “herd immunity” with heterogeneity in activity, as illustrated in the following plot (red = total fraction of susceptibles, black is among more active “influencers”, and black dashed is among less active individuals):

Show[ListPlot[Table[{Var,(OverBar[p])/.(OverBar[p])->freq[1]p1+freq[2]p2/.p2->(p1^(b2/b1))/.subb/.(OverBar[b])->1/.(Subscript[V, b])->Var/.getme[2.5,Var,1/2]/.R0->5/2/.freq[1]->1/2/.freq[2]->1/2},{Var,0,1,1/100}],Joined->True,PlotStyle->Red],
ListPlot[Table[{Var,p1/.p2->(p1^(b2/b1))/.subb/.(OverBar[b])->1/.(Subscript[V, b])->Var/.getme[2.5,Var,1/2]/.R0->5/2/.freq[1]->1/2/.freq[2]->1/2},{Var,0,1,1/100}],Joined->True,PlotStyle->Black],ListPlot[Table[{Var,p2/.p2->(p1^(b2/b1))/.subb/.(OverBar[b])->1/.(Subscript[V, b])->Var/.getme[2.5,Var,1/2]/.R0->5/2/.freq[1]->1/2/.freq[2]->1/2},{Var,0,1,1/100}],Joined->True,PlotStyle->{Black,Dashed}],PlotRange->{{0,1},{0,1}},AxesOrigin->{0,0},AxesLabel->{"(Subscript[V, b])","Freq Susceptibles"}

]

Here we note that the fraction of individuals that are susceptible in the different activity classes satisfies (pi^(1/bi))=constant even with more activity classes.

We can verify that p2->(p1^(b2/b1)) has the correct form for the solution in that the two equations defined by Britton et al. become identical, if we rearrange and then take the two sides to the power of b1/b2:

(e^((b2 α βI (b1 (-1+p1) freq[1]+b2 (-1+(p1^(b2/b1))) freq[2]))/(κI)))==(p1^(b2/b1))

(e^((b1 α βI (b1 (-1+p1) freq[1]+b2 (-1+(p1^(b2/b1))) freq[2]))/(κI)))==p1

This also true with more heterogeneity classes, e.g., with three heterogeneity classes:

eqn1three={p1-Exp[-α freq[1] βI (b1^2) (1-p1)/κI-α freq[2] βI b1 b2 (1-p2)/κI-α freq[3] βI b1 b3 (1-p3)/κI],
p2-Exp[-α freq[1] βI b1 b2 (1-p1)/κI-α freq[2] βI (b2^2) (1-p2)/κI-α freq[3] βI b2 b3 (1-p3)/κI],
p3-Exp[-α freq[1] βI b1 b3 (1-p1)/κI-α freq[2] βI b2 b3 (1-p2)/κI-α freq[3] βI (b3^2) (1-p3)/κI]}//Simplify

A solution has the form (p1^(1/b1))==(p2^(1/b2))==(p3^(1/b3))

Simplify[eqn1three/.p2->(p1^(b2/b1))/.p3->(p1^(b3/b1))]

as can be confirmed again by rearranging the last two equations with the exponent on one side and the (pi^(bi/b1)) term on the other, taking both sides to the power of b1/bi, and noting that we regain the first equation.  The same is true regardless of the number of activity classes, as can be shown by plugging in pi=(p1^(bi/b1)) into equation (1) of Britton et al. pi=Exp[-α (∑_(j=1)^classes)freq[i] βI bi bj (1-pj)/κI].

Exploring why α remains roughly the same, regardless of heterogeneity in activity

Returning to why α is nearly the same value regardless ...

From the first equation in equation (1), we solved for α above.  Plugging in the solution for {p1,p2}, we see that α is nearly constant regardless of the variance:

Table[((b2 ((b1^2) freq[1]+(b2^2) freq[2]) Log[p1])/(b1 (-b1 b2 (-1+p1) R0 (-1+freq[2])+(b1^2) (-1+p1 R0) (-1+freq[2])-(b2^2) (-1+R0) freq[2])))/.p2->(p1^(b2/b1))/.subb/.(OverBar[b])->1/.(Subscript[V, b])->Var/.getme[2.5,Var,1/2]/.R0->5/2/.freq[1]->1/2/.freq[2]->1/2,{Var,0,1,1/10}]

ListPlot[Table[{Var,((b2 ((b1^2) freq[1]+(b2^2) freq[2]) Log[p1])/(b1 (-b1 b2 (-1+p1) R0 (-1+freq[2])+(b1^2) (-1+p1 R0) (-1+freq[2])-(b2^2) (-1+R0) freq[2])))/.(OverBar[p])->freq[1]p1+freq[2]p2/.p2->(p1^(b2/b1))/.subb/.(OverBar[b])->1/.(Subscript[V, b])->Var/.getme[2.5,Var,1/2]/.R0->5/2/.freq[1]->1/2/.freq[2]->1/2},{Var,0,99/100,1/100}],Joined->True,PlotStyle->Blue,AxesLabel->{"(Subscript[V, b])","α"}]

This slight change in the necessary level of mitigation is real (checked and is not a numerical inaccuracy), but does it matter that α isn’t exactly constant at the level α=((Log[p1])/(1-R0)) predicted in the one-class case?

Here we obtain the predicted values of p1=((s1[∞])/(s1[0])) and p2=((s2[∞])/(s2[0])) from equation 1 of Britton et al. assuming that α was given by the mitigation level needed to prevent a second wave in a homogeneous population, α=((Log[1/R0])/(1-R0)):

Simplify[eqn1/.α->((Log[1/R0])/(1-R0))/.subβR,{freq[1]+freq[2]==1}]

Again, after the first wave, the two equations are equivalent if p2=(p1^(b2/b1)), leaving:

Simplify[%[[1]]/.p2->(p1^(b2/b1))]

Clear[getme2]
getme2[R0_,Var_,fr_]:=getme2[R0,Var,fr]=FindRoot[(p1-((((1/R0)))^(-((b1 R0 (b1 (-1+p1) freq[1]+b2 (-1+(p1^(b2/b1))) freq[2]))/((-1+R0) ((b1^2) freq[1]+(b2^2) freq[2])))))==0/.subb/.(OverBar[b])->1/.(Subscript[V, b])->Var/.freq[1]->fr/.freq[2]->1-fr),{p1,0.2,0,1}]

Show[ListPlot[Table[{Var,(OverBar[p])/.(OverBar[p])->freq[1]p1+freq[2]p2/.p2->(p1^(b2/b1))/.subb/.(OverBar[b])->1/.(Subscript[V, b])->Var/.getme[2.5,Var,1/2]/.R0->5/2/.freq[1]->1/2/.freq[2]->1/2},{Var,0,1,1/100}],Joined->True,PlotStyle->Red],
ListPlot[Table[{Var,p1/.p2->(p1^(b2/b1))/.subb/.(OverBar[b])->1/.(Subscript[V, b])->Var/.getme[2.5,Var,1/2]/.R0->5/2/.freq[1]->1/2/.freq[2]->1/2},{Var,0,1,1/100}],Joined->True,PlotStyle->Black],ListPlot[Table[{Var,p2/.p2->(p1^(b2/b1))/.subb/.(OverBar[b])->1/.(Subscript[V, b])->Var/.getme[2.5,Var,1/2]/.R0->5/2/.freq[1]->1/2/.freq[2]->1/2},{Var,0,1,1/100}],Joined->True,PlotStyle->{Black,Dashed}],
ListPlot[Table[{Var,(OverBar[p])/.(OverBar[p])->freq[1]p1+freq[2]p2/.p2->(p1^(b2/b1))/.subb/.(OverBar[b])->1/.(Subscript[V, b])->Var/.getme2[2.5,Var,1/2]/.R0->5/2/.freq[1]->1/2/.freq[2]->1/2},{Var,0,1,1/100}],Joined->True,PlotStyle->Purple],
PlotRange->{{0,1},{0,1}},AxesOrigin->{0,0},AxesLabel->{"(Subscript[V, b])","Freq Susceptibles"}

]

There is a slight but negligible difference between the two choices (setting α to the maximum possible value to avoid the second wave knowing the heterogeneity in red, and not knowing the heterogeneity in purple), the difference being:

ListPlot[Table[{Var,((OverBar[p])/.(OverBar[p])->freq[1]p1+freq[2]p2/.p2->(p1^(b2/b1))/.subb/.(OverBar[b])->1/.(Subscript[V, b])->Var/.getme[2.5,Var,1/2])-((OverBar[p])/.(OverBar[p])->freq[1]p1+freq[2]p2/.p2->(p1^(b2/b1))/.subb/.(OverBar[b])->1/.(Subscript[V, b])->Var/.getme2[2.5,Var,1/2])/.R0->5/2/.freq[1]->1/2/.freq[2]->1/2},{Var,0,1,1/100}],Joined->True,PlotStyle->Pink,AxesLabel->{"(Subscript[V, b])","Diff"}]

While not obvious, this has an important implication: one could set α to the level in an unstructured population that would just avoid a second wave, without knowing the underlying heterogeneity in the population, and yet also just avoid a second wave if the population were structured.  The implications for the total number of cases and total number of deaths would, however, be large.

Optimal solution

While having a first wave run through a population can lead to fewer total cases and deaths, the resulting fraction of {p1,p2} is not optimal.  The optimal solution involves immunizing the most influential spreaders, then the next, until 1/R0 is reached.  This is the logic behind first vaccinating health care workers, teachers, and others who have a lot of contact.  What does this optimal strategy look like and how many fewer people would have to be immune?

With two activity classes, this would imply vaccinating the most active class first and not the least active class (p2=1) until the threshold needed to prevent disease spread (((b1^2) p1 freq[1]+(b2^2) p2 freq[2])/((b1^2) freq[1]+(b2^2) freq[2]))==(1/R0) or p1=0 is reached. Then, in the latter case (i.e., when (((b1^2) p1 freq[1]+(b2^2) p2 freq[2])/((b1^2) freq[1]+(b2^2) freq[2]))≥(1/R0) even when p1=0 and p2=1), vaccinating the least active class until the threshold is reached.

Solve[((((b1^2) p1 freq[1]+(b2^2) p2 freq[2])/((b1^2) freq[1]+(b2^2) freq[2]))/.p2->1)==(1/R0),p1]//Flatten

Solve[((((b1^2) p1 freq[1]+(b2^2) p2 freq[2])/((b1^2) freq[1]+(b2^2) freq[2]))/.p1->0)==(1/R0),p2]//Flatten

With equal frequencies of the two classes, only the more active class needs to be vaccinated as long as the variance in mobility ≥ 49-20 (Sqrt[6])~ 0.01:

Solve[((((b2^2) freq[2])/((b1^2) freq[1]+(b2^2) freq[2]))/.subb/.(OverBar[b])->1/.freq[1]->1/2/.freq[2]->1/2)==(1/R0),(Subscript[V, b])];
%/.R0->5/2//Simplify
%//N

An intentional vaccination strategy could thus lead to even fewer cases and deaths (green curve), although the difference isn’t that large if activity is highly variable:

Show[ListPlot[Table[{Var,(OverBar[p])/.(OverBar[p])->freq[1]p1+freq[2]p2/.If[Var≥49-20 (Sqrt[6]),{p1->(((b1^2) freq[1]+(b2^2) freq[2]-(b2^2) R0 freq[2])/((b1^2) R0 freq[1])),p2->1},{p1->0,p2->(((b1^2) freq[1]+(b2^2) freq[2])/((b2^2) R0 freq[2]))}]/.subb/.(OverBar[b])->1/.(Subscript[V, b])->Var/.getme[2.5,Var,1/2]/.R0->5/2/.freq[1]->1/2/.freq[2]->1/2},{Var,0,1,1/200}],Joined->True,PlotStyle->Green,PlotRange->All],
ListPlot[Table[{Var,p1/.If[Var≥49-20 (Sqrt[6]),{p1->(((b1^2) freq[1]+(b2^2) freq[2]-(b2^2) R0 freq[2])/((b1^2) R0 freq[1])),p2->1},{p1->0,p2->(((b1^2) freq[1]+(b2^2) freq[2])/((b2^2) R0 freq[2]))}]/.subb/.(OverBar[b])->1/.(Subscript[V, b])->Var/.getme[2.5,Var,1/2]/.R0->5/2/.freq[1]->1/2/.freq[2]->1/2},{Var,00,1,1/200}],Joined->True,PlotStyle->Black],ListPlot[Table[{Var,p2/.If[Var≥49-20 (Sqrt[6]),{p1->(((b1^2) freq[1]+(b2^2) freq[2]-(b2^2) R0 freq[2])/((b1^2) R0 freq[1])),p2->1},{p1->0,p2->(((b1^2) freq[1]+(b2^2) freq[2])/((b2^2) R0 freq[2]))}]/.subb/.(OverBar[b])->1/.(Subscript[V, b])->Var/.getme[2.5,Var,1/2]/.R0->5/2/.freq[1]->1/2/.freq[2]->1/2},{Var,0,1,1/200}],Joined->True,PlotStyle->{Black,Dashed}],
ListPlot[Table[{Var,(OverBar[p])/.(OverBar[p])->freq[1]p1+freq[2]p2/.p2->(p1^(b2/b1))/.subb/.(OverBar[b])->1/.(Subscript[V, b])->Var/.getme[2.5,Var,1/2]/.R0->5/2/.freq[1]->1/2/.freq[2]->1/2},{Var,0,1,1/100}],Joined->True,PlotStyle->Red],

PlotRange->{{0,1},{0,1}},AxesOrigin->{0,0},AxesLabel->{"(Subscript[V, b])","Freq Susceptibles"}

]

Show[ListPlot[Table[{Var,1-(OverBar[p])/.(OverBar[p])->freq[1]p1+freq[2]p2/.If[Var≥49-20 (Sqrt[6]),{p1->(((b1^2) freq[1]+(b2^2) freq[2]-(b2^2) R0 freq[2])/((b1^2) R0 freq[1])),p2->1},{p1->0,p2->(((b1^2) freq[1]+(b2^2) freq[2])/((b2^2) R0 freq[2]))}]/.subb/.(OverBar[b])->1/.(Subscript[V, b])->Var/.getme[2.5,Var,1/2]/.R0->5/2/.freq[1]->1/2/.freq[2]->1/2},{Var,0,1,1/200}],Joined->True,PlotStyle->Green,PlotRange->All],
ListPlot[Table[{Var,0.6},{Var,00,1,1/200}],Joined->True,PlotStyle->{Dashed,Blue}],ListPlot[Table[{Var,1-(OverBar[p])/.(OverBar[p])->freq[1]p1+freq[2]p2/.p2->(p1^(b2/b1))/.subb/.(OverBar[b])->1/.(Subscript[V, b])->Var/.getme[2.5,Var,1/2]/.R0->5/2/.freq[1]->1/2/.freq[2]->1/2},{Var,0,1,1/100}],Joined->True,PlotStyle->Red],
Graphics[{Text[Style["Blue: everybody equally immune",Blue],{0.74,0.68}],Text[Style["Red: mitigation preventing a second wave",Red],{0.67,0.48}],
Text[Style["Best case immunization program",Green],{0.74,0.25}]}],
PlotRange->{{0,1},{0,1}},AxesOrigin->{0,0},AxesLabel->{"(Subscript[V, b])",None},PlotLabel->"Immunity needed to prevent disease spread"

]

Numerical checks

Here we confirm the main points above. Note that in the plots below, the number of currently infectious individuals is multiplied by five for visibility.

Clear[numsol]
numsol[βI_,κI_,κE_,ϵ_,s10_,s20_,b1_,b2_,α_]:=numsol[βI,κI,κE,ϵ,s10,s20,b1,b2,α]=
NDSolve[{D[s1[t],t]== s1[t] (-(b1^2) α i1[t] βI-b1 b2 α i2[t] βI),
D[e1[t],t]==s1[t] ((b1^2) α i1[t] βI+b1 b2 α i2[t] βI)-e1[t] κE,
D[i1[t],t]==e1[t] κE-i1[t] κI,
D[s2[t],t]==s2[t] (-b1 b2 α i1[t] βI-(b2^2) α i2[t] βI),
D[e2[t],t]==s2[t] (b1 b2 α i1[t] βI+(b2^2) α i2[t] βI)-e2[t] κE,
D[i2[t],t]==e2[t] κE-i2[t] κI,
s1[0]==s10(1-ϵ),e1[0]==s10 ϵ,i1[0]==0,
s2[0]==s20(1-ϵ),e2[0]==s20 ϵ,i2[0]==0},{s1[t],e1[t],i1[t],s2[t],e2[t],i2[t]},{t,0,500}]

NO HETEROGENEITY: Fixing R0=2.5 (βI = 1/2, κI=1/5, b1=b2=1), the first plot shows the dynamics without mitigation, the second with mitigation using α=((Log[(1/R0)])/(1-R0))=0.61, and the third after mitigation is lifted (α=1) using the final fractions of susceptibles after the first wave (restart), confirming that this level of social mitigation prevents a second wave:

restart=Flatten[{s1[t],s2[t]}/.numsol[1/2,1/5,1,1/10000,1/2,1/2,1,1,((Log[(2/5)])/(1-5/2))]/.t->500]

GraphicsGrid[{{
Plot[Evaluate[{s1[t],5i1[t]}/.numsol[1/2,1/5,1,1/10000,1/2,1/2,1,1,1]],{t,0,500},PlotRange->All],
Plot[Evaluate[{s1[t],5i1[t]}/.numsol[1/2,1/5,1,1/10000,1/2,1/2,1,1,((Log[(2/5)])/(1-5/2))]],{t,0,500},PlotRange->All],
Plot[Evaluate[{s1[t],5i1[t]}/.numsol[1/2,1/5,1,1/10000,restart[[1]],restart[[2]],1,1,1]],{t,0,500},PlotRange->All]
}}]

Total number of cases in the unmitigated vs mitigated + relaxed (first and last graph):

{Evaluate[1-(s1[t]+s2[t])/.numsol[1/2,1/5,1,1/10000,1/2,1/2,1,1,1]/.t->500],Evaluate[1-(s1[t]+s2[t])/.numsol[1/2,1/5,1,1/10000,restart[[1]],restart[[2]],1,1,1]/.t->500]}

This illustrates that if α were slightly smaller (0.55) there would be a second wave (last panel):

restart=Flatten[{s1[t],s2[t]}/.numsol[1/2,1/5,1,1/10000,1/2,1/2,1,1,0.55]/.t->500];
GraphicsGrid[{{
Plot[Evaluate[{s1[t],5i1[t]}/.numsol[1/2,1/5,1,1/10000,1/2,1/2,1,1,1]],{t,0,500},PlotRange->All],
Plot[Evaluate[{s1[t],5i1[t]}/.numsol[1/2,1/5,1,1/10000,1/2,1/2,1,1,0.55]],{t,0,500},PlotRange->All],
Plot[Evaluate[{s1[t],5i1[t]}/.numsol[1/2,1/5,1,1/10000,restart[[1]],restart[[2]],1,1,1]],{t,0,500},PlotRange->All]
}}]

Total number of cases in the unmitigated vs mitigated + relaxed (first and last graph):

{Evaluate[1-(s1[t]+s2[t])/.numsol[1/2,1/5,1,1/10000,1/2,1/2,1,1,1]/.t->500],Evaluate[1-(s1[t]+s2[t])/.numsol[1/2,1/5,1,1/10000,restart[[1]],restart[[2]],1,1,1]/.t->500]}

WITH HETEROGENEITY: The same plots but with b1=1.5 and b2=0.5 ((OverBar[b])=1, (Subscript[V, b])=1/4, βI now has to be 2 κI), the following confirms that α=((Log[(1/R0)])/(1-R0)) still prevents the second wave.

((R0 κI)/((b1^2) freq[1]+(b2^2) freq[2]))/.b1->3/2/.b2->1/2/.freq[1]->1/2/.freq[2]->1/2/.R0->5/2

restart=Flatten[{s1[t],s2[t]}/.numsol[2/5,1/5,1,1/10000,1/2,1/2,1.5,0.5,((Log[(2/5)])/(1-5/2))]/.t->500];
GraphicsGrid[{{
Plot[Evaluate[{s1[t],5 i1[t]}/.numsol[2/5,1/5,1,1/10000,1/2,1/2,1.5,0.5,1]],{t,0,500},PlotRange->All],
Plot[Evaluate[{s1[t],5 i1[t]}/.numsol[2/5,1/5,1,1/10000,1/2,1/2,1.5,0.5,((Log[(2/5)])/(1-5/2))]],{t,0,500},PlotRange->All],
Plot[Evaluate[{s1[t],5i1[t]}/.numsol[2/5,1/5,1,1/10000,restart[[1]],restart[[2]],1.5,0.5,1]],{t,0,500},PlotRange->All]
}}]

Total number of cases in the unmitigated vs mitigated + relaxed (first and last graph):

{Evaluate[1-(s1[t]+s2[t])/.numsol[2/5,1/5,1,1/10000,1/2,1/2,1.5,0.5,1]/.t->500],Evaluate[1-(s1[t]+s2[t])/.numsol[2/5,1/5,1,1/10000,restart[[1]],restart[[2]],1.5,0.5,1]/.t->500]}

(Note that heterogeneity reduces the total number of cases even without mitigation, because the individuals remaining as Rt drops are more likely the least active.)

This plot illustrates that if α were slightly smaller (0.55) there would be a second wave (last panel):

restart=Flatten[{s1[t],s2[t]}/.numsol[2/5,1/5,1,1/10000,1/2,1/2,1.5,0.5,0.55]/.t->500];
GraphicsGrid[{{
Plot[Evaluate[{s1[t],5 i1[t]}/.numsol[2/5,1/5,1,1/10000,1/2,1/2,1.5,0.5,1]],{t,0,500},PlotRange->All],
Plot[Evaluate[{s1[t],5 i1[t]}/.numsol[2/5,1/5,1,1/10000,1/2,1/2,1.5,0.5,0.55]],{t,0,500},PlotRange->All],
Plot[Evaluate[{s1[t],5i1[t]}/.numsol[2/5,1/5,1,1/10000,restart[[1]],restart[[2]],1.5,0.5,1]],{t,0,500},PlotRange->All]
}}]

Total number of cases in the unmitigated vs mitigated + relaxed (first and last graph):

{Evaluate[1-(s1[t]+s2[t])/.numsol[2/5,1/5,1,1/10000,1/2,1/2,1.5,0.5,1]/.t->500],Evaluate[1-(s1[t]+s2[t])/.numsol[2/5,1/5,1,1/10000,restart[[1]],restart[[2]],1.5,0.5,1]/.t->500]}