Exploring transients around steady state: (Subscript[SIR, n]) model

Model includes n waning classes, corresponding to reducing antibody levels over time.

Variant evolution

Figure 2B: persistently immune evasive (n=5) (increases mean waning rate by 67%)
Can infect last two classes (m=2)

maxtime=500;
maxy=8;

Parameters:

TRYp=TRYp; (*No heterogeneity, so p not relevant now.*)
TRYf=99/100;(*Starting fraction of resident*)
TRYolder=13/100; (*~ Age 70+ from StatsCan 5022509/38929902 *)
	(*https://www150.statcan.gc.ca/t1/tbl1/en/cv.action?pid=1710000501*)
TRYκ=1/5; (*Five day infectious period*)
TRYδ=1/125; (*Waning rate based on Menegale et al. (2023)*) 
TRYδO=1/125; (*No clear age difference noted in Menegale et al. (2023)*) 
TRYconvert=10/10; (*Only this fraction of cases seroconvert or boost immunity*)

i->((δ (β-κ))/(β (δ+q κ)))

Solve[((δ (β-κ))/(β (δ+q κ)))==inf,β]

TRYinf=1/50; (*Assumed initial endemic frequency of infected individuals*)
TRYβ=((TRYδ TRYκ)/(TRYδ-TRYinf TRYδ-TRYconvert TRYinf TRYκ));

((TRYβ)/(TRYκ))//N

Variant has the same transmissibility:

TRYβV=TRYβ;

but infects earlier during waning.  Specifically, we assume five waning classes and the variant infects two earlier.

TRYn=5;
TRYm=2;

Thus, the mean waning rate (once fixed) increases by 67%:

((TRYδ+Δδ)/(TRYδ))/.Δδ->(TRYm/(TRYn-TRYm))TRYδ//N

while the mean waiting time decreases by 40%:

1/%

Starting equilibrium would be (r0 includes all of the n resistance classes):
{s->(κ/β),i->((δ (β-κ))/(β (δ+q κ)))}

start={s0->((TRYκ)/(TRYβ)),i0->((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))),r0->1-((TRYκ)/(TRYβ))-((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert)))}

Waning distribution for a population of younger and older individuals:

Show[
Plot[TRYconvert*(1-CDF[GammaDistribution[TRYn,1/(TRYδ TRYn)],t]),{t,0,365},PlotRange->{All,{0,1}},PlotStyle->RGBColor[0.4,0.4,0.8]],
Plot[TRYconvert*(1-CDF[GammaDistribution[TRYn,1/(TRYδO TRYn)],t]),{t,0,365},PlotRange->{All,{0,1}},PlotStyle->RGBColor[0.4,0.6,0.4]]
]

Now the susceptible and resistant classes depend only on age (not whether they carried a variant before):
[NOTE: For the ease of colouring, the first stage class is the variant.]

{ages=2,stages=2,n=TRYn};

Clear[solution]
solution[{p_,f_,older_,β_,βV_,δY_,δO_,κ_,q_},start]:=solution[{p,f,older,β,βV,δY,δO,κ,q},start]=
Block[{ages=2,stages=2,n=TRYn},
dsdt[j_,t_]=(1-q)*Sum[κ[j,jj]*i[j,jj,t],{jj,1,stages}]+n*δ[j]*r[n,j,t]-Sum[β[k,kk]*i[k,kk,t],{k,1,ages},{kk,1,stages}]*s[j,t] ;didt[j_,1,t_]=Sum[β[k,1]*i[k,1,t],{k,1,ages}]*Sum[r[nn,j,t],{nn,n+1-TRYm,n}]+Sum[β[k,1]*i[k,1,t],{k,1,ages}]*s[j,t] -κ[j,1]*i[j,1,t];didt[j_,2,t_]=Sum[β[k,2]*i[k,2,t],{k,1,ages}]*s[j,t] -κ[j,2]*i[j,2,t];drdt[1,j_,t_]=q*Sum[κ[j,jj]*i[j,jj,t],{jj,1,stages}]-n*δ[j]*r[1,j,t];For[nn=2,nn≤n-TRYm,nn++,
drdt[nn,j_,t_]=n*δ[j]*r[nn-1,j,t]-n*δ[j]*r[nn,j,t];
];
For[nn=n+1-TRYm,nn≤n,nn++,
drdt[nn,j_,t_]=n*δ[j]*r[nn-1,j,t]-Sum[β[k,1]*i[k,1,t],{k,1,ages}]*r[nn,j,t]-n*δ[j]*r[nn,j,t];
];
pars={β[j_,1]->βV,β[j_,2]->β,δ[1]->δY,δ[2]->δO,κ[j_,jj_]->κ}; 
frac[1,1]=(1-older)(1-f);
frac[2,1]=older(1-f);
frac[1,2]=(1-older)f;
frac[2,2]=older f;
nvars=Drop[Flatten[Table[{s[j,t],Table[i[j,jj,t],{jj,1,stages}],Table[r[nn,j,t],{nn,1,n}]},{j,1,ages}]],-1];
r[n,ages,t_]=1-Total[nvars];
neqns=Drop[Flatten[Table[{D[s[j,t],t]==dsdt[j,t],Table[D[i[j,jj,t],t]==didt[j,jj,t],{jj,1,stages}],Table[D[r[nn,j,t],t]==drdt[nn,j,t],{nn,1,n}]},{j,1,ages}]],-1];
nstart=Drop[Flatten[Table[{
Table[i[j,jj,0]==frac[j,jj] i0,{jj,1,stages}],
s[j,0]==(frac[j,1]+frac[j,2])*s0,Table[r[nn,j,0]==(frac[j,1]+frac[j,2]) (r0/n),{nn,1,n}]},{j,1,ages}]]/.start,-1];
NDSolve[Flatten[{nstart/.pars,neqns/.pars}],nvars,{t,0,maxtime}]
]

all=Flatten[Table[{j,jj},{j,1,ages},{jj,1,stages}],1]

labels={"<70,Res","<70,Mutant","70+,Res","70+,Mutant"};

colours={RGBColor[0.4,0.4,0.8],ColorData["Pastel",1],RGBColor[0.4,0.6,0.4], ColorData["Pastel",3/4]}

coltab=Join[{{1->{0,colours[[1]]}}},Table[{i->{{i-1},colours[[i]]}},{i,2,Length[all]}]];
Table[{labels[[i]],colours[[i]]},{i,1,Length[all]}]//MatrixForm

partemp={TRYp,TRYf,TRYolder,TRYβ,TRYβV,TRYδ,TRYδO,TRYκ,TRYconvert};

solution[partemp,start];

Showing the % of the population infectious:

plot1 =Show[
Plot[Evaluate@Table[ Sum[100*i[all[[c,1]],all[[c,2]],t]/.solution[partemp,start],{c,1,jcum}],{jcum,1,Length[all]}],{t,0,maxtime},PlotRange->{Automatic,{0,maxy}},
PlotStyle->None,Filling->{{1->{0,}},{2->{{1},}},{3->{{2},}},{4->{{3},}}}],
Plot[100((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black,Dashed}],
Plot[100(((TRYn/(TRYn-TRYm))TRYδ (TRYβV-TRYκ))/(TRYβV ((TRYn/(TRYn-TRYm))TRYδ+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black}],
ImageSize->250
]

Mutant fraction change between t=0 and t=50:

((Sum[i[j,1,t],{j,1,ages}])/(Sum[i[j,jj,t],{j,1,ages},{jj,1,stages}]))/.solution[partemp,start]/.t->50

Solve[Exp[s 50] ==( ((p50/(1-p50))/(p0/(1-p0)))/.p50->%/.p0->1-TRYf),s][[1]]

This is consistent with the predicted selection coefficient:

((TRYβV-TRYβ)/(TRYβ))TRYκ+(TRYm/TRYn)(((TRYκ TRYconvert(TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))))TRYβV//N

New equilibrium:

((((TRYδ+Δδ)(TRYβV-TRYκ))/(TRYβV ((TRYδ+Δδ)+TRYκ TRYconvert))))/.Δδ->(TRYm/(TRYn-TRYm))TRYδ//N

% rise in the endemic equilibrium:

100*(((%-TRYinf))/TRYinf)

The following illustrates the relatively minor effect caused by reducing seroconversion (to q = 80%), as long as δ/q is held constant:

partemp={TRYp,TRYf,TRYolder,TRYβ,TRYβV,(8/10)TRYδ,(8/10)TRYδO,TRYκ,(8/10)};

solution[partemp,start];

Showing the % of the population infectious:

plot2=Show[plot1,
Plot[Evaluate@Table[ Sum[100*i[all[[c,1]],all[[c,2]],t]/.solution[partemp,start],{c,1,jcum}],{jcum,1,Length[all]}],{t,0,maxtime},PlotRange->{Automatic,{0,maxy}},
PlotStyle->{,,,}],
Plot[100((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black,Dashed}],
Plot[100(((TRYn/(TRYn-TRYm))TRYδ (TRYβV-TRYκ))/(TRYβV ((TRYn/(TRYn-TRYm))TRYδ+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black}],
ImageSize->250
]

Seroconversion affects the dynamics (curves compared to the previous shaded plot) but not the initial speed of spread of the variant (described by the selection coefficient) or the equilibrium.

Clear[solution, i, s, r, ages, stages, n, neqns, nstart, nvars, didt, dsdt, drdt];

Figure 2A: transiently immune evasive (n=5) (increases mean waning rate by 67%)
Can infect last two classes (m=2)

maxtime=500;
maxy=8;

Parameters:

TRYp=TRYp; (*No heterogeneity, so p not relevant now.*)
TRYf=99/100;(*Starting fraction of resident*)
TRYolder=13/100; (*~ Age 70+ from StatsCan 5022509/38929902 *)
	(*https://www150.statcan.gc.ca/t1/tbl1/en/cv.action?pid=1710000501*)
TRYκ=1/5; (*Five day infectious period*)
TRYδ=1/125; (*Waning rate based on Menegale et al. (2023)*) 
TRYδO=1/125; (*No clear age difference noted in Menegale et al. (2023)*) 
TRYconvert=10/10; (*Only this fraction of cases seroconvert or boost immunity*)

i->((δ (β-κ))/(β (δ+q κ)))

Solve[((δ (β-κ))/(β (δ+q κ)))==inf,β]

TRYinf=1/50; (*Assumed initial endemic frequency of infected individuals*)
TRYβ=((TRYδ TRYκ)/(TRYδ-TRYinf TRYδ-TRYconvert TRYinf TRYκ));

((TRYβ)/(TRYκ))//N

Variant has the same transmissibility:

TRYβV=TRYβ;

but infects earlier during waning.  Specifically, we assume five waning classes and the variant infects two earlier (40% reduction in mean waning time):

TRYn=5;
TRYm=2;

Starting equilibrium would be (r0 includes all of the n resistance classes):
{s->(κ/β),i->((δ (β-κ))/(β (δ+q κ)))}

start={s0->((TRYκ)/(TRYβ)),i0->((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))),r0->1-((TRYκ)/(TRYβ))-((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert)))}

Waning distribution for a population of younger and older individuals:

Show[
Plot[TRYconvert*(1-CDF[GammaDistribution[TRYn,1/(TRYδ TRYn)],t]),{t,0,365},PlotRange->{All,{0,1}},PlotStyle->RGBColor[0.4,0.4,0.8]],
Plot[TRYconvert*(1-CDF[GammaDistribution[TRYn,1/(TRYδO TRYn)],t]),{t,0,365},PlotRange->{All,{0,1}},PlotStyle->RGBColor[0.4,0.6,0.4]]
]

Now the susceptible and resistant classes depend only on age (not whether they carried a variant before):
[NOTE: For the ease of colouring, the first stage class is the variant.]

{ages=2,stages=2,exposures=2,n=TRYn};

Clear[solution]
solution[{p_,f_,older_,β_,βV_,δY_,δO_,κ_,q_},start]:=solution[{p,f,older,β,βV,δY,δO,κ,q},start]=
Block[{ages=2,stages=2,exposures=2,n=TRYn},
dsdt[j_,zz_,t_]=(1-q)*Sum[κ[j,jj]*i[j,jj,zz,t],{jj,1,stages}]+n*δ[j]*r[n,j,zz,t]-Sum[β[k,kk]*i[k,kk,kkk,t],{k,1,ages},{kk,1,stages},{kkk,1,exposures}]*s[j,zz,t] ;
(*For the susceptible and recovered classes, we introduce zz to measure whether individuals have ever been exposed to the new variant (zz=1) or not (zz=2), this allows us to model transient immunity. We assume that those who do not seroconvert, with probability 1-c, remain in their zz class.*)didt[j_,1,1,t_]=Sum[β[k,1]*i[k,1,kkk,t],{k,1,ages},{kkk,1,exposures}]*Sum[r[nn,j,2,t],{nn,n+1-TRYm,n}]+Sum[β[k,1]*i[k,1,kkk,t],{k,1,ages},{kkk,1,exposures}]*s[j,1,t] -κ[j,1]*i[j,1,1,t];
(*This is the variant [stage 1], but it can only infect the recovering individuals who have never been exposed to a variant infection [stage 2]. Once infected by the variant these will move into the first type of recovering individual [zz=1] and are no longer available*)
didt[j_,1,2,t_]=Sum[β[k,1]*i[k,1,kkk,t],{k,1,ages},{kkk,1,exposures}]*s[j,2,t] -κ[j,1]*i[j,1,2,t];(*Although no infection with the variant lacks exposure to the variant, we allow this in case some fail to seroconvert, with probability 1-c.*)
didt[j_,2,zz_,t_]=Sum[β[k,2]*i[k,2,kkk,t],{k,1,ages},{kkk,1,exposures}]*s[j,zz,t]-κ[j,2]*i[j,2,zz,t];drdt[1,j_,1,t_]=q*(κ[j,1]*i[j,1,1,t]+κ[j,1]*i[j,1,2,t]+κ[j,2]*i[j,2,1,t])-n*δ[j]*r[1,j,1,t];
(*All infections with the variant move into this first resistance class only, as do those exposed in prior infections (i[j,2,1,t]).*)
drdt[1,j_,2,t_]=q*(κ[j,2]*i[j,2,2,t])-n*δ[j]*r[1,j,2,t];
For[nn=2,nn≤n-TRYm,nn++,
drdt[nn,j_,zz_,t_]=n*δ[j]*r[nn-1,j,zz,t]-n*δ[j]*r[nn,j,zz,t];
];
For[nn=n+1-TRYm,nn≤n,nn++,
drdt[nn,j_,1,t_]=n*δ[j]*r[nn-1,j,1,t]-n*δ[j]*r[nn,j,1,t];
drdt[nn,j_,2,t_]=n*δ[j]*r[nn-1,j,2,t]-Sum[β[k,1]*i[k,1,kkk,t],{k,1,ages},{kkk,1,exposures}]*r[nn,j,2,t]-n*δ[j]*r[nn,j,2,t];
(*Only those that have not been infected before by the new variant [stage 2] can become infected earlier by the new variant [stage 1]*)];
pars={β[j_,1]->βV,β[j_,2]->β,δ[1]->δY,δ[2]->δO,κ[j_,jj_]->κ}; 
frac[1,1]=(1-older)(1-f);
frac[2,1]=older(1-f);
frac[1,2]=(1-older)f;
frac[2,2]=older f;
nvars=Drop[Flatten[Table[{Table[s[j,zz,t],{zz,1,exposures}],
Table[i[j,jj,zz,t],{jj,1,stages},{zz,1,exposures}],
Table[r[nn,j,zz,t],{nn,1,n},{zz,1,exposures}]},{j,1,ages}]],-1];
r[n,ages,exposures,t_]=1-Total[nvars];
neqns=Drop[Flatten[Table[{Table[D[s[j,zz,t],t]==dsdt[j,zz,t],{zz,1,exposures}],Table[D[i[j,jj,zz,t],t]==didt[j,jj,zz,t],{jj,1,stages},{zz,1,exposures}],Table[D[r[nn,j,zz,t],t]==drdt[nn,j,zz,t],{nn,1,n},{zz,1,exposures}]},{j,1,ages}]],-1];
nstart=Drop[Flatten[Table[{
Table[i[j,jj,1,0]==0,{jj,1,stages}],
Table[i[j,jj,2,0]==frac[j,jj] i0,{jj,1,stages}],
s[j,1,0]==0,(*Nobody starts exposed to the variant*)
s[j,2,0]==(frac[j,1]+frac[j,2])*s0,
Table[r[nn,j,1,0]==0,{nn,1,n}],
Table[r[nn,j,2,0]==(frac[j,1]+frac[j,2]) (r0/n),{nn,1,n}]},{j,1,ages}]]/.start,-1];(*All susceptible and resistant individuals initially susceptible to the new variant [in stage 2] *)
NDSolve[Flatten[{nstart/.pars,neqns/.pars}],nvars,{t,0,maxtime}]
]

all=Flatten[Table[{j,jj},{j,1,ages},{jj,1,stages}],1]

labels={"<70,Res","<70,Mutant","70+,Res","70+,Mutant"};

colours={RGBColor[0.4,0.4,0.8],ColorData["Pastel",1],RGBColor[0.4,0.6,0.4], ColorData["Pastel",3/4]}

coltab=Join[{{1->{0,colours[[1]]}}},Table[{i->{{i-1},colours[[i]]}},{i,2,Length[all]}]];
Table[{labels[[i]],colours[[i]]},{i,1,Length[all]}]//MatrixForm

partemp={TRYp,TRYf,TRYolder,TRYβ,TRYβV,TRYδ,TRYδO,TRYκ,TRYconvert};

solution[partemp,start];

Showing the % of the population infectious:

plot1=Show[
Plot[Evaluate@Table[ Sum[100i[all[[c,1]],all[[c,2]],1,t]+100i[all[[c,1]],all[[c,2]],2,t]/.solution[partemp,start],{c,1,jcum}],{jcum,1,Length[all]}],{t,0,maxtime},PlotRange->{Automatic,{0,maxy}},
PlotStyle->None,Filling->{{1->{0,}},{2->{{1},}},{3->{{2},}},{4->{{3},}}}],
Plot[100((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black,Dashed}],
Plot[100(((TRYn/(TRYn-0))TRYδ (TRYβV-TRYκ))/(TRYβV ((TRYn/(TRYn-0))TRYδ+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black}],
ImageSize->250
]

Mutant fraction change between t=0 and t=50:

((Sum[i[j,1,zz,t],{j,1,ages},{zz,1,exposures}])/(Sum[i[j,jj,zz,t],{j,1,ages},{jj,1,stages},{zz,1,exposures}]))/.solution[partemp,start]/.t->50

Solve[Exp[s 50] ==( ((p50/(1-p50))/(p0/(1-p0)))/.p50->%/.p0->1-TRYf),s][[1]]

This is consistent with the predicted selection coefficient:

((TRYβV-TRYβ)/(TRYβ))TRYκ+(TRYm/TRYn)(((TRYκ TRYconvert(TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))))TRYβV//N

New equilibrium:

((((TRYδ+Δδ)(TRYβV-TRYκ))/(TRYβV ((TRYδ+Δδ)+TRYκ TRYconvert))))/.Δδ->(0/(TRYn-0))TRYδ//N

% rise in the endemic equilibrium:

100*(((%-TRYinf))/TRYinf)

The following illustrates the relatively minor effect caused by reducing seroconversion (to q = 80%), as long as δ/q is held constant:

partemp={TRYp,TRYf,TRYolder,TRYβ,TRYβV,(8/10)TRYδ,(8/10)TRYδO,TRYκ,(8/10)};

solution[partemp,start];

Showing the % of the population infectious:

plot2=Show[plot1,
Plot[Evaluate@Table[ Sum[100*i[all[[c,1]],all[[c,2]],t]/.solution[partemp,start],{c,1,jcum}],{jcum,1,Length[all]}],{t,0,maxtime},PlotRange->{Automatic,{0,maxy}},
PlotStyle->{,,,}],
Plot[100((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black,Dashed}],
Plot[100(((TRYn/(TRYn-0))TRYδ (TRYβV-TRYκ))/(TRYβV ((TRYn/(TRYn-0))TRYδ+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black}],
ImageSize->250
]

Seroconversion affects the dynamics (curves compared to the previous shaded plot) but not the initial speed of spread of the variant (described by the selection coefficient) or the equilibrium.

Clear[solution, i, s, r, ages, stages, n, neqns, nstart, nvars, didt, dsdt, drdt];

Figure 3 & S2B: more transmissible (n=5)

maxtime=500;
maxy=8;

Parameters:

TRYp=TRYp; (*No heterogeneity, so p not relevant now.*)
TRYf=99/100;(*Starting fraction of resident*)
TRYolder=13/100; (*~ Age 70+ from StatsCan 5022509/38929902 *)
	(*https://www150.statcan.gc.ca/t1/tbl1/en/cv.action?pid=1710000501*)
TRYκ=1/5; (*Five day infectious period*)
TRYδ=1/125; (*Waning rate based on Menegale et al. (2023)*) 
TRYδO=1/125; (*No clear age difference noted in Menegale et al. (2023)*) 
TRYconvert=10/10; (*Only this fraction of cases seroconvert or boost immunity*)

i->((δ (β-κ))/(β (δ+q κ)))

Solve[((δ (β-κ))/(β (δ+q κ)))==inf,β]

TRYinf=1/50; (*Assumed initial endemic frequency of infected individuals*)
TRYβ=((TRYδ TRYκ)/(TRYδ-TRYinf TRYδ-TRYconvert TRYinf TRYκ));

((TRYβ)/(TRYκ))//N

Variant has a higher transmissibility (set to have the same selection coefficient as for the immune evasive variant):

Flatten[Solve[((βV-TRYβ)/(TRYβ))TRYκ==0.0833,βV]]

TRYβV=βV/.%;

% increase in β:

100*((TRYβV-TRYβ)/(TRYβ))

and now has no immune evasive properties.  Specifically, we assume five waning classes and the variant infects only susceptibles (m=0):

TRYn=5;
TRYm=0;

Starting equilibrium would be (r0 includes all of the n resistance classes):
{s->(κ/β),i->((δ (β-κ))/(β (δ+q κ)))}

start={s0->((TRYκ)/(TRYβ)),i0->((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))),r0->1-((TRYκ)/(TRYβ))-((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert)))}

Waning distribution for a population of younger and older individuals:

Show[
Plot[TRYconvert*(1-CDF[GammaDistribution[TRYn,1/(TRYδ TRYn)],t]),{t,0,365},PlotRange->{All,{0,1}},PlotStyle->RGBColor[0.4,0.4,0.8]],
Plot[TRYconvert*(1-CDF[GammaDistribution[TRYn,1/(TRYδO TRYn)],t]),{t,0,365},PlotRange->{All,{0,1}},PlotStyle->RGBColor[0.4,0.6,0.4]]
]

Now the susceptible and resistant classes depend only on age (not whether they carried a variant before):
[NOTE: For the ease of colouring, the first stage class is the variant.]

{ages=2,stages=2,n=TRYn};

Clear[solution]
solution[{p_,f_,older_,β_,βV_,δY_,δO_,κ_,q_},start]:=solution[{p,f,older,β,βV,δY,δO,κ,q},start]=
Block[{ages=2,stages=2,n=TRYn},
dsdt[j_,t_]=(1-q)*Sum[κ[j,jj]*i[j,jj,t],{jj,1,stages}]+n*δ[j]*r[n,j,t]-Sum[β[k,kk]*i[k,kk,t],{k,1,ages},{kk,1,stages}]*s[j,t] ;didt[j_,1,t_]=Sum[β[k,1]*i[k,1,t],{k,1,ages}]*Sum[r[nn,j,t],{nn,n+1-TRYm,n}]+Sum[β[k,1]*i[k,1,t],{k,1,ages}]*s[j,t] -κ[j,1]*i[j,1,t];didt[j_,2,t_]=Sum[β[k,2]*i[k,2,t],{k,1,ages}]*s[j,t] -κ[j,2]*i[j,2,t];drdt[1,j_,t_]=q*Sum[κ[j,jj]*i[j,jj,t],{jj,1,stages}]-n*δ[j]*r[1,j,t];For[nn=2,nn≤n-TRYm,nn++,
drdt[nn,j_,t_]=n*δ[j]*r[nn-1,j,t]-n*δ[j]*r[nn,j,t];
];
For[nn=n+1-TRYm,nn≤n,nn++,
drdt[nn,j_,t_]=n*δ[j]*r[nn-1,j,t]-Sum[β[k,1]*i[k,1,t],{k,1,ages}]*r[nn,j,t]-n*δ[j]*r[nn,j,t];
];
pars={β[j_,1]->βV,β[j_,2]->β,δ[1]->δY,δ[2]->δO,κ[j_,jj_]->κ}; 
frac[1,1]=(1-older)(1-f);
frac[2,1]=older(1-f);
frac[1,2]=(1-older)f;
frac[2,2]=older f;
nvars=Drop[Flatten[Table[{s[j,t],Table[i[j,jj,t],{jj,1,stages}],Table[r[nn,j,t],{nn,1,n}]},{j,1,ages}]],-1];
r[n,ages,t_]=1-Total[nvars];
neqns=Drop[Flatten[Table[{D[s[j,t],t]==dsdt[j,t],Table[D[i[j,jj,t],t]==didt[j,jj,t],{jj,1,stages}],Table[D[r[nn,j,t],t]==drdt[nn,j,t],{nn,1,n}]},{j,1,ages}]],-1];
nstart=Drop[Flatten[Table[{
Table[i[j,jj,0]==frac[j,jj] i0,{jj,1,stages}],
s[j,0]==(frac[j,1]+frac[j,2])*s0,Table[r[nn,j,0]==(frac[j,1]+frac[j,2]) (r0/n),{nn,1,n}]},{j,1,ages}]]/.start,-1];
NDSolve[Flatten[{nstart/.pars,neqns/.pars}],nvars,{t,0,maxtime}]
]

all=Flatten[Table[{j,jj},{j,1,ages},{jj,1,stages}],1]

labels={"<70,Res","<70,Mutant","70+,Res","70+,Mutant"};

colours={RGBColor[0.4,0.4,0.8],ColorData["Pastel",1],RGBColor[0.4,0.6,0.4], ColorData["Pastel",3/4]}

coltab=Join[{{1->{0,colours[[1]]}}},Table[{i->{{i-1},colours[[i]]}},{i,2,Length[all]}]];
Table[{labels[[i]],colours[[i]]},{i,1,Length[all]}]//MatrixForm

partemp={TRYp,TRYf,TRYolder,TRYβ,TRYβV,TRYδ,TRYδO,TRYκ,TRYconvert};

solution[partemp,start];

Showing the % of the population infectious:

plot1=Show[
Plot[Evaluate@Table[ Sum[100 i[all[[c,1]],all[[c,2]],t]/.solution[partemp,start],{c,1,jcum}],{jcum,1,Length[all]}],{t,0,maxtime},PlotRange->{Automatic,{0,maxy}},
PlotStyle->None,Filling->{{1->{0,}},{2->{{1},}},{3->{{2},}},{4->{{3},}}}],
Plot[100((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black,Dashed}],
Plot[100(((TRYn/(TRYn-TRYm))TRYδ (TRYβV-TRYκ))/(TRYβV ((TRYn/(TRYn-TRYm))TRYδ+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black}],
ImageSize->250
]

Mutant fraction change between t=0 and t=50:

((Sum[i[j,1,t],{j,1,ages}])/(Sum[i[j,jj,t],{j,1,ages},{jj,1,stages}]))/.solution[partemp,start]/.t->50

Solve[Exp[s 50] ==( ((p50/(1-p50))/(p0/(1-p0)))/.p50->%/.p0->1-TRYf),s][[1]]

This is consistent with the predicted selection coefficient:

((TRYβV-TRYβ)/(TRYβ))TRYκ+(TRYm/TRYn)(((TRYκ TRYconvert(TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))))TRYβV//N

New equilibrium:

((((TRYδ+Δδ)(TRYβV-TRYκ))/(TRYβV ((TRYδ+Δδ)+TRYκ TRYconvert))))/.Δδ->(TRYm/(TRYn-TRYm))TRYδ//N

% rise in the endemic equilibrium:

100*(((%-TRYinf))/TRYinf)

This considers the number of susceptibles, with the black line at 1/(Subscript[(OverTilde[R]), 0]) (the fraction of susceptibles that causes the disease to stop growing before the variant) and the dashed black line at 1/(Subscript[(OverTilde[R]), 0]^( * )) (the fraction of susceptibles that causes the disease to stop growing once the variant predominates):

Show[
Plot[Evaluate@Table[ Sum[100 s[c,t]/.solution[partemp,start],{c,1,jcum}],{jcum,1,2}],{t,0,maxtime},PlotRange->{Automatic,{0,100}},
PlotStyle->None,Filling->{{1->{0,}},{2->{{1},}}}],
Plot[100((TRYκ)/(TRYβ)),{t,0,maxtime},PlotStyle->{Black}],
Plot[100((TRYκ)/(TRYβV)),{t,0,maxtime},PlotStyle->{Black,Dashed}],
ImageSize->250
]

The first time that the number of susceptibles falls below ((TRYκ)/(TRYβV))is when we expect the number of the new variant to start declining (roughly, given that waning changes the dynamics):

(1+(((TRYβV-TRYβ))/(TRYβ)))

droptime=First[Select[Table[Flatten[{t, Sum[100 s[c,t]/.solution[partemp,start],{c,1,2}]}],{t,0,maxtime}],#[[2]]<100((TRYκ)/(TRYβV))&]]

Show[
plot1,
ListPlot[{{droptime[[1]],0},{droptime[[1]],maxy}},Joined->True,PlotStyle->{Gray}],
ImageSize->250
]

The following illustrates the relatively minor effect caused by reducing seroconversion (to q = 80%), as long as δ/q is held constant:

partemp={TRYp,TRYf,TRYolder,TRYβ,TRYβV,(8/10)TRYδ,(8/10)TRYδO,TRYκ,(8/10)};

solution[partemp,start];

Showing the % of the population infectious:

plot2=Show[plot1,
Plot[Evaluate@Table[ Sum[100*i[all[[c,1]],all[[c,2]],t]/.solution[partemp,start],{c,1,jcum}],{jcum,1,Length[all]}],{t,0,maxtime},PlotRange->{Automatic,{0,maxy}},
PlotStyle->{,,,}],
Plot[100((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black,Dashed}],
Plot[100(((TRYn/(TRYn-0))TRYδ (TRYβV-TRYκ))/(TRYβV ((TRYn/(TRYn-0))TRYδ+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black}],
ImageSize->250
]

Seroconversion affects the dynamics (curves compared to the previous shaded plot) but not the initial speed of spread of the variant (described by the selection coefficient) or the equilibrium.

Clear[solution, i, s, r, ages, stages, n, neqns, nstart, nvars, didt, dsdt, drdt];

Figure S1A: both persistently immune evasive and more transmissible (n=5)
Can infect last class (m=1)

maxtime=500;
maxy=8;

Parameters:

TRYp=TRYp; (*No heterogeneity, so p not relevant now.*)
TRYf=99/100;(*Starting fraction of resident*)
TRYolder=13/100; (*~ Age 70+ from StatsCan 5022509/38929902 *)
	(*https://www150.statcan.gc.ca/t1/tbl1/en/cv.action?pid=1710000501*)
TRYκ=1/5; (*Five day infectious period*)
TRYδ=1/125; (*Waning rate based on Menegale et al. (2023)*) 
TRYδO=1/125; (*No clear age difference noted in Menegale et al. (2023)*) 
TRYconvert=10/10; (*Only this fraction of cases seroconvert or boost immunity*)

i->((δ (β-κ))/(β (δ+q κ)))

Solve[((δ (β-κ))/(β (δ+q κ)))==inf,β]

TRYinf=1/50; (*Assumed initial endemic frequency of infected individuals*)
TRYβ=((TRYδ TRYκ)/(TRYδ-TRYinf TRYδ-TRYconvert TRYinf TRYκ));

((TRYβ)/(TRYκ))//N

Now we allow partial immune evasiveness, coupled with higher transmission.  Specifically, we assume five waning classes and the variant infects one earlier resistant class (m=1):

TRYn=5;
TRYm=1;

Variant has a higher transmissibility (set to have the same selection coefficient as for the immune evasive variant):

Flatten[Solve[((βV-TRYβ)/(TRYβ))TRYκ+(TRYm/TRYn)(((TRYκ TRYconvert(TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))))βV==0.0833,βV]]

TRYβV=βV/.%;

TRYβ//N

TRYβV/TRYβ//N

Starting equilibrium would be (r0 includes all of the n resistance classes):
{s->(κ/β),i->((δ (β-κ))/(β (δ+q κ)))}

start={s0->((TRYκ)/(TRYβ)),i0->((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))),r0->1-((TRYκ)/(TRYβ))-((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert)))}

Waning distribution for a population of younger and older individuals:

Show[
Plot[TRYconvert*(1-CDF[GammaDistribution[TRYn,1/(TRYδ TRYn)],t]),{t,0,365},PlotRange->{All,{0,1}},PlotStyle->RGBColor[0.4,0.4,0.8]],
Plot[TRYconvert*(1-CDF[GammaDistribution[TRYn,1/(TRYδO TRYn)],t]),{t,0,365},PlotRange->{All,{0,1}},PlotStyle->RGBColor[0.4,0.6,0.4]]
]

Now the susceptible and resistant classes depend only on age (not whether they carried a variant before):
[NOTE: For the ease of colouring, the first stage class is the variant.]

{ages=2,stages=2,n=TRYn};

Clear[solution]
solution[{p_,f_,older_,β_,βV_,δY_,δO_,κ_,q_},start]:=solution[{p,f,older,β,βV,δY,δO,κ,q},start]=
Block[{ages=2,stages=2,n=TRYn},
dsdt[j_,t_]=(1-q)*Sum[κ[j,jj]*i[j,jj,t],{jj,1,stages}]+n*δ[j]*r[n,j,t]-Sum[β[k,kk]*i[k,kk,t],{k,1,ages},{kk,1,stages}]*s[j,t] ;didt[j_,1,t_]=Sum[β[k,1]*i[k,1,t],{k,1,ages}]*Sum[r[nn,j,t],{nn,n+1-TRYm,n}]+Sum[β[k,1]*i[k,1,t],{k,1,ages}]*s[j,t] -κ[j,1]*i[j,1,t];didt[j_,2,t_]=Sum[β[k,2]*i[k,2,t],{k,1,ages}]*s[j,t] -κ[j,2]*i[j,2,t];drdt[1,j_,t_]=q*Sum[κ[j,jj]*i[j,jj,t],{jj,1,stages}]-n*δ[j]*r[1,j,t];For[nn=2,nn≤n-TRYm,nn++,
drdt[nn,j_,t_]=n*δ[j]*r[nn-1,j,t]-n*δ[j]*r[nn,j,t];
];
For[nn=n+1-TRYm,nn≤n,nn++,
drdt[nn,j_,t_]=n*δ[j]*r[nn-1,j,t]-Sum[β[k,1]*i[k,1,t],{k,1,ages}]*r[nn,j,t]-n*δ[j]*r[nn,j,t];
];
pars={β[j_,1]->βV,β[j_,2]->β,δ[1]->δY,δ[2]->δO,κ[j_,jj_]->κ}; 
frac[1,1]=(1-older)(1-f);
frac[2,1]=older(1-f);
frac[1,2]=(1-older)f;
frac[2,2]=older f;
nvars=Drop[Flatten[Table[{s[j,t],Table[i[j,jj,t],{jj,1,stages}],Table[r[nn,j,t],{nn,1,n}]},{j,1,ages}]],-1];
r[n,ages,t_]=1-Total[nvars];
neqns=Drop[Flatten[Table[{D[s[j,t],t]==dsdt[j,t],Table[D[i[j,jj,t],t]==didt[j,jj,t],{jj,1,stages}],Table[D[r[nn,j,t],t]==drdt[nn,j,t],{nn,1,n}]},{j,1,ages}]],-1];
nstart=Drop[Flatten[Table[{
Table[i[j,jj,0]==frac[j,jj] i0,{jj,1,stages}],
s[j,0]==(frac[j,1]+frac[j,2])*s0,Table[r[nn,j,0]==(frac[j,1]+frac[j,2]) (r0/n),{nn,1,n}]},{j,1,ages}]]/.start,-1];
NDSolve[Flatten[{nstart/.pars,neqns/.pars}],nvars,{t,0,maxtime}]
]

all=Flatten[Table[{j,jj},{j,1,ages},{jj,1,stages}],1]

labels={"<70,Res","<70,Mutant","70+,Res","70+,Mutant"};

colours={RGBColor[0.4,0.4,0.8],ColorData["Pastel",1],RGBColor[0.4,0.6,0.4], ColorData["Pastel",3/4]}

coltab=Join[{{1->{0,colours[[1]]}}},Table[{i->{{i-1},colours[[i]]}},{i,2,Length[all]}]];
Table[{labels[[i]],colours[[i]]},{i,1,Length[all]}]//MatrixForm

partemp={TRYp,TRYf,TRYolder,TRYβ,TRYβV,TRYδ,TRYδO,TRYκ,TRYconvert};

solution[partemp,start];

plot1=Show[
Plot[Evaluate@Table[100 Sum[i[all[[c,1]],all[[c,2]],t]/.solution[partemp,start],{c,1,jcum}],{jcum,1,Length[all]}],{t,0,maxtime},PlotRange->{Automatic,{0,maxy}},
PlotStyle->None,Filling->{{1->{0,}},{2->{{1},}},{3->{{2},}},{4->{{3},}}}],
Plot[100 ((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black,Dashed}],
Plot[100 (((TRYn/(TRYn-TRYm))TRYδ (TRYβV-TRYκ))/(TRYβV ((TRYn/(TRYn-TRYm))TRYδ+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black}],
ImageSize->250
]

Mutant fraction change between t=0 and t=50:

((Sum[i[j,1,t],{j,1,ages}])/(Sum[i[j,jj,t],{j,1,ages},{jj,1,stages}]))/.solution[partemp,start]/.t->50

Solve[Exp[s 50] ==( ((p50/(1-p50))/(p0/(1-p0)))/.p50->%/.p0->1-TRYf),s][[1]]

This is consistent with the predicted selection coefficient:

((TRYβV-TRYβ)/(TRYβ))TRYκ+(TRYm/TRYn)(((TRYκ TRYconvert(TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))))TRYβV//N

New equilibrium:

((((TRYδ+Δδ)(TRYβV-TRYκ))/(TRYβV ((TRYδ+Δδ)+TRYκ TRYconvert))))/.Δδ->(TRYm/(TRYn-TRYm))TRYδ//N

% rise in the endemic equilibrium:

100*(((%-TRYinf))/TRYinf)

The following illustrates the relatively minor effect caused by reducing seroconversion (to q = 80%), as long as δ/q is held constant:

partemp={TRYp,TRYf,TRYolder,TRYβ,TRYβV,(8/10)TRYδ,(8/10)TRYδO,TRYκ,(8/10)};

solution[partemp,start];

Showing the % of the population infectious:

plot2=Show[plot1,
Plot[Evaluate@Table[ Sum[100*i[all[[c,1]],all[[c,2]],t]/.solution[partemp,start],{c,1,jcum}],{jcum,1,Length[all]}],{t,0,maxtime},PlotRange->{Automatic,{0,maxy}},
PlotStyle->{,,,}],
Plot[100((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black,Dashed}],
Plot[100(((TRYn/(TRYn-TRYm))TRYδ (TRYβV-TRYκ))/(TRYβV ((TRYn/(TRYn-TRYm))TRYδ+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black}],
ImageSize->250
]

Seroconversion affects the dynamics (curves compared to the previous shaded plot) but not the initial speed of spread of the variant (described by the selection coefficient) or the equilibrium.

Clear[solution, i, s, r, ages, stages, n, neqns, nstart, nvars, didt, dsdt, drdt];

Figure S1B: strongly immune evasive (transient) but less transmissible (n=5)
-> leads to a decline in cases due to the spread of the variant, but eventually the previous variant takes over because of its transmission advantage
Can infect last three classes (m=3)

maxtime=500;
maxy=8;

Parameters:

TRYp=TRYp; (*No heterogeneity, so p not relevant now.*)
TRYf=99/100;(*Starting fraction of resident*)
TRYolder=13/100; (*~ Age 70+ from StatsCan 5022509/38929902 *)
	(*https://www150.statcan.gc.ca/t1/tbl1/en/cv.action?pid=1710000501*)
TRYκ=1/5; (*Five day infectious period*)
TRYδ=1/125; (*Waning rate based on Menegale et al. (2023)*) 
TRYδO=1/125; (*No clear age difference noted in Menegale et al. (2023)*) 
TRYconvert=10/10; (*Only this fraction of cases seroconvert or boost immunity*)

i->((δ (β-κ))/(β (δ+q κ)))

Solve[((δ (β-κ))/(β (δ+q κ)))==inf,β]

TRYinf=1/50; (*Assumed initial endemic frequency of infected individuals*)
TRYβ=((TRYδ TRYκ)/(TRYδ-TRYinf TRYδ-TRYconvert TRYinf TRYκ));

((TRYβ)/(TRYκ))//N

Now we allow stronger immune evasiveness, coupled with lower transmission.  Specifically, we assume five waning classes and the variant infects three earlier resistant classes (m=3):

TRYn=5;
TRYm=3;

Variant has a lower transmissibility (set to have the same selection coefficient as for the immune evasive variant):

Flatten[Solve[((βV-TRYβ)/(TRYβ))TRYκ+(TRYm/TRYn)(((TRYκ TRYconvert(TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))))βV==0.0833,βV]]

TRYβV=βV/.%;

TRYβ//N

Starting equilibrium would be (r0 includes all of the n resistance classes):
{s->(κ/β),i->((δ (β-κ))/(β (δ+q κ)))}

start={s0->((TRYκ)/(TRYβ)),i0->((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))),r0->1-((TRYκ)/(TRYβ))-((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert)))}

Waning distribution for a population of younger and older individuals:

Show[
Plot[TRYconvert*(1-CDF[GammaDistribution[TRYn,1/(TRYδ TRYn)],t]),{t,0,365},PlotRange->{All,{0,1}},PlotStyle->RGBColor[0.4,0.4,0.8]],
Plot[TRYconvert*(1-CDF[GammaDistribution[TRYn,1/(TRYδO TRYn)],t]),{t,0,365},PlotRange->{All,{0,1}},PlotStyle->RGBColor[0.4,0.6,0.4]]
]

Now the susceptible and resistant classes depend only on age (not whether they carried a variant before):
[NOTE: For the ease of colouring, the first stage class is the variant.]

{ages=2,stages=2,exposures=2,n=TRYn};

Clear[solution]
solution[{p_,f_,older_,β_,βV_,δY_,δO_,κ_,q_},start]:=solution[{p,f,older,β,βV,δY,δO,κ,q},start]=
Block[{ages=2,stages=2,exposures=2,n=TRYn},
dsdt[j_,zz_,t_]=(1-q)*Sum[κ[j,jj]*i[j,jj,zz,t],{jj,1,stages}]+n*δ[j]*r[n,j,zz,t]-Sum[β[k,kk]*i[k,kk,kkk,t],{k,1,ages},{kk,1,stages},{kkk,1,exposures}]*s[j,zz,t] ;
(*For the susceptible and recovered classes, we introduce zz to measure whether individuals have ever been exposed to the new variant (zz=1) or not (zz=2), this allows us to model transient immunity. We assume that those who do not seroconvert, with probability 1-c, remain in their zz class.*)didt[j_,1,1,t_]=Sum[β[k,1]*i[k,1,kkk,t],{k,1,ages},{kkk,1,exposures}]*Sum[r[nn,j,2,t],{nn,n+1-TRYm,n}]+Sum[β[k,1]*i[k,1,kkk,t],{k,1,ages},{kkk,1,exposures}]*s[j,1,t] -κ[j,1]*i[j,1,1,t];
(*This is the variant [stage 1], but it can only infect the recovering individuals who have never been exposed to a variant infection [stage 2]. Once infected by the variant these will move into the first type of recovering individual [zz=1] and are no longer available*)
didt[j_,1,2,t_]=Sum[β[k,1]*i[k,1,kkk,t],{k,1,ages},{kkk,1,exposures}]*s[j,2,t] -κ[j,1]*i[j,1,2,t];(*Although no infection with the variant lacks exposure to the variant, we allow this in case some fail to seroconvert, with probability 1-c.*)
didt[j_,2,zz_,t_]=Sum[β[k,2]*i[k,2,kkk,t],{k,1,ages},{kkk,1,exposures}]*s[j,zz,t]-κ[j,2]*i[j,2,zz,t];drdt[1,j_,1,t_]=q*(κ[j,1]*i[j,1,1,t]+κ[j,1]*i[j,1,2,t]+κ[j,2]*i[j,2,1,t])-n*δ[j]*r[1,j,1,t];
(*All infections with the variant move into this first resistance class only, as do those exposed in prior infections (i[j,2,1,t]).*)
drdt[1,j_,2,t_]=q*(κ[j,2]*i[j,2,2,t])-n*δ[j]*r[1,j,2,t];
For[nn=2,nn≤n-TRYm,nn++,
drdt[nn,j_,zz_,t_]=n*δ[j]*r[nn-1,j,zz,t]-n*δ[j]*r[nn,j,zz,t];
];
For[nn=n+1-TRYm,nn≤n,nn++,
drdt[nn,j_,1,t_]=n*δ[j]*r[nn-1,j,1,t]-n*δ[j]*r[nn,j,1,t];
drdt[nn,j_,2,t_]=n*δ[j]*r[nn-1,j,2,t]-Sum[β[k,1]*i[k,1,kkk,t],{k,1,ages},{kkk,1,exposures}]*r[nn,j,2,t]-n*δ[j]*r[nn,j,2,t];
(*Only those that have not been infected before by the new variant [stage 2] can become infected earlier by the new variant [stage 1]*)];
pars={β[j_,1]->βV,β[j_,2]->β,δ[1]->δY,δ[2]->δO,κ[j_,jj_]->κ}; 
frac[1,1]=(1-older)(1-f);
frac[2,1]=older(1-f);
frac[1,2]=(1-older)f;
frac[2,2]=older f;
nvars=Drop[Flatten[Table[{Table[s[j,zz,t],{zz,1,exposures}],
Table[i[j,jj,zz,t],{jj,1,stages},{zz,1,exposures}],
Table[r[nn,j,zz,t],{nn,1,n},{zz,1,exposures}]},{j,1,ages}]],-1];
r[n,ages,exposures,t_]=1-Total[nvars];
neqns=Drop[Flatten[Table[{Table[D[s[j,zz,t],t]==dsdt[j,zz,t],{zz,1,exposures}],Table[D[i[j,jj,zz,t],t]==didt[j,jj,zz,t],{jj,1,stages},{zz,1,exposures}],Table[D[r[nn,j,zz,t],t]==drdt[nn,j,zz,t],{nn,1,n},{zz,1,exposures}]},{j,1,ages}]],-1];
nstart=Drop[Flatten[Table[{
Table[i[j,jj,1,0]==0,{jj,1,stages}],
Table[i[j,jj,2,0]==frac[j,jj] i0,{jj,1,stages}],
s[j,1,0]==0,(*Nobody starts exposed to the variant*)
s[j,2,0]==(frac[j,1]+frac[j,2])*s0,
Table[r[nn,j,1,0]==0,{nn,1,n}],
Table[r[nn,j,2,0]==(frac[j,1]+frac[j,2]) (r0/n),{nn,1,n}]},{j,1,ages}]]/.start,-1];(*All susceptible and resistant individuals initially susceptible to the new variant [in stage 2] *)
NDSolve[Flatten[{nstart/.pars,neqns/.pars}],nvars,{t,0,maxtime}]
]

all=Flatten[Table[{j,jj},{j,1,ages},{jj,1,stages}],1]

labels={"<70,Res","<70,Mutant","70+,Res","70+,Mutant"};

colours={RGBColor[0.4,0.4,0.8],ColorData["Pastel",1],RGBColor[0.4,0.6,0.4], ColorData["Pastel",3/4]}

coltab=Join[{{1->{0,colours[[1]]}}},Table[{i->{{i-1},colours[[i]]}},{i,2,Length[all]}]];
Table[{labels[[i]],colours[[i]]},{i,1,Length[all]}]//MatrixForm

partemp={TRYp,TRYf,TRYolder,TRYβ,TRYβV,TRYδ,TRYδO,TRYκ,TRYconvert};

solution[partemp,start];

Showing the % of the population infectious:

plot1=Show[
Plot[Evaluate@Table[ Sum[100i[all[[c,1]],all[[c,2]],1,t]+100i[all[[c,1]],all[[c,2]],2,t]/.solution[partemp,start],{c,1,jcum}],{jcum,1,Length[all]}],{t,0,maxtime},PlotRange->{Automatic,{0,maxy}},
PlotStyle->None,Filling->{{1->{0,}},{2->{{1},}},{3->{{2},}},{4->{{3},}}}],
Plot[100((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black,Dashed}],
Plot[100(((TRYn/(TRYn-0))TRYδ (TRYβV-TRYκ))/(TRYβV ((TRYn/(TRYn-0))TRYδ+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black}],
ImageSize->250
]

Mutant fraction change between t=0 and t=50:

((Sum[i[j,1,zz,t],{j,1,ages},{zz,1,exposures}])/(Sum[i[j,jj,zz,t],{j,1,ages},{jj,1,stages},{zz,1,exposures}]))/.solution[partemp,start]/.t->50

Solve[Exp[s 50] ==( ((p50/(1-p50))/(p0/(1-p0)))/.p50->%/.p0->1-TRYf),s][[1]]

This is consistent with the predicted selection coefficient:

((TRYβV-TRYβ)/(TRYβ))TRYκ+(TRYm/TRYn)(((TRYκ TRYconvert(TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))))TRYβV//N

New equilibrium (if resident fixed):

((((TRYδ+Δδ)(TRYβV-TRYκ))/(TRYβV ((TRYδ+Δδ)+TRYκ TRYconvert))))/.Δδ->(0/(TRYn-0))TRYδ//N

% rise in the endemic equilibrium:

100*(((%-TRYinf))/TRYinf)

The following illustrates the relatively minor effect caused by reducing seroconversion (to q = 80%), as long as δ/q is held constant:

partemp={TRYp,TRYf,TRYolder,TRYβ,TRYβV,(8/10)TRYδ,(8/10)TRYδO,TRYκ,(8/10)};

solution[partemp,start];

Showing the % of the population infectious:

plot2=Show[plot1,
Plot[Evaluate@Table[ Sum[100*i[all[[c,1]],all[[c,2]],t]/.solution[partemp,start],{c,1,jcum}],{jcum,1,Length[all]}],{t,0,maxtime},PlotRange->{Automatic,{0,maxy}},
PlotStyle->{,,,}],
Plot[100((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black,Dashed}],
Plot[100(((TRYn/(TRYn-0))TRYδ (TRYβV-TRYκ))/(TRYβV ((TRYn/(TRYn-0))TRYδ+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black}],
ImageSize->250
]

Seroconversion affects the dynamics (curves compared to the previous shaded plot) but not the initial speed of spread of the variant (described by the selection coefficient) or the equilibrium.

Clear[solution, i, s, r, ages, stages, n, neqns, nstart, nvars, didt, dsdt, drdt];

Figure S2A: more transmissible (n=1)
-> Less oscillatory

maxtime=500;
maxy=8;

Parameters:

TRYp=TRYp; (*No heterogeneity, so p not relevant now.*)
TRYf=99/100;(*Starting fraction of resident*)
TRYolder=13/100; (*~ Age 70+ from StatsCan 5022509/38929902 *)
	(*https://www150.statcan.gc.ca/t1/tbl1/en/cv.action?pid=1710000501*)
TRYκ=1/5; (*Five day infectious period*)
TRYδ=1/125; (*Waning rate based on Menegale et al. (2023)*) 
TRYδO=1/125; (*No clear age difference noted in Menegale et al. (2023)*) 
TRYconvert=10/10; (*Only this fraction of cases seroconvert or boost immunity*)

i->((δ (β-κ))/(β (δ+q κ)))

Solve[((δ (β-κ))/(β (δ+q κ)))==inf,β]

TRYinf=1/50; (*Assumed initial endemic frequency of infected individuals*)
TRYβ=((TRYδ TRYκ)/(TRYδ-TRYinf TRYδ-TRYconvert TRYinf TRYκ));

((TRYβ)/(TRYκ))//N

Variant has a higher transmissibility (set to have the same selection coefficient as for the immune evasive variant):

Flatten[Solve[((βV-TRYβ)/(TRYβ))TRYκ==0.0833,βV]]

TRYβV=βV/.%;

% increase in β:

100*((TRYβV-TRYβ)/(TRYβ))

and now has no immune evasive properties.  Specifically, we assume five waning classes and the variant infects only susceptibles (m=0):

TRYn=1;
TRYm=0;

Starting equilibrium would be (r0 includes all of the n resistance classes):
{s->(κ/β),i->((δ (β-κ))/(β (δ+q κ)))}

start={s0->((TRYκ)/(TRYβ)),i0->((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))),r0->1-((TRYκ)/(TRYβ))-((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert)))}

Waning distribution for a population of younger and older individuals:

Show[
Plot[TRYconvert*(1-CDF[GammaDistribution[TRYn,1/(TRYδ TRYn)],t]),{t,0,365},PlotRange->{All,{0,1}},PlotStyle->RGBColor[0.4,0.4,0.8]],
Plot[TRYconvert*(1-CDF[GammaDistribution[TRYn,1/(TRYδO TRYn)],t]),{t,0,365},PlotRange->{All,{0,1}},PlotStyle->RGBColor[0.4,0.6,0.4]]
]

Now the susceptible and resistant classes depend only on age (not whether they carried a variant before):
[NOTE: For the ease of colouring, the first stage class is the variant.]

{ages=2,stages=2,n=TRYn};

Clear[solution]
solution[{p_,f_,older_,β_,βV_,δY_,δO_,κ_,q_},start]:=solution[{p,f,older,β,βV,δY,δO,κ,q},start]=
Block[{ages=2,stages=2,n=TRYn},
dsdt[j_,t_]=(1-q)*Sum[κ[j,jj]*i[j,jj,t],{jj,1,stages}]+n*δ[j]*r[n,j,t]-Sum[β[k,kk]*i[k,kk,t],{k,1,ages},{kk,1,stages}]*s[j,t] ;didt[j_,1,t_]=Sum[β[k,1]*i[k,1,t],{k,1,ages}]*Sum[r[nn,j,t],{nn,n+1-TRYm,n}]+Sum[β[k,1]*i[k,1,t],{k,1,ages}]*s[j,t] -κ[j,1]*i[j,1,t];didt[j_,2,t_]=Sum[β[k,2]*i[k,2,t],{k,1,ages}]*s[j,t] -κ[j,2]*i[j,2,t];drdt[1,j_,t_]=q*Sum[κ[j,jj]*i[j,jj,t],{jj,1,stages}]-n*δ[j]*r[1,j,t];For[nn=2,nn≤n-TRYm,nn++,
drdt[nn,j_,t_]=n*δ[j]*r[nn-1,j,t]-n*δ[j]*r[nn,j,t];
];
For[nn=n+1-TRYm,nn≤n,nn++,
drdt[nn,j_,t_]=n*δ[j]*r[nn-1,j,t]-Sum[β[k,1]*i[k,1,t],{k,1,ages}]*r[nn,j,t]-n*δ[j]*r[nn,j,t];
];
pars={β[j_,1]->βV,β[j_,2]->β,δ[1]->δY,δ[2]->δO,κ[j_,jj_]->κ}; 
frac[1,1]=(1-older)(1-f);
frac[2,1]=older(1-f);
frac[1,2]=(1-older)f;
frac[2,2]=older f;
nvars=Drop[Flatten[Table[{s[j,t],Table[i[j,jj,t],{jj,1,stages}],Table[r[nn,j,t],{nn,1,n}]},{j,1,ages}]],-1];
r[n,ages,t_]=1-Total[nvars];
neqns=Drop[Flatten[Table[{D[s[j,t],t]==dsdt[j,t],Table[D[i[j,jj,t],t]==didt[j,jj,t],{jj,1,stages}],Table[D[r[nn,j,t],t]==drdt[nn,j,t],{nn,1,n}]},{j,1,ages}]],-1];
nstart=Drop[Flatten[Table[{
Table[i[j,jj,0]==frac[j,jj] i0,{jj,1,stages}],
s[j,0]==(frac[j,1]+frac[j,2])*s0,Table[r[nn,j,0]==(frac[j,1]+frac[j,2]) (r0/n),{nn,1,n}]},{j,1,ages}]]/.start,-1];
NDSolve[Flatten[{nstart/.pars,neqns/.pars}],nvars,{t,0,maxtime}]
]

all=Flatten[Table[{j,jj},{j,1,ages},{jj,1,stages}],1]

labels={"<70,Res","<70,Mutant","70+,Res","70+,Mutant"};

colours={RGBColor[0.4,0.4,0.8],ColorData["Pastel",1],RGBColor[0.4,0.6,0.4], ColorData["Pastel",3/4]}

coltab=Join[{{1->{0,colours[[1]]}}},Table[{i->{{i-1},colours[[i]]}},{i,2,Length[all]}]];
Table[{labels[[i]],colours[[i]]},{i,1,Length[all]}]//MatrixForm

partemp={TRYp,TRYf,TRYolder,TRYβ,TRYβV,TRYδ,TRYδO,TRYκ,TRYconvert};

solution[partemp,start];

Showing the % of the population infectious:

plot1=Show[
Plot[Evaluate@Table[ Sum[100 i[all[[c,1]],all[[c,2]],t]/.solution[partemp,start],{c,1,jcum}],{jcum,1,Length[all]}],{t,0,maxtime},PlotRange->{Automatic,{0,maxy}},
PlotStyle->None,Filling->{{1->{0,}},{2->{{1},}},{3->{{2},}},{4->{{3},}}}],
Plot[100((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black,Dashed}],
Plot[100(((TRYn/(TRYn-TRYm))TRYδ (TRYβV-TRYκ))/(TRYβV ((TRYn/(TRYn-TRYm))TRYδ+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black}],
ImageSize->250
]

Mutant fraction change between t=0 and t=50:

((Sum[i[j,1,t],{j,1,ages}])/(Sum[i[j,jj,t],{j,1,ages},{jj,1,stages}]))/.solution[partemp,start]/.t->50

Solve[Exp[s 50] ==( ((p50/(1-p50))/(p0/(1-p0)))/.p50->%/.p0->1-TRYf),s][[1]]

This is consistent with the predicted selection coefficient:

((TRYβV-TRYβ)/(TRYβ))TRYκ+(TRYm/TRYn)(((TRYκ TRYconvert(TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))))TRYβV//N

New equilibrium:

((((TRYδ+Δδ)(TRYβV-TRYκ))/(TRYβV ((TRYδ+Δδ)+TRYκ TRYconvert))))/.Δδ->(TRYm/(TRYn-TRYm))TRYδ//N

% rise in the endemic equilibrium:

100*(((%-TRYinf))/TRYinf)

This considers the number of susceptibles, with the black line at 1/(Subscript[(OverTilde[R]), 0]) (the fraction of susceptibles that causes the disease to stop growing before the variant) and the dashed black line at 1/(Subscript[(OverTilde[R]), 0]^( * )) (the fraction of susceptibles that causes the disease to stop growing once the variant predominates):

Show[
Plot[Evaluate@Table[ Sum[100 s[c,t]/.solution[partemp,start],{c,1,jcum}],{jcum,1,2}],{t,0,maxtime},PlotRange->{Automatic,{0,100}},
PlotStyle->None,Filling->{{1->{0,}},{2->{{1},}}}],
Plot[100((TRYκ)/(TRYβ)),{t,0,maxtime},PlotStyle->{Black}],
Plot[100((TRYκ)/(TRYβV)),{t,0,maxtime},PlotStyle->{Black,Dashed}],
ImageSize->250
]

The first time that the number of susceptibles falls below ((TRYκ)/(TRYβV))is when we expect the number of the new variant to start declining (roughly, given that waning changes the dynamics):

(1+(((TRYβV-TRYβ))/(TRYβ)))

droptime=First[Select[Table[Flatten[{t, Sum[100 s[c,t]/.solution[partemp,start],{c,1,2}]}],{t,0,maxtime}],#[[2]]<100((TRYκ)/(TRYβV))&]]

Show[
plot1,
ListPlot[{{droptime[[1]],0},{droptime[[1]],maxy}},Joined->True,PlotStyle->{Gray}],
ImageSize->250
]

The following illustrates the relatively minor effect caused by reducing seroconversion (to q = 80%), as long as δ/q is held constant:

partemp={TRYp,TRYf,TRYolder,TRYβ,TRYβV,(8/10)TRYδ,(8/10)TRYδO,TRYκ,(8/10)};

solution[partemp,start];

Showing the % of the population infectious:

plot2=Show[plot1,
Plot[Evaluate@Table[ Sum[100*i[all[[c,1]],all[[c,2]],t]/.solution[partemp,start],{c,1,jcum}],{jcum,1,Length[all]}],{t,0,maxtime},PlotRange->{Automatic,{0,maxy}},
PlotStyle->{,,,}],
Plot[100((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black,Dashed}],
Plot[100(((TRYn/(TRYn-0))TRYδ (TRYβV-TRYκ))/(TRYβV ((TRYn/(TRYn-0))TRYδ+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black}],
ImageSize->250
]

Seroconversion affects the dynamics (curves compared to the previous shaded plot) but not the initial speed of spread of the variant (described by the selection coefficient) or the equilibrium.

Clear[solution, i, s, r, ages, stages, n, neqns, nstart, nvars, didt, dsdt, drdt];

Figure S2C: more transmissible (n=10)
-> More oscillatory

maxtime=500;
maxy=8;

Parameters:

TRYp=TRYp; (*No heterogeneity, so p not relevant now.*)
TRYf=99/100;(*Starting fraction of resident*)
TRYolder=13/100; (*~ Age 70+ from StatsCan 5022509/38929902 *)
	(*https://www150.statcan.gc.ca/t1/tbl1/en/cv.action?pid=1710000501*)
TRYκ=1/5; (*Five day infectious period*)
TRYδ=1/125; (*Waning rate based on Menegale et al. (2023)*) 
TRYδO=1/125; (*No clear age difference noted in Menegale et al. (2023)*) 
TRYconvert=10/10; (*Only this fraction of cases seroconvert or boost immunity*)

i->((δ (β-κ))/(β (δ+q κ)))

Solve[((δ (β-κ))/(β (δ+q κ)))==inf,β]

TRYinf=1/50; (*Assumed initial endemic frequency of infected individuals*)
TRYβ=((TRYδ TRYκ)/(TRYδ-TRYinf TRYδ-TRYconvert TRYinf TRYκ));

((TRYβ)/(TRYκ))//N

Variant has a higher transmissibility (set to have the same selection coefficient as for the immune evasive variant):

Flatten[Solve[((βV-TRYβ)/(TRYβ))TRYκ==0.0833,βV]]

TRYβV=βV/.%;

% increase in β:

100*((TRYβV-TRYβ)/(TRYβ))

and now has no immune evasive properties.  Specifically, we assume five waning classes and the variant infects only susceptibles (m=0):

TRYn=10;
TRYm=0;

Starting equilibrium would be (r0 includes all of the n resistance classes):
{s->(κ/β),i->((δ (β-κ))/(β (δ+q κ)))}

start={s0->((TRYκ)/(TRYβ)),i0->((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))),r0->1-((TRYκ)/(TRYβ))-((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert)))}

Waning distribution for a population of younger and older individuals:

Show[
Plot[TRYconvert*(1-CDF[GammaDistribution[TRYn,1/(TRYδ TRYn)],t]),{t,0,365},PlotRange->{All,{0,1}},PlotStyle->RGBColor[0.4,0.4,0.8]],
Plot[TRYconvert*(1-CDF[GammaDistribution[TRYn,1/(TRYδO TRYn)],t]),{t,0,365},PlotRange->{All,{0,1}},PlotStyle->RGBColor[0.4,0.6,0.4]]
]

Now the susceptible and resistant classes depend only on age (not whether they carried a variant before):
[NOTE: For the ease of colouring, the first stage class is the variant.]

{ages=2,stages=2,n=TRYn};

Clear[solution]
solution[{p_,f_,older_,β_,βV_,δY_,δO_,κ_,q_},start]:=solution[{p,f,older,β,βV,δY,δO,κ,q},start]=
Block[{ages=2,stages=2,n=TRYn},
dsdt[j_,t_]=(1-q)*Sum[κ[j,jj]*i[j,jj,t],{jj,1,stages}]+n*δ[j]*r[n,j,t]-Sum[β[k,kk]*i[k,kk,t],{k,1,ages},{kk,1,stages}]*s[j,t] ;didt[j_,1,t_]=Sum[β[k,1]*i[k,1,t],{k,1,ages}]*Sum[r[nn,j,t],{nn,n+1-TRYm,n}]+Sum[β[k,1]*i[k,1,t],{k,1,ages}]*s[j,t] -κ[j,1]*i[j,1,t];didt[j_,2,t_]=Sum[β[k,2]*i[k,2,t],{k,1,ages}]*s[j,t] -κ[j,2]*i[j,2,t];drdt[1,j_,t_]=q*Sum[κ[j,jj]*i[j,jj,t],{jj,1,stages}]-n*δ[j]*r[1,j,t];For[nn=2,nn≤n-TRYm,nn++,
drdt[nn,j_,t_]=n*δ[j]*r[nn-1,j,t]-n*δ[j]*r[nn,j,t];
];
For[nn=n+1-TRYm,nn≤n,nn++,
drdt[nn,j_,t_]=n*δ[j]*r[nn-1,j,t]-Sum[β[k,1]*i[k,1,t],{k,1,ages}]*r[nn,j,t]-n*δ[j]*r[nn,j,t];
];
pars={β[j_,1]->βV,β[j_,2]->β,δ[1]->δY,δ[2]->δO,κ[j_,jj_]->κ}; 
frac[1,1]=(1-older)(1-f);
frac[2,1]=older(1-f);
frac[1,2]=(1-older)f;
frac[2,2]=older f;
nvars=Drop[Flatten[Table[{s[j,t],Table[i[j,jj,t],{jj,1,stages}],Table[r[nn,j,t],{nn,1,n}]},{j,1,ages}]],-1];
r[n,ages,t_]=1-Total[nvars];
neqns=Drop[Flatten[Table[{D[s[j,t],t]==dsdt[j,t],Table[D[i[j,jj,t],t]==didt[j,jj,t],{jj,1,stages}],Table[D[r[nn,j,t],t]==drdt[nn,j,t],{nn,1,n}]},{j,1,ages}]],-1];
nstart=Drop[Flatten[Table[{
Table[i[j,jj,0]==frac[j,jj] i0,{jj,1,stages}],
s[j,0]==(frac[j,1]+frac[j,2])*s0,Table[r[nn,j,0]==(frac[j,1]+frac[j,2]) (r0/n),{nn,1,n}]},{j,1,ages}]]/.start,-1];
NDSolve[Flatten[{nstart/.pars,neqns/.pars}],nvars,{t,0,maxtime}]
]

all=Flatten[Table[{j,jj},{j,1,ages},{jj,1,stages}],1]

labels={"<70,Res","<70,Mutant","70+,Res","70+,Mutant"};

colours={RGBColor[0.4,0.4,0.8],ColorData["Pastel",1],RGBColor[0.4,0.6,0.4], ColorData["Pastel",3/4]}

coltab=Join[{{1->{0,colours[[1]]}}},Table[{i->{{i-1},colours[[i]]}},{i,2,Length[all]}]];
Table[{labels[[i]],colours[[i]]},{i,1,Length[all]}]//MatrixForm

partemp={TRYp,TRYf,TRYolder,TRYβ,TRYβV,TRYδ,TRYδO,TRYκ,TRYconvert};

solution[partemp,start];

Showing the % of the population infectious:

plot1=Show[
Plot[Evaluate@Table[ Sum[100 i[all[[c,1]],all[[c,2]],t]/.solution[partemp,start],{c,1,jcum}],{jcum,1,Length[all]}],{t,0,maxtime},PlotRange->{Automatic,{0,maxy}},
PlotStyle->None,Filling->{{1->{0,}},{2->{{1},}},{3->{{2},}},{4->{{3},}}}],
Plot[100((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black,Dashed}],
Plot[100(((TRYn/(TRYn-TRYm))TRYδ (TRYβV-TRYκ))/(TRYβV ((TRYn/(TRYn-TRYm))TRYδ+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black}],
ImageSize->250
]

Mutant fraction change between t=0 and t=50:

((Sum[i[j,1,t],{j,1,ages}])/(Sum[i[j,jj,t],{j,1,ages},{jj,1,stages}]))/.solution[partemp,start]/.t->50

Solve[Exp[s 50] ==( ((p50/(1-p50))/(p0/(1-p0)))/.p50->%/.p0->1-TRYf),s][[1]]

This is consistent with the predicted selection coefficient:

((TRYβV-TRYβ)/(TRYβ))TRYκ+(TRYm/TRYn)(((TRYκ TRYconvert(TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))))TRYβV//N

New equilibrium:

((((TRYδ+Δδ)(TRYβV-TRYκ))/(TRYβV ((TRYδ+Δδ)+TRYκ TRYconvert))))/.Δδ->(TRYm/(TRYn-TRYm))TRYδ//N

% rise in the endemic equilibrium:

100*(((%-TRYinf))/TRYinf)

This considers the number of susceptibles, with the black line at 1/(Subscript[(OverTilde[R]), 0]) (the fraction of susceptibles that causes the disease to stop growing before the variant) and the dashed black line at 1/(Subscript[(OverTilde[R]), 0]^( * )) (the fraction of susceptibles that causes the disease to stop growing once the variant predominates):

Show[
Plot[Evaluate@Table[ Sum[100 s[c,t]/.solution[partemp,start],{c,1,jcum}],{jcum,1,2}],{t,0,maxtime},PlotRange->{Automatic,{0,100}},
PlotStyle->None,Filling->{{1->{0,}},{2->{{1},}}}],
Plot[100((TRYκ)/(TRYβ)),{t,0,maxtime},PlotStyle->{Black}],
Plot[100((TRYκ)/(TRYβV)),{t,0,maxtime},PlotStyle->{Black,Dashed}],
ImageSize->250
]

The first time that the number of susceptibles falls below ((TRYκ)/(TRYβV))is when we expect the number of the new variant to start declining (roughly, given that waning changes the dynamics):

(1+(((TRYβV-TRYβ))/(TRYβ)))

droptime=First[Select[Table[Flatten[{t, Sum[100 s[c,t]/.solution[partemp,start],{c,1,2}]}],{t,0,maxtime}],#[[2]]<100((TRYκ)/(TRYβV))&]]

Show[
plot1,
ListPlot[{{droptime[[1]],0},{droptime[[1]],maxy}},Joined->True,PlotStyle->{Gray}],
ImageSize->250
]

The following illustrates the relatively minor effect caused by reducing seroconversion (to q = 80%), as long as δ/q is held constant:

partemp={TRYp,TRYf,TRYolder,TRYβ,TRYβV,(8/10)TRYδ,(8/10)TRYδO,TRYκ,(8/10)};

solution[partemp,start];

Showing the % of the population infectious:

plot2=Show[plot1,
Plot[Evaluate@Table[ Sum[100*i[all[[c,1]],all[[c,2]],t]/.solution[partemp,start],{c,1,jcum}],{jcum,1,Length[all]}],{t,0,maxtime},PlotRange->{Automatic,{0,maxy}},
PlotStyle->{,,,}],
Plot[100((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black,Dashed}],
Plot[100(((TRYn/(TRYn-0))TRYδ (TRYβV-TRYκ))/(TRYβV ((TRYn/(TRYn-0))TRYδ+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black}],
ImageSize->250
]

Seroconversion affects the dynamics (curves compared to the previous shaded plot) but not the initial speed of spread of the variant (described by the selection coefficient) or the equilibrium.

Clear[solution, i, s, r, ages, stages, n, neqns, nstart, nvars, didt, dsdt, drdt];

Alt Figure 2B: persistently immune evasive (n=10) (increases mean waning rate by 67%)
Can infect last four classes (m=4)
-> More oscillatory

maxtime=500;
maxy=8;

Parameters:

TRYp=TRYp; (*No heterogeneity, so p not relevant now.*)
TRYf=99/100;(*Starting fraction of resident*)
TRYolder=13/100; (*~ Age 70+ from StatsCan 5022509/38929902 *)
	(*https://www150.statcan.gc.ca/t1/tbl1/en/cv.action?pid=1710000501*)
TRYκ=1/5; (*Five day infectious period*)
TRYδ=1/125; (*Waning rate based on Menegale et al. (2023)*) 
TRYδO=1/125; (*No clear age difference noted in Menegale et al. (2023)*) 
TRYconvert=10/10; (*Only this fraction of cases seroconvert or boost immunity*)

i->((δ (β-κ))/(β (δ+q κ)))

Solve[((δ (β-κ))/(β (δ+q κ)))==inf,β]

TRYinf=1/50; (*Assumed initial endemic frequency of infected individuals*)
TRYβ=((TRYδ TRYκ)/(TRYδ-TRYinf TRYδ-TRYconvert TRYinf TRYκ));

((TRYβ)/(TRYκ))//N

Variant has the same transmissibility:

TRYβV=TRYβ;

but infects earlier during waning.  Specifically, we assume five waning classes and the variant infects two earlier.

TRYn=10;
TRYm=4;

Thus, the mean waning rate (once fixed) increases by 67%:

((TRYδ+Δδ)/(TRYδ))/.Δδ->(TRYm/(TRYn-TRYm))TRYδ//N

while the mean waiting time decreases by 40%:

1/%

Starting equilibrium would be (r0 includes all of the n resistance classes):
{s->(κ/β),i->((δ (β-κ))/(β (δ+q κ)))}

start={s0->((TRYκ)/(TRYβ)),i0->((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))),r0->1-((TRYκ)/(TRYβ))-((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert)))}

Waning distribution for a population of younger and older individuals:

Show[
Plot[TRYconvert*(1-CDF[GammaDistribution[TRYn,1/(TRYδ TRYn)],t]),{t,0,365},PlotRange->{All,{0,1}},PlotStyle->RGBColor[0.4,0.4,0.8]],
Plot[TRYconvert*(1-CDF[GammaDistribution[TRYn,1/(TRYδO TRYn)],t]),{t,0,365},PlotRange->{All,{0,1}},PlotStyle->RGBColor[0.4,0.6,0.4]]
]

Now the susceptible and resistant classes depend only on age (not whether they carried a variant before):
[NOTE: For the ease of colouring, the first stage class is the variant.]

{ages=2,stages=2,n=TRYn};

Clear[solution]
solution[{p_,f_,older_,β_,βV_,δY_,δO_,κ_,q_},start]:=solution[{p,f,older,β,βV,δY,δO,κ,q},start]=
Block[{ages=2,stages=2,n=TRYn},
dsdt[j_,t_]=(1-q)*Sum[κ[j,jj]*i[j,jj,t],{jj,1,stages}]+n*δ[j]*r[n,j,t]-Sum[β[k,kk]*i[k,kk,t],{k,1,ages},{kk,1,stages}]*s[j,t] ;didt[j_,1,t_]=Sum[β[k,1]*i[k,1,t],{k,1,ages}]*Sum[r[nn,j,t],{nn,n+1-TRYm,n}]+Sum[β[k,1]*i[k,1,t],{k,1,ages}]*s[j,t] -κ[j,1]*i[j,1,t];didt[j_,2,t_]=Sum[β[k,2]*i[k,2,t],{k,1,ages}]*s[j,t] -κ[j,2]*i[j,2,t];drdt[1,j_,t_]=q*Sum[κ[j,jj]*i[j,jj,t],{jj,1,stages}]-n*δ[j]*r[1,j,t];For[nn=2,nn≤n-TRYm,nn++,
drdt[nn,j_,t_]=n*δ[j]*r[nn-1,j,t]-n*δ[j]*r[nn,j,t];
];
For[nn=n+1-TRYm,nn≤n,nn++,
drdt[nn,j_,t_]=n*δ[j]*r[nn-1,j,t]-Sum[β[k,1]*i[k,1,t],{k,1,ages}]*r[nn,j,t]-n*δ[j]*r[nn,j,t];
];
pars={β[j_,1]->βV,β[j_,2]->β,δ[1]->δY,δ[2]->δO,κ[j_,jj_]->κ}; 
frac[1,1]=(1-older)(1-f);
frac[2,1]=older(1-f);
frac[1,2]=(1-older)f;
frac[2,2]=older f;
nvars=Drop[Flatten[Table[{s[j,t],Table[i[j,jj,t],{jj,1,stages}],Table[r[nn,j,t],{nn,1,n}]},{j,1,ages}]],-1];
r[n,ages,t_]=1-Total[nvars];
neqns=Drop[Flatten[Table[{D[s[j,t],t]==dsdt[j,t],Table[D[i[j,jj,t],t]==didt[j,jj,t],{jj,1,stages}],Table[D[r[nn,j,t],t]==drdt[nn,j,t],{nn,1,n}]},{j,1,ages}]],-1];
nstart=Drop[Flatten[Table[{
Table[i[j,jj,0]==frac[j,jj] i0,{jj,1,stages}],
s[j,0]==(frac[j,1]+frac[j,2])*s0,Table[r[nn,j,0]==(frac[j,1]+frac[j,2]) (r0/n),{nn,1,n}]},{j,1,ages}]]/.start,-1];
NDSolve[Flatten[{nstart/.pars,neqns/.pars}],nvars,{t,0,maxtime}]
]

all=Flatten[Table[{j,jj},{j,1,ages},{jj,1,stages}],1]

labels={"<70,Res","<70,Mutant","70+,Res","70+,Mutant"};

colours={RGBColor[0.4,0.4,0.8],ColorData["Pastel",1],RGBColor[0.4,0.6,0.4], ColorData["Pastel",3/4]}

coltab=Join[{{1->{0,colours[[1]]}}},Table[{i->{{i-1},colours[[i]]}},{i,2,Length[all]}]];
Table[{labels[[i]],colours[[i]]},{i,1,Length[all]}]//MatrixForm

partemp={TRYp,TRYf,TRYolder,TRYβ,TRYβV,TRYδ,TRYδO,TRYκ,TRYconvert};

solution[partemp,start];

Showing the % of the population infectious:

plot1 =Show[
Plot[Evaluate@Table[ Sum[100*i[all[[c,1]],all[[c,2]],t]/.solution[partemp,start],{c,1,jcum}],{jcum,1,Length[all]}],{t,0,maxtime},PlotRange->{Automatic,{0,maxy}},
PlotStyle->None,Filling->{{1->{0,}},{2->{{1},}},{3->{{2},}},{4->{{3},}}}],
Plot[100((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black,Dashed}],
Plot[100(((TRYn/(TRYn-TRYm))TRYδ (TRYβV-TRYκ))/(TRYβV ((TRYn/(TRYn-TRYm))TRYδ+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black}],
ImageSize->250
]

Mutant fraction change between t=0 and t=50:

((Sum[i[j,1,t],{j,1,ages}])/(Sum[i[j,jj,t],{j,1,ages},{jj,1,stages}]))/.solution[partemp,start]/.t->50

Solve[Exp[s 50] ==( ((p50/(1-p50))/(p0/(1-p0)))/.p50->%/.p0->1-TRYf),s][[1]]

This is consistent with the predicted selection coefficient:

((TRYβV-TRYβ)/(TRYβ))TRYκ+(TRYm/TRYn)(((TRYκ TRYconvert(TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))))TRYβV//N

New equilibrium:

((((TRYδ+Δδ)(TRYβV-TRYκ))/(TRYβV ((TRYδ+Δδ)+TRYκ TRYconvert))))/.Δδ->(TRYm/(TRYn-TRYm))TRYδ//N

% rise in the endemic equilibrium:

100*(((%-TRYinf))/TRYinf)

The following illustrates the relatively minor effect caused by reducing seroconversion (to q = 80%), as long as δ/q is held constant:

partemp={TRYp,TRYf,TRYolder,TRYβ,TRYβV,(8/10)TRYδ,(8/10)TRYδO,TRYκ,(8/10)};

solution[partemp,start];

Showing the % of the population infectious:

plot2=Show[plot1,
Plot[Evaluate@Table[ Sum[100*i[all[[c,1]],all[[c,2]],t]/.solution[partemp,start],{c,1,jcum}],{jcum,1,Length[all]}],{t,0,maxtime},PlotRange->{Automatic,{0,maxy}},
PlotStyle->{,,,}],
Plot[100((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black,Dashed}],
Plot[100(((TRYn/(TRYn-TRYm))TRYδ (TRYβV-TRYκ))/(TRYβV ((TRYn/(TRYn-TRYm))TRYδ+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black}],
ImageSize->250
]

Seroconversion affects the dynamics (curves compared to the previous shaded plot) but not the initial speed of spread of the variant (described by the selection coefficient) or the equilibrium.

Clear[solution, i, s, r, ages, stages, n, neqns, nstart, nvars, didt, dsdt, drdt];

Vaccines

Figure 6: Change in vaccinations

We consider two vaccination rates, 0.012% corresponding to the low summer rate (May-July 2023) and  0.174% corresponding to the high fall rate (Sept-Dec 2022) in Canada (see vaccination-administration.xlsx from https://health-infobase.canada.ca/covid-19/vaccine-administration/).

Below we will explore the impact of each level of vaccination, as well as a policy limiting vaccination to older individuals.

TRYvLOW=1/8333; (* 0.00012 vaccinations per person per day*)
TRYvHIGH=1/575; (* 0.00174 vaccinations per person per day*)

maxtime=500;
maxy=3;

Parameters:

TRYp=TRYp; (*No heterogeneity, so p not relevant now.*)
TRYf=1;(*Starting fraction of resident*)
TRYolder=13/100; (*~ Age 70+ from StatsCan 5022509/38929902 *)
	(*https://www150.statcan.gc.ca/t1/tbl1/en/cv.action?pid=1710000501*)
TRYκ=1/5; (*Five day infectious period*)
TRYδ=1/125; (*Waning rate based on Menegale et al. (2023)*) 
TRYδO=1/125; (*No clear age difference noted in Menegale et al. (2023)*) 
TRYconvert=10/10; (*Only this fraction of cases seroconvert or boost immunity*)

i->((δ (β-κ))/(β (δ+q κ)))

Solve[((δ (β-κ))/(β (δ+q κ)))==inf,β]

Incidence is expected to decline by:

(v/(δ+κ))/.v->{ TRYvLOW,TRYvHIGH}/.δ->TRYδ/.κ->TRYκ//N

TRYinf=1/50; (*Assumed initial endemic frequency of infected individuals*)
TRYβ=((TRYδ TRYκ)/(TRYδ-TRYinf TRYδ-TRYconvert TRYinf TRYκ));

((TRYβ)/(TRYκ))//N

We assume five waning classes and no variant (m=0):

TRYn=5;
TRYm=0;

TRYβV=TRYβ;

Starting equilibrium would be (r0 includes all of the n resistance classes):
{s->(κ/β),i->((δ (β-κ))/(β (δ+q κ)))}

start={s0->((TRYκ)/(TRYβ)),i0->((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))),r0->1-((TRYκ)/(TRYβ))-((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert)))}

[NOTE: For the ease of colouring, the first stage class is the variant and the second the resident (here there are no variants).]

{ages=2,stages=2,n=TRYn};

Clear[solution]
solution[{p_,f_,older_,β_,βV_,δY_,δO_,vY_,vO_,κ_,q_},start]:=
solution[{p,f,older,β,βV,δY,δO,vY,vO,κ,q},start]=
Block[{ages=2,stages=2,n=TRYn},
dsdt[j_,t_]=(1-q)*Sum[κ[j,jj]*i[j,jj,t],{jj,1,stages}]+n*δ[j]*r[n,j,t]-Sum[β[k,kk]*i[k,kk,t],{k,1,ages},{kk,1,stages}]*s[j,t] -q*v[j];didt[j_,1,t_]=Sum[β[k,1]*i[k,1,t],{k,1,ages}]*Sum[r[nn,j,t],{nn,n+1-TRYm,n}]+Sum[β[k,1]*i[k,1,t],{k,1,ages}]*s[j,t] -κ[j,1]*i[j,1,t];didt[j_,2,t_]=Sum[β[k,2]*i[k,2,t],{k,1,ages}]*s[j,t] -κ[j,2]*i[j,2,t];drdt[1,j_,t_]=q*v[j]+q*Sum[κ[j,jj]*i[j,jj,t],{jj,1,stages}]-n*δ[j]*r[1,j,t];For[nn=2,nn≤n-TRYm,nn++,
drdt[nn,j_,t_]=n*δ[j]*r[nn-1,j,t]-n*δ[j]*r[nn,j,t];
];
For[nn=n+1-TRYm,nn≤n,nn++,
drdt[nn,j_,t_]=n*δ[j]*r[nn-1,j,t]-Sum[β[k,1]*i[k,1,t],{k,1,ages}]*r[nn,j,t]-n*δ[j]*r[nn,j,t];
];
pars={β[j_,1]->βV,β[j_,2]->β,δ[1]->δY,δ[2]->δO,κ[j_,jj_]->κ,v[1]->vY,v[2]->vO}; 
frac[1,1]=(1-older)(1-f);
frac[2,1]=older(1-f);
frac[1,2]=(1-older)f;
frac[2,2]=older f;
nvars=Drop[Flatten[Table[{s[j,t],Table[i[j,jj,t],{jj,1,stages}],Table[r[nn,j,t],{nn,1,n}]},{j,1,ages}]],-1];
r[n,ages,t_]=1-Total[nvars];
neqns=Drop[Flatten[Table[{D[s[j,t],t]==dsdt[j,t],Table[D[i[j,jj,t],t]==didt[j,jj,t],{jj,1,stages}],Table[D[r[nn,j,t],t]==drdt[nn,j,t],{nn,1,n}]},{j,1,ages}]],-1];
nstart=Drop[Flatten[Table[{
Table[i[j,jj,0]==frac[j,jj] i0,{jj,1,stages}],
s[j,0]==(frac[j,1]+frac[j,2])*s0,Table[r[nn,j,0]==(frac[j,1]+frac[j,2]) (r0/n),{nn,1,n}]},{j,1,ages}]]/.start,-1];
NDSolve[Flatten[{nstart/.pars,neqns/.pars}],nvars,{t,0,maxtime}]
]

all=Flatten[Table[{j,jj},{j,1,ages},{jj,1,stages}],1]

labels={"<70,Res","<70,Mutant","70+,Res","70+,Mutant"};

colours={RGBColor[0.4,0.4,0.8],ColorData["Pastel",1],RGBColor[0.4,0.6,0.4], ColorData["Pastel",3/4]}

coltab=Join[{{1->{0,colours[[1]]}}},Table[{i->{{i-1},colours[[i]]}},{i,2,Length[all]}]];
Table[{labels[[i]],colours[[i]]},{i,1,Length[all]}]//MatrixForm

No vaccinations:

TRYvY=(1-TRYolder)*0; (*Fraction of susceptibles getting vaccinated per day*) 
TRYvO=TRYolder*0; (*Fraction of older susceptibles getting vaccinated per day*)

partemp={TRYp,TRYf,TRYolder,TRYβ,TRYβV,TRYδ,TRYδO,TRYvY,TRYvO,TRYκ,TRYconvert};

solution[partemp,start];

Showing the % of the population infectious:

Show[
Plot[Evaluate@Table[ Sum[100 i[all[[c,1]],all[[c,2]],t]/.solution[partemp,start],{c,1,jcum}],{jcum,1,Length[all]}],{t,0,maxtime},PlotRange->{Automatic,{0,maxy}},
PlotStyle->None,Filling->{{1->{0,}},{2->{{1},}},{3->{{2},}},{4->{{3},}}}],
Plot[100((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black,Dashed}],
Plot[100(((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert)))-((TRYvY+TRYvO)/(TRYδ+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black}],
ImageSize->250
]

The mean infection rate for <70 (top) and 70+ individuals (bottom) in the first versus last 10 time steps:

{{Mean[Flatten[Table[(( i[1,1,t]+i[1,2,t])/(1-TRYolder))/.solution[partemp,start],{t,1,10}]]],Mean[Flatten[Table[ ((i[1,1,t]+i[1,2,t])/(1-TRYolder))/.solution[partemp,start],{t,maxtime-9,maxtime}]]]},
{Mean[Flatten[Table[ ((i[2,1,t]+i[2,2,t])/TRYolder)/.solution[partemp,start],{t,1,10}]]],Mean[Flatten[Table[ ((i[2,1,t]+i[2,2,t])/TRYolder)/.solution[partemp,start],{t,maxtime-9,maxtime}]]]}}//MatrixForm

TRYvY=(1-TRYolder)TRYvLOW; (*Fraction of younger susceptibles getting vaccinated per day*) 
TRYvO=TRYolder TRYvLOW; (*Fraction of older susceptibles getting vaccinated per day*)

partemp={TRYp,TRYf,TRYolder,TRYβ,TRYβV,TRYδ,TRYδO,TRYvY,TRYvO,TRYκ,TRYconvert};

solution[partemp,start];

Showing the % of the population infectious:

Show[
Plot[Evaluate@Table[ Sum[100 i[all[[c,1]],all[[c,2]],t]/.solution[partemp,start],{c,1,jcum}],{jcum,1,Length[all]}],{t,0,maxtime},PlotRange->{Automatic,{0,maxy}},
PlotStyle->None,Filling->{{1->{0,}},{2->{{1},}},{3->{{2},}},{4->{{3},}}}],
Plot[100((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black,Dashed}],
Plot[100(((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert)))-((TRYvY+TRYvO)/(TRYδ+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black}],
ImageSize->250
]

The mean infection rate for <70 (top) and 70+ individuals (bottom) in the first versus last 10 time steps:

{{Mean[Flatten[Table[(( i[1,1,t]+i[1,2,t])/(1-TRYolder))/.solution[partemp,start],{t,1,10}]]],Mean[Flatten[Table[ ((i[1,1,t]+i[1,2,t])/(1-TRYolder))/.solution[partemp,start],{t,maxtime-9,maxtime}]]]},
{Mean[Flatten[Table[ ((i[2,1,t]+i[2,2,t])/TRYolder)/.solution[partemp,start],{t,1,10}]]],Mean[Flatten[Table[ ((i[2,1,t]+i[2,2,t])/TRYolder)/.solution[partemp,start],{t,maxtime-9,maxtime}]]]}}//MatrixForm

TRYvY=(1-TRYolder)TRYvHIGH; (*Fraction of susceptibles getting vaccinated per day*) 
TRYvO=TRYolder TRYvHIGH; (*Fraction of older susceptibles getting vaccinated per day*)

partemp={TRYp,TRYf,TRYolder,TRYβ,TRYβV,TRYδ,TRYδO,TRYvY,TRYvO,TRYκ,TRYconvert};

solution[partemp,start];

Showing the % of the population infectious:

Show[
Plot[Evaluate@Table[ Sum[100 i[all[[c,1]],all[[c,2]],t]/.solution[partemp,start],{c,1,jcum}],{jcum,1,Length[all]}],{t,0,maxtime},PlotRange->{Automatic,{0,maxy}},
PlotStyle->None,Filling->{{1->{0,}},{2->{{1},}},{3->{{2},}},{4->{{3},}}}],
Plot[100((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black,Dashed}],
Plot[100(((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert)))-((TRYvY+TRYvO)/(TRYδ+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black}],
ImageSize->250
]

{neweq,oldeq}/.oldeq->((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert)))/.neweq->(((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert)))-((TRYvY+TRYvO)/(TRYδ+TRYκ TRYconvert)))//N

100((neweq-oldeq)/oldeq)/.oldeq->((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert)))/.neweq->(((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert)))-((TRYvY+TRYvO)/(TRYδ+TRYκ TRYconvert)))//N

The mean infection rate for <70 (top) and 70+ individuals (bottom) in the first versus last 10 time steps:

{{Mean[Flatten[Table[(( i[1,1,t]+i[1,2,t])/(1-TRYolder))/.solution[partemp,start],{t,1,10}]]],Mean[Flatten[Table[ ((i[1,1,t]+i[1,2,t])/(1-TRYolder))/.solution[partemp,start],{t,maxtime-9,maxtime}]]]},
{Mean[Flatten[Table[ ((i[2,1,t]+i[2,2,t])/TRYolder)/.solution[partemp,start],{t,1,10}]]],Mean[Flatten[Table[ ((i[2,1,t]+i[2,2,t])/TRYolder)/.solution[partemp,start],{t,maxtime-9,maxtime}]]]}}//MatrixForm

TRYvY=0; (*Fraction of susceptibles getting vaccinated per day*) 
TRYvO=TRYolder TRYvHIGH; (*Fraction of older susceptibles getting vaccinated per day*)

partemp={TRYp,TRYf,TRYolder,TRYβ,TRYβV,TRYδ,TRYδO,TRYvY,TRYvO,TRYκ,TRYconvert};

solution[partemp,start];

Showing the % of the population infectious:

Show[
Plot[Evaluate@Table[ Sum[100 i[all[[c,1]],all[[c,2]],t]/.solution[partemp,start],{c,1,jcum}],{jcum,1,Length[all]}],{t,0,maxtime},PlotRange->{Automatic,{0,maxy}},
PlotStyle->None,Filling->{{1->{0,}},{2->{{1},}},{3->{{2},}},{4->{{3},}}}],
Plot[100((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black,Dashed}],
Plot[100(((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert)))-((TRYvY+TRYvO)/(TRYδ+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black}],
ImageSize->250
]

The mean infection rate per capita for <70 (top) and 70+ individuals (bottom) in the first versus last 10 time steps:

{{Mean[Flatten[Table[(( i[1,1,t]+i[1,2,t])/(1-TRYolder))/.solution[partemp,start],{t,1,10}]]],Mean[Flatten[Table[ ((i[1,1,t]+i[1,2,t])/(1-TRYolder))/.solution[partemp,start],{t,maxtime-9,maxtime}]]]},
{Mean[Flatten[Table[ ((i[2,1,t]+i[2,2,t])/TRYolder)/.solution[partemp,start],{t,1,10}]]],Mean[Flatten[Table[ ((i[2,1,t]+i[2,2,t])/TRYolder)/.solution[partemp,start],{t,maxtime-9,maxtime}]]]}}//MatrixForm

Increasing maxtime=5000 for a more accurate estimate of the long-termimpact gives

(({{0.019987408476232836`, 0.019465007705094916`}, {0.019814068933413358`, 0.015233496283498422`}}))

so the decline in cases among <70 is 2.67% and in 70+ is 23.83%.

Clear[solution, i, s, r, ages, stages, n, neqns, nstart, nvars, didt, dsdt, drdt];

NPIs

NPIs (f=0.5, p=0.25)
Considers impact of changing NPI from no adherence (f=0) to 50% (f=0.5) instantaneously at time 0

Parameters:

maxtime=500;
maxy=8;

Parameters:

TRYp=1/4; (*Efficacy of masking per interaction*)
TRYf=1/2;(*Additional fraction masking*)
TRYolder=13/100; (*~ Age 70+ from StatsCan 5022509/38929902 *)
	(*https://www150.statcan.gc.ca/t1/tbl1/en/cv.action?pid=1710000501*)
TRYκ=1/5; (*Five day infectious period*)
TRYδ=1/125; (*Waning rate based on Menegale et al. (2023)*) 
TRYδO=1/125; (*No clear age difference noted in Menegale et al. (2023)*) 
TRYconvert=10/10; (*Only this fraction of cases seroconvert or boost immunity*)

i->((δ (β-κ))/(β (δ+q κ)))

Solve[((δ (β-κ))/(β (δ+q κ)))==inf,β]

TRYinf=1/50; (*Assumed initial endemic frequency of infected individuals*)
TRYβ=((TRYδ TRYκ)/(TRYδ-TRYinf TRYδ-TRYconvert TRYinf TRYκ));

((TRYβ)/(TRYκ))//N

Starting with no masking in the population and everybody starts masking.  Starting equilibrium would be (r0 includes all of the n resistance classes):

{s->(κ/β),i->((δ (β-κ))/(β (δ+κ)))}

start={s0->((TRYκ)/(TRYβ)),i0->((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))),r0->1-((TRYκ)/(TRYβ))-((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert)))}

{ages=2,stages=2,n=5};

Waning distribution

Show[
Plot[TRYconvert*(1-CDF[GammaDistribution[n,1/(TRYδ n)],t]),{t,0,365},PlotRange->{All,{0,1}},PlotStyle->RGBColor[0.4,0.4,0.8]],
Plot[TRYconvert*(1-CDF[GammaDistribution[n,1/(TRYδO n)],t]),{t,0,365},PlotRange->{All,{0,1}},PlotStyle->RGBColor[0.4,0.6,0.4]]
]

The two stages in this case correspond to those who do (jj=2) or do not (jj=1) engage in the protective NPI measures.

Clear[solution]
solution[{p_,f_,older_,β_,δY_,δO_,κ_,q_},start]:=solution[{p,f,older,β,δY,δO,κ,q},start]=
Block[{ages=2,stages=2,n=5},
dsdt[j_,jj_,t_]=(1-q)*κ[j,jj]*i[j,jj,t]+n*δ[j,jj]*r[n,j,jj,t]-Sum[m[jj,kk]β[jj,kk]*i[k,kk,t],{k,1,ages},{kk,1,stages}]*s[j,jj,t] ;didt[j_,jj_,t_]=Sum[m[jj,kk]β[jj,kk]*i[k,kk,t],{k,1,ages},{kk,1,stages}]*s[j,jj,t] -κ[j,jj]*i[j,jj,t];drdt[1,j_,jj_,t_]=q*κ[j,jj]*i[j,jj,t]-n*δ[j,jj]*r[1,j,jj,t];For[nn=2,nn≤n,nn++,
drdt[nn,j_,jj_,t_]=n*δ[j,jj]*r[nn-1,j,jj,t]-n*δ[j,jj]*r[nn,j,jj,t];
];
pars={m[jj_,kk_]->(((1-p))^(jj+kk-2)),β[j_,k_]->β,δ[1,jj_]->δY,δ[2,jj_]->δO,κ[j_,jj_]->κ}; 
frac[1,1]=(1-older)(1-f);
frac[2,1]=older(1-f);
frac[1,2]=(1-older)f;
frac[2,2]=older f;
nvars=Drop[Flatten[Table[{s[j,jj,t],i[j,jj,t],Table[r[nn,j,jj,t],{nn,1,n}]},{j,1,ages},{jj,1,stages}]],-1];
r[n,ages,stages,t_]=1-Total[nvars];
neqns=Drop[Flatten[Table[{D[s[j,jj,t],t]==dsdt[j,jj,t],D[i[j,jj,t],t]==didt[j,jj,t],Table[D[r[nn,j,jj,t],t]==drdt[nn,j,jj,t],{nn,1,n}]},{j,1,ages},{jj,1,stages}]],-1];
nstart=Drop[Flatten[Table[{s[j,jj,0]==frac[j,jj] s0,
i[j,jj,0]==frac[j,jj] i0,
Table[r[nn,j,jj,0]==frac[j,jj] (r0/n),{nn,1,n}]},{j,1,ages},{jj,1,stages}]]/.start,-1];
NDSolve[Flatten[{nstart/.pars,neqns/.pars}],nvars,{t,0,maxtime}]
]

all=Flatten[Table[{j,jj},{j,1,ages},{jj,1,stages}],1]

labels={"<70,None","<70,NPI","70+,None","70+,NPI"};

colours={RGBColor[0.4,0.4,0.8],ColorData["Pastel",1],RGBColor[0.4,0.6,0.4], ColorData["Pastel",3/4]}

coltab=Join[{{1->{0,colours[[1]]}}},Table[{i->{{i-1},colours[[i]]}},{i,2,Length[all]}]];
Table[{labels[[i]],colours[[i]]},{i,1,Length[all]}]//MatrixForm

Equilibrium with a heterogeneous population (allowing some fraction of seroconversion, q):

neweq[{p_,f_,older_,β_,δY_,δO_,κ_,q_}]=(1/(2 (((-1+p))^2) β (δ+q κ)))δ ((1+f (-2+p)) (-1+p) p β-(2+(-2+p) p) κ+2 (Sqrt[4 (-1+f) (-1+p) p β κ+((((-1+p) (-1+f p) β-p κ))^2)])-p (Sqrt[4 (-1+f) (-1+p) p β κ+((((-1+p) (-1+f p) β-p κ))^2)]))/.δ->δY (*Doesn't consider age differences here in waning*);

partemp={TRYp,TRYf,TRYolder,TRYβ,TRYδ,TRYδO,TRYκ,TRYconvert};

Equilibrium number of cases:

100 neweq[partemp]//N

solution[partemp,start];

plot1=Show[
Plot[Evaluate@Table[ 100 Sum[i[all[[c,1]],all[[c,2]],t]/.solution[partemp,start],{c,1,jcum}],{jcum,1,Length[all]}],{t,0,maxtime},PlotRange->{All,{0,maxy}},
PlotStyle->None,Filling->{{1->{0,}},{2->{{1},}},{3->{{2},}},{4->{{3},}}}],
Plot[100((TRYδ ((1-p)TRYβ-TRYκ))/((1-p)TRYβ (TRYδ+TRYκ TRYconvert)))/.p->0,{t,0,maxtime},PlotStyle->{Black,Dashed}],
Plot[100 neweq[partemp]/.p->TRYp/.f->TRYf/.β->TRYβ/.κ->TRYκ/.δ->TRYδ/.q->TRYconvert,{t,0,maxtime},PlotStyle->{Black}]
]

The following illustrates the relatively minor effect caused by reducing seroconversion (to q = 80%), as long as δ/q is held constant:

partemp={TRYp,TRYf,TRYolder,TRYβ,(8/10)TRYδ,(8/10)TRYδO,TRYκ,(8/10)};

Equilibrium number of cases:

100 neweq[partemp]//N

solution[partemp,start];

Showing the % of the population infectious:

plot2=Show[plot1,
Plot[Evaluate@Table[ Sum[100*i[all[[c,1]],all[[c,2]],t]/.solution[partemp,start],{c,1,jcum}],{jcum,1,Length[all]}],{t,0,maxtime},PlotRange->{Automatic,{0,maxy}},
PlotStyle->{,,,}],
Plot[100((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black,Dashed}],
Plot[100 neweq[partemp]/.p->TRYp/.f->TRYf/.β->TRYβ/.κ->TRYκ/.δ->(8/10)TRYδ/.q->(8/10)TRYconvert,{t,0,maxtime},PlotStyle->{Black}],
ImageSize->250
]

Clear[solution,i,s,r,ages,stages,n,neqns,nstart,nvars,didt,dsdt,drdt];

NPIs (f=0.5, p=0.25) - moves between classes
Allows switching between masking and non-masking (1% per day), with initial population not masking

Parameters:

maxtime=500;
maxy=8;

Parameters:

TRYp=1/4; (*Efficacy of masking per interaction*)
TRYf=1/2;(*Additional fraction masking*)
TRYolder=13/100; (*~ Age 70+ from StatsCan 5022509/38929902 *)
	(*https://www150.statcan.gc.ca/t1/tbl1/en/cv.action?pid=1710000501*)
TRYκ=1/5; (*Five day infectious period*)
TRYδ=1/125; (*Waning rate based on Menegale et al. (2023)*) 
TRYδO=1/125; (*No clear age difference noted in Menegale et al. (2023)*) 
TRYconvert=10/10; (*Only this fraction of cases seroconvert or boost immunity*)

i->((δ (β-κ))/(β (δ+q κ)))

Adding in a rate of switching classes:

TRYrate=1/100; (*Rate of switching NPI classes*)

Solve[((δ (β-κ))/(β (δ+q κ)))==inf,β]

TRYinf=1/50; (*Assumed initial endemic frequency of infected individuals*)
TRYβ=((TRYδ TRYκ)/(TRYδ-TRYinf TRYδ-TRYconvert TRYinf TRYκ));

((TRYβ)/(TRYκ))//N

Starting with no masking in the population and everybody starts masking.  Starting equilibrium would be (r0 includes all of the n resistance classes):

{s->(κ/β),i->((δ (β-κ))/(β (δ+κ)))}

start={s0->((TRYκ)/(TRYβ)),i0->((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))),r0->1-((TRYκ)/(TRYβ))-((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert)))}

{ages=2,stages=2,n=5};

Waning distribution

Show[
Plot[TRYconvert*(1-CDF[GammaDistribution[n,1/(TRYδ n)],t]),{t,0,365},PlotRange->{All,{0,1}},PlotStyle->RGBColor[0.4,0.4,0.8]],
Plot[TRYconvert*(1-CDF[GammaDistribution[n,1/(TRYδO n)],t]),{t,0,365},PlotRange->{All,{0,1}},PlotStyle->RGBColor[0.4,0.6,0.4]]
]

The two stages in this case correspond to those who do (jj=2) or do not (jj=1) engage in the protective NPI measures.

Clear[solution]
solution[{p_,f_,older_,β_,δY_,δO_,κ_,q_,rate_},start]:=solution[{p,f,older,β,δY,δO,κ,q,rate},start]=
Block[{ages=2,stages=2,n=5},
dsdt[j_,jj_,t_]=(1-q)*κ[j,jj]*i[j,jj,t]+n*δ[j,jj]*r[n,j,jj,t]-Sum[m[jj,kk]β[jj,kk]*i[k,kk,t],{k,1,ages},{kk,1,stages}]*s[j,jj,t] -
If[jj==1,f,1-f] rate s[j,jj,t]+If[jj==1,1-f,f] rate s[j,Mod[jj,2]+1,t];didt[j_,jj_,t_]=Sum[m[jj,kk]β[jj,kk]*i[k,kk,t],{k,1,ages},{kk,1,stages}]*s[j,jj,t] -κ[j,jj]*i[j,jj,t]-
If[jj==1,f,1-f] rate i[j,jj,t]+If[jj==1,1-f,f] rate i[j,Mod[jj,2]+1,t];drdt[1,j_,jj_,t_]=q*κ[j,jj]*i[j,jj,t]-n*δ[j,jj]*r[1,j,jj,t]-
If[jj==1,f,1-f] rate r[1,j,jj,t]+If[jj==1,1-f,f] rate r[1,j,Mod[jj,2]+1,t];For[nn=2,nn≤n,nn++,
drdt[nn,j_,jj_,t_]=n*δ[j,jj]*r[nn-1,j,jj,t]-n*δ[j,jj]*r[nn,j,jj,t]-
If[jj==1,f,1-f] rate r[nn,j,jj,t]+If[jj==1,1-f,f] rate r[nn,j,Mod[jj,2]+1,t];;
];
pars={m[jj_,kk_]->(((1-p))^(jj+kk-2)),β[j_,k_]->β,δ[1,jj_]->δY,δ[2,jj_]->δO,κ[j_,jj_]->κ}; 
frac[1,1]=(1-older);
frac[2,1]=older;
frac[1,2]=0;
frac[2,2]=0;
nvars=Drop[Flatten[Table[{s[j,jj,t],i[j,jj,t],Table[r[nn,j,jj,t],{nn,1,n}]},{j,1,ages},{jj,1,stages}]],-1];
r[n,ages,stages,t_]=1-Total[nvars];
neqns=Drop[Flatten[Table[{D[s[j,jj,t],t]==dsdt[j,jj,t],D[i[j,jj,t],t]==didt[j,jj,t],Table[D[r[nn,j,jj,t],t]==drdt[nn,j,jj,t],{nn,1,n}]},{j,1,ages},{jj,1,stages}]],-1];
nstart=Drop[Flatten[Table[{s[j,jj,0]==frac[j,jj] s0,
i[j,jj,0]==frac[j,jj] i0,
Table[r[nn,j,jj,0]==frac[j,jj] (r0/n),{nn,1,n}]},{j,1,ages},{jj,1,stages}]]/.start,-1];
NDSolve[Flatten[{nstart/.pars,neqns/.pars}],nvars,{t,0,maxtime}]
]

all=Flatten[Table[{j,jj},{j,1,ages},{jj,1,stages}],1]

labels={"<70,None","<70,NPI","70+,None","70+,NPI"};

colours={RGBColor[0.4,0.4,0.8],ColorData["Pastel",1],RGBColor[0.4,0.6,0.4], ColorData["Pastel",3/4]}

coltab=Join[{{1->{0,colours[[1]]}}},Table[{i->{{i-1},colours[[i]]}},{i,2,Length[all]}]];
Table[{labels[[i]],colours[[i]]},{i,1,Length[all]}]//MatrixForm

Equilibrium with a heterogeneous population (allowing some fraction of seroconversion, q):

neweq[{p_,f_,older_,β_,δY_,δO_,κ_,q_,rate_}]=(1/(2 (((-1+p))^2) β (δ+q κ)))δ ((1+f (-2+p)) (-1+p) p β-(2+(-2+p) p) κ+2 (Sqrt[4 (-1+f) (-1+p) p β κ+((((-1+p) (-1+f p) β-p κ))^2)])-p (Sqrt[4 (-1+f) (-1+p) p β κ+((((-1+p) (-1+f p) β-p κ))^2)]))/.δ->δY (*Doesn't consider age differences here in waning*);

partemp={TRYp,TRYf,TRYolder,TRYβ,TRYδ,TRYδO,TRYκ,TRYconvert,TRYrate};

100 neweq[partemp]//N

solution[partemp,start];

plot1=Show[
Plot[Evaluate@Table[ 100 Sum[i[all[[c,1]],all[[c,2]],t]/.solution[partemp,start],{c,1,jcum}],{jcum,1,Length[all]}],{t,0,maxtime},PlotRange->{All,{0,maxy}},
PlotStyle->None,Filling->{{1->{0,}},{2->{{1},}},{3->{{2},}},{4->{{3},}}}],
Plot[100((TRYδ ((1-p)TRYβ-TRYκ))/((1-p)TRYβ (TRYδ+TRYκ TRYconvert)))/.p->0,{t,0,maxtime},PlotStyle->{Black,Dashed}],
Plot[100 neweq[partemp],{t,0,maxtime},PlotStyle->{Black}]
]

The following illustrates the relatively minor effect caused by reducing seroconversion (to q = 80%), as long as δ/q is held constant:

partemp={TRYp,TRYf,TRYolder,TRYβ,(8/10)TRYδ,(8/10)TRYδO,TRYκ,(8/10),TRYrate};

100 neweq[partemp]//N

solution[partemp,start];

Showing the % of the population infectious:

plot2=Show[plot1,
Plot[Evaluate@Table[ Sum[100*i[all[[c,1]],all[[c,2]],t]/.solution[partemp,start],{c,1,jcum}],{jcum,1,Length[all]}],{t,0,maxtime},PlotRange->{Automatic,{0,maxy}},
PlotStyle->{,,,}],
Plot[100((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black,Dashed}],
Plot[100 neweq[partemp],{t,0,maxtime},PlotStyle->{Black}]
]

Clear[solution,i,s,r,ages,stages,n,neqns,nstart,nvars,didt,dsdt,drdt];

NPIs (f=0.25, p=0.5) - moves between classes (similar result since f*p is constant)
Allows switching between masking and non-masking (1% per day), with initial population not masking

Parameters:

maxtime=500;
maxy=8;

Parameters:

TRYp=1/2; (*Efficacy of masking per interaction*)
TRYf=1/4;(*Additional fraction masking*)
TRYolder=13/100; (*~ Age 70+ from StatsCan 5022509/38929902 *)
	(*https://www150.statcan.gc.ca/t1/tbl1/en/cv.action?pid=1710000501*)
TRYκ=1/5; (*Five day infectious period*)
TRYδ=1/125; (*Waning rate based on Menegale et al. (2023)*) 
TRYδO=1/125; (*No clear age difference noted in Menegale et al. (2023)*) 
TRYconvert=10/10; (*Only this fraction of cases seroconvert or boost immunity*)

i->((δ (β-κ))/(β (δ+q κ)))

Adding in a rate of switching classes:

TRYrate=1/100; (*Rate of switching NPI classes*)

Solve[((δ (β-κ))/(β (δ+q κ)))==inf,β]

TRYinf=1/50; (*Assumed initial endemic frequency of infected individuals*)
TRYβ=((TRYδ TRYκ)/(TRYδ-TRYinf TRYδ-TRYconvert TRYinf TRYκ));

((TRYβ)/(TRYκ))//N

Starting with no masking in the population and everybody starts masking.  Starting equilibrium would be (r0 includes all of the n resistance classes):

{s->(κ/β),i->((δ (β-κ))/(β (δ+κ)))}

start={s0->((TRYκ)/(TRYβ)),i0->((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))),r0->1-((TRYκ)/(TRYβ))-((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert)))}

{ages=2,stages=2,n=5};

Waning distribution

Show[
Plot[TRYconvert*(1-CDF[GammaDistribution[n,1/(TRYδ n)],t]),{t,0,365},PlotRange->{All,{0,1}},PlotStyle->RGBColor[0.4,0.4,0.8]],
Plot[TRYconvert*(1-CDF[GammaDistribution[n,1/(TRYδO n)],t]),{t,0,365},PlotRange->{All,{0,1}},PlotStyle->RGBColor[0.4,0.6,0.4]]
]

The two stages in this case correspond to those who do (jj=2) or do not (jj=1) engage in the protective NPI measures.

Clear[solution]
solution[{p_,f_,older_,β_,δY_,δO_,κ_,q_,rate_},start]:=solution[{p,f,older,β,δY,δO,κ,q,rate},start]=
Block[{ages=2,stages=2,n=5},
dsdt[j_,jj_,t_]=(1-q)*κ[j,jj]*i[j,jj,t]+n*δ[j,jj]*r[n,j,jj,t]-Sum[m[jj,kk]β[jj,kk]*i[k,kk,t],{k,1,ages},{kk,1,stages}]*s[j,jj,t] -
If[jj==1,f,1-f] rate s[j,jj,t]+If[jj==1,1-f,f] rate s[j,Mod[jj,2]+1,t];didt[j_,jj_,t_]=Sum[m[jj,kk]β[jj,kk]*i[k,kk,t],{k,1,ages},{kk,1,stages}]*s[j,jj,t] -κ[j,jj]*i[j,jj,t]-
If[jj==1,f,1-f] rate i[j,jj,t]+If[jj==1,1-f,f] rate i[j,Mod[jj,2]+1,t];drdt[1,j_,jj_,t_]=q*κ[j,jj]*i[j,jj,t]-n*δ[j,jj]*r[1,j,jj,t]-
If[jj==1,f,1-f] rate r[1,j,jj,t]+If[jj==1,1-f,f] rate r[1,j,Mod[jj,2]+1,t];For[nn=2,nn≤n,nn++,
drdt[nn,j_,jj_,t_]=n*δ[j,jj]*r[nn-1,j,jj,t]-n*δ[j,jj]*r[nn,j,jj,t]-
If[jj==1,f,1-f] rate r[nn,j,jj,t]+If[jj==1,1-f,f] rate r[nn,j,Mod[jj,2]+1,t];;
];
pars={m[jj_,kk_]->(((1-p))^(jj+kk-2)),β[j_,k_]->β,δ[1,jj_]->δY,δ[2,jj_]->δO,κ[j_,jj_]->κ}; 
frac[1,1]=(1-older);
frac[2,1]=older;
frac[1,2]=0;
frac[2,2]=0;
nvars=Drop[Flatten[Table[{s[j,jj,t],i[j,jj,t],Table[r[nn,j,jj,t],{nn,1,n}]},{j,1,ages},{jj,1,stages}]],-1];
r[n,ages,stages,t_]=1-Total[nvars];
neqns=Drop[Flatten[Table[{D[s[j,jj,t],t]==dsdt[j,jj,t],D[i[j,jj,t],t]==didt[j,jj,t],Table[D[r[nn,j,jj,t],t]==drdt[nn,j,jj,t],{nn,1,n}]},{j,1,ages},{jj,1,stages}]],-1];
nstart=Drop[Flatten[Table[{s[j,jj,0]==frac[j,jj] s0,
i[j,jj,0]==frac[j,jj] i0,
Table[r[nn,j,jj,0]==frac[j,jj] (r0/n),{nn,1,n}]},{j,1,ages},{jj,1,stages}]]/.start,-1];
NDSolve[Flatten[{nstart/.pars,neqns/.pars}],nvars,{t,0,maxtime}]
]

all=Flatten[Table[{j,jj},{j,1,ages},{jj,1,stages}],1]

labels={"<70,None","<70,NPI","70+,None","70+,NPI"};

colours={RGBColor[0.4,0.4,0.8],ColorData["Pastel",1],RGBColor[0.4,0.6,0.4], ColorData["Pastel",3/4]}

coltab=Join[{{1->{0,colours[[1]]}}},Table[{i->{{i-1},colours[[i]]}},{i,2,Length[all]}]];
Table[{labels[[i]],colours[[i]]},{i,1,Length[all]}]//MatrixForm

Equilibrium with a heterogeneous population (allowing some fraction of seroconversion, q):

neweq[{p_,f_,older_,β_,δY_,δO_,κ_,q_,rate_}]=(1/(2 (((-1+p))^2) β (δ+q κ)))δ ((1+f (-2+p)) (-1+p) p β-(2+(-2+p) p) κ+2 (Sqrt[4 (-1+f) (-1+p) p β κ+((((-1+p) (-1+f p) β-p κ))^2)])-p (Sqrt[4 (-1+f) (-1+p) p β κ+((((-1+p) (-1+f p) β-p κ))^2)]))/.δ->δY (*Doesn't consider age differences here in waning*);

partemp={TRYp,TRYf,TRYolder,TRYβ,TRYδ,TRYδO,TRYκ,TRYconvert,TRYrate};

100 neweq[partemp]//N

solution[partemp,start];

plot1=Show[
Plot[Evaluate@Table[ 100 Sum[i[all[[c,1]],all[[c,2]],t]/.solution[partemp,start],{c,1,jcum}],{jcum,1,Length[all]}],{t,0,maxtime},PlotRange->{All,{0,maxy}},
PlotStyle->None,Filling->{{1->{0,}},{2->{{1},}},{3->{{2},}},{4->{{3},}}}],
Plot[100((TRYδ ((1-p)TRYβ-TRYκ))/((1-p)TRYβ (TRYδ+TRYκ TRYconvert)))/.p->0,{t,0,maxtime},PlotStyle->{Black,Dashed}],
Plot[100 neweq[partemp],{t,0,maxtime},PlotStyle->{Black}]
]

The following illustrates the relatively minor effect caused by reducing seroconversion (to q = 80%), as long as δ/q is held constant:

partemp={TRYp,TRYf,TRYolder,TRYβ,(8/10)TRYδ,(8/10)TRYδO,TRYκ,(8/10),TRYrate};

100 neweq[partemp]//N

solution[partemp,start];

Showing the % of the population infectious:

plot2=Show[plot1,
Plot[Evaluate@Table[ Sum[100*i[all[[c,1]],all[[c,2]],t]/.solution[partemp,start],{c,1,jcum}],{jcum,1,Length[all]}],{t,0,maxtime},PlotRange->{Automatic,{0,maxy}},
PlotStyle->{,,,}],
Plot[100((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black,Dashed}],
Plot[100 neweq[partemp],{t,0,maxtime},PlotStyle->{Black}]
]

Clear[solution,i,s,r,ages,stages,n,neqns,nstart,nvars,didt,dsdt,drdt];

Alternative: (Subscript[SIR, n]) model (with leaky variants)
-> Consistent results (accounting for impact on resistant classes)

Alt Figure 3: leaky immunity & more transmissible (n=5)

maxtime=500;
maxy=8;

Parameters:

TRYp=TRYp; (*No heterogeneity, so p not relevant now.*)
TRYf=99/100;(*Starting fraction of resident*)
TRYolder=13/100; (*~ Age 70+ from StatsCan 5022509/38929902 *)
	(*https://www150.statcan.gc.ca/t1/tbl1/en/cv.action?pid=1710000501*)
TRYκ=1/5; (*Five day infectious period*)
TRYδ=1/125; (*Waning rate based on Menegale et al. (2023)*) 
TRYδO=1/125; (*No clear age difference noted in Menegale et al. (2023)*) 
TRYconvert=10/10; (*Only this fraction of cases seroconvert or boost immunity*)
TRYleak=1/10;

TRYinf=1/50; (*Assumed initial endemic frequency of infected individuals*)

Now the R[j,t] classes are no longer equal in frequency at equilibrium and must be numerically calculated:

Flatten[{0==5 δ R[5,t]-β1 I1[t] S[t],
0==-κ I1[t]+leak β1 I1[t] Sum[R[j,t],{j,1,5}]+β1 I1[t] S[t],
0==κ I1[t]-5 δ R[1,t]-leak β1 I1[t] R[1,t],
Table[0==5 δ R[-1+j,t]-5 δ R[j,t]-leak β1 I1[t] R[j,t],{j,2,5}]}]

Solve[%,Flatten[{S[t],I1[t],Table[R[j,t],{j,1,5}]}]]//Simplify

startme=%[[2]];

Finding the value of β that gives the target incidence (TRYinf):

Solve[(S[t]+I1[t]+Sum[R[j,t],{j,1,5}]/.startme/.I1[t]->TRYinf/.δ->TRYδ/.κ->TRYκ/.leak->TRYleak)==1,β1]

TRYβ=β1/.%[[2]]

(*TRYβ=((TRYδ TRYκ)/(TRYδ-TRYinf TRYδ-TRYconvert TRYinf TRYκ));*)

((TRYβ)/(TRYκ))//N

Variant has a higher transmissibility (set to have the same selection coefficient as for the immune evasive variant):

Flatten[Solve[((βV-TRYβ)/(TRYβ))TRYκ==0.0833,βV]]

TRYβV=βV/.%;

% increase in β:

100*((TRYβV-TRYβ)/(TRYβ))

and now has no immune evasive properties.  Specifically, we assume five waning classes and the variant infects only susceptibles (m=0):

TRYn=5;
TRYm=0;

Equilibrium once variant fixes:

Sort[Solve[(S[t]+I1[t]+Sum[R[j,t],{j,1,TRYn-TRYm}]/.startme/.δ->TRYδ/.κ->TRYκ/.leak->TRYleak/.β1->TRYβV)==1,I1[t]]//N]
neweq=I1[t]/.Last[%]

Starting equilibrium would be (r0 includes all of the n resistance classes and is just an approximate start):

start={s0->S[t],i0->TRYinf,r0->1-S[t]-TRYinf}/.startme/.I1[t]->TRYinf/.δ->TRYδ/.κ->TRYκ/.leak->TRYleak/.β1->TRYβ;

Waning distribution for a population of younger and older individuals:

Show[
Plot[TRYconvert*(1-CDF[GammaDistribution[TRYn,1/(TRYδ TRYn)],t]),{t,0,365},PlotRange->{All,{0,1}},PlotStyle->RGBColor[0.4,0.4,0.8]],
Plot[TRYconvert*(1-CDF[GammaDistribution[TRYn,1/(TRYδO TRYn)],t]),{t,0,365},PlotRange->{All,{0,1}},PlotStyle->RGBColor[0.4,0.6,0.4]]
]

Now the susceptible and resistant classes depend only on age (not whether they carried a variant before):
[NOTE: For the ease of colouring, the first stage class is the variant.]

{ages=2,stages=2,n=TRYn};

Clear[solution]
solution[{p_,f_,older_,β_,βV_,δY_,δO_,κ_,q_},start]:=solution[{p,f,older,β,βV,δY,δO,κ,q},start]=
Block[{ages=2,stages=2,n=TRYn},
dsdt[j_,t_]=(1-q)*Sum[κ[j,jj]*i[j,jj,t],{jj,1,stages}]+n*δ[j]*r[n,j,t]-Sum[β[k,kk]*i[k,kk,t],{k,1,ages},{kk,1,stages}]*s[j,t] ;didt[j_,1,t_]=Sum[β[k,1]*i[k,1,t],{k,1,ages}]*Sum[r[nn,j,t],{nn,n+1-TRYm,n}]+Sum[β[k,1]*i[k,1,t],{k,1,ages}]*s[j,t] +TRYleak*Sum[β[k,1]*i[k,1,t],{k,1,ages}]*Sum[r[nn,j,t],{nn,1,TRYn}]-κ[j,1]*i[j,1,t];didt[j_,2,t_]=Sum[β[k,2]*i[k,2,t],{k,1,ages}]*s[j,t]+TRYleak*Sum[β[k,2]*i[k,2,t],{k,1,ages}]*Sum[r[nn,j,t],{nn,1,TRYn}] -κ[j,2]*i[j,2,t];drdt[1,j_,t_]=q*Sum[κ[j,jj]*i[j,jj,t],{jj,1,stages}]-n*δ[j]*r[1,j,t]-TRYleak*Sum[β[k,kk]*i[k,kk,t],{k,1,ages},{kk,1,stages}]*r[1,j,t];For[nn=2,nn≤n-TRYm,nn++,
drdt[nn,j_,t_]=n*δ[j]*r[nn-1,j,t]-n*δ[j]*r[nn,j,t]-TRYleak*Sum[β[k,kk]*i[k,kk,t],{k,1,ages},{kk,1,stages}]*r[nn,j,t];
];
For[nn=n+1-TRYm,nn≤n,nn++,
drdt[nn,j_,t_]=n*δ[j]*r[nn-1,j,t]-Sum[β[k,1]*i[k,1,t],{k,1,ages}]*r[nn,j,t]-TRYleak Sum[β[k,kk]*i[k,kk,t],{k,1,ages},{kk,1,stages}]*r[nn,j,t]-n*δ[j]*r[nn,j,t];
];
pars={β[j_,1]->βV,β[j_,2]->β,δ[1]->δY,δ[2]->δO,κ[j_,jj_]->κ}; 
frac[1,1]=(1-older)(1-f);
frac[2,1]=older(1-f);
frac[1,2]=(1-older)f;
frac[2,2]=older f;
nvars=Drop[Flatten[Table[{s[j,t],Table[i[j,jj,t],{jj,1,stages}],Table[r[nn,j,t],{nn,1,n}]},{j,1,ages}]],-1];
r[n,ages,t_]=1-Total[nvars];
neqns=Drop[Flatten[Table[{D[s[j,t],t]==dsdt[j,t],Table[D[i[j,jj,t],t]==didt[j,jj,t],{jj,1,stages}],Table[D[r[nn,j,t],t]==drdt[nn,j,t],{nn,1,n}]},{j,1,ages}]],-1];
nstart=Drop[Flatten[Table[{
Table[i[j,jj,0]==frac[j,jj] i0,{jj,1,stages}],
s[j,0]==(frac[j,1]+frac[j,2])*s0,Table[r[nn,j,0]==(frac[j,1]+frac[j,2]) (r0/n),{nn,1,n}]},{j,1,ages}]]/.start,-1];
NDSolve[Flatten[{nstart/.pars,neqns/.pars}],nvars,{t,0,maxtime}]
]

all=Flatten[Table[{j,jj},{j,1,ages},{jj,1,stages}],1]

labels={"<70,Res","<70,Mutant","70+,Res","70+,Mutant"};

colours={RGBColor[0.4,0.4,0.8],ColorData["Pastel",1],RGBColor[0.4,0.6,0.4], ColorData["Pastel",3/4]}

coltab=Join[{{1->{0,colours[[1]]}}},Table[{i->{{i-1},colours[[i]]}},{i,2,Length[all]}]];
Table[{labels[[i]],colours[[i]]},{i,1,Length[all]}]//MatrixForm

partemp={TRYp,TRYf,TRYolder,TRYβ,TRYβV,TRYδ,TRYδO,TRYκ,TRYconvert};

solution[partemp,start];

Showing the % of the population infectious:

plot1=Show[
Plot[Evaluate@Table[ Sum[100 i[all[[c,1]],all[[c,2]],t]/.solution[partemp,start],{c,1,jcum}],{jcum,1,Length[all]}],{t,0,maxtime},PlotRange->{Automatic,{0,maxy}},
PlotStyle->None,Filling->{{1->{0,}},{2->{{1},}},{3->{{2},}},{4->{{3},}}}],
Plot[100TRYinf,{t,0,maxtime},PlotStyle->{Black,Dashed}],
Plot[100neweq,{t,0,maxtime},PlotStyle->{Black}],
ImageSize->250
]

Mutant fraction change between t=0 and t=50:

((Sum[i[j,1,t],{j,1,ages}])/(Sum[i[j,jj,t],{j,1,ages},{jj,1,stages}]))/.solution[partemp,start]/.t->50

Solve[Exp[s 50] ==( ((p50/(1-p50))/(p0/(1-p0)))/.p50->%/.p0->1-TRYf),s][[1]]

This is consistent with the predicted selection coefficient (updating the last part for the leaky model):

((TRYβV-TRYβ)/(TRYβ))TRYκ+Sum[R[j,t],{j,1+TRYn-TRYm,TRYn}]TRYβV/.startme/.I1[t]->TRYinf/.δ->TRYδ/.κ->TRYκ/.leak->TRYleak/.β1->TRYβV//N

New equilibrium:

neweq//N

% rise in the endemic equilibrium:

100*(((%-TRYinf))/TRYinf)

Clear[solution, i, s, r, ages, stages, n, neqns, nstart, nvars, didt, dsdt, drdt];

Alt Figure 3: leaky immunity & more permanently immune evasive (m=2, n=5)

maxtime=500;
maxy=8;

Parameters:

TRYp=TRYp; (*No heterogeneity, so p not relevant now.*)
TRYf=99/100;(*Starting fraction of resident*)
TRYolder=13/100; (*~ Age 70+ from StatsCan 5022509/38929902 *)
	(*https://www150.statcan.gc.ca/t1/tbl1/en/cv.action?pid=1710000501*)
TRYκ=1/5; (*Five day infectious period*)
TRYδ=1/125; (*Waning rate based on Menegale et al. (2023)*) 
TRYδO=1/125; (*No clear age difference noted in Menegale et al. (2023)*) 
TRYconvert=10/10; (*Only this fraction of cases seroconvert or boost immunity*)
TRYleak=1/10;

TRYinf=1/50; (*Assumed initial endemic frequency of infected individuals*)

Now the R[j,t] classes are no longer equal in frequency at equilibrium and must be numerically calculated:

Flatten[{0==5 δ R[5,t]-β1 I1[t] S[t],
0==-κ I1[t]+leak β1 I1[t] Sum[R[j,t],{j,1,5}]+β1 I1[t] S[t],
0==κ I1[t]-5 δ R[1,t]-leak β1 I1[t] R[1,t],
Table[0==5 δ R[-1+j,t]-5 δ R[j,t]-leak β1 I1[t] R[j,t],{j,2,5}]}]

Solve[%,Flatten[{S[t],I1[t],Table[R[j,t],{j,1,5}]}]]//Simplify

startme=%[[2]];

Finding the value of β that gives the target incidence (TRYinf):

Solve[(S[t]+I1[t]+Sum[R[j,t],{j,1,5}]/.startme/.I1[t]->TRYinf/.δ->TRYδ/.κ->TRYκ/.leak->TRYleak)==1,β1]

TRYβ=β1/.%[[2]]

(*TRYβ=((TRYδ TRYκ)/(TRYδ-TRYinf TRYδ-TRYconvert TRYinf TRYκ));*)

((TRYβ)/(TRYκ))//N

Variant has the same transmissibility:

TRYβV=TRYβ;

but is more immune evasive:

TRYn=5;
TRYm=2;

Equilibrium once variant fixes:

Sort[Solve[(S[t]+I1[t]+Sum[R[j,t],{j,1,TRYn-TRYm}]/.startme/.δ->TRYδ/.κ->TRYκ/.leak->TRYleak/.β1->TRYβV)==1,I1[t]]//N]
neweq=I1[t]/.Last[%]

Starting equilibrium would be (r0 includes all of the n resistance classes and is just an approximate start):

start={s0->S[t],i0->TRYinf,r0->1-S[t]-TRYinf}/.startme/.I1[t]->TRYinf/.δ->TRYδ/.κ->TRYκ/.leak->TRYleak/.β1->TRYβ;

Waning distribution for a population of younger and older individuals:

Show[
Plot[TRYconvert*(1-CDF[GammaDistribution[TRYn,1/(TRYδ TRYn)],t]),{t,0,365},PlotRange->{All,{0,1}},PlotStyle->RGBColor[0.4,0.4,0.8]],
Plot[TRYconvert*(1-CDF[GammaDistribution[TRYn,1/(TRYδO TRYn)],t]),{t,0,365},PlotRange->{All,{0,1}},PlotStyle->RGBColor[0.4,0.6,0.4]]
]

Now the susceptible and resistant classes depend only on age (not whether they carried a variant before):
[NOTE: For the ease of colouring, the first stage class is the variant.]

{ages=2,stages=2,n=TRYn};

Clear[solution]
solution[{p_,f_,older_,β_,βV_,δY_,δO_,κ_,q_},start]:=solution[{p,f,older,β,βV,δY,δO,κ,q},start]=
Block[{ages=2,stages=2,n=TRYn},
dsdt[j_,t_]=(1-q)*Sum[κ[j,jj]*i[j,jj,t],{jj,1,stages}]+n*δ[j]*r[n,j,t]-Sum[β[k,kk]*i[k,kk,t],{k,1,ages},{kk,1,stages}]*s[j,t] ;didt[j_,1,t_]=Sum[β[k,1]*i[k,1,t],{k,1,ages}]*Sum[r[nn,j,t],{nn,n+1-TRYm,n}]+Sum[β[k,1]*i[k,1,t],{k,1,ages}]*s[j,t] +TRYleak*Sum[β[k,1]*i[k,1,t],{k,1,ages}]*Sum[r[nn,j,t],{nn,1,TRYn}]-κ[j,1]*i[j,1,t];didt[j_,2,t_]=Sum[β[k,2]*i[k,2,t],{k,1,ages}]*s[j,t]+TRYleak*Sum[β[k,2]*i[k,2,t],{k,1,ages}]*Sum[r[nn,j,t],{nn,1,TRYn}] -κ[j,2]*i[j,2,t];drdt[1,j_,t_]=q*Sum[κ[j,jj]*i[j,jj,t],{jj,1,stages}]-n*δ[j]*r[1,j,t]-TRYleak*Sum[β[k,kk]*i[k,kk,t],{k,1,ages},{kk,1,stages}]*r[1,j,t];For[nn=2,nn≤n-TRYm,nn++,
drdt[nn,j_,t_]=n*δ[j]*r[nn-1,j,t]-n*δ[j]*r[nn,j,t]-TRYleak*Sum[β[k,kk]*i[k,kk,t],{k,1,ages},{kk,1,stages}]*r[nn,j,t];
];
For[nn=n+1-TRYm,nn≤n,nn++,
drdt[nn,j_,t_]=n*δ[j]*r[nn-1,j,t]-Sum[β[k,1]*i[k,1,t],{k,1,ages}]*r[nn,j,t]-TRYleak Sum[β[k,kk]*i[k,kk,t],{k,1,ages},{kk,1,stages}]*r[nn,j,t]-n*δ[j]*r[nn,j,t];
];
pars={β[j_,1]->βV,β[j_,2]->β,δ[1]->δY,δ[2]->δO,κ[j_,jj_]->κ}; 
frac[1,1]=(1-older)(1-f);
frac[2,1]=older(1-f);
frac[1,2]=(1-older)f;
frac[2,2]=older f;
nvars=Drop[Flatten[Table[{s[j,t],Table[i[j,jj,t],{jj,1,stages}],Table[r[nn,j,t],{nn,1,n}]},{j,1,ages}]],-1];
r[n,ages,t_]=1-Total[nvars];
neqns=Drop[Flatten[Table[{D[s[j,t],t]==dsdt[j,t],Table[D[i[j,jj,t],t]==didt[j,jj,t],{jj,1,stages}],Table[D[r[nn,j,t],t]==drdt[nn,j,t],{nn,1,n}]},{j,1,ages}]],-1];
nstart=Drop[Flatten[Table[{
Table[i[j,jj,0]==frac[j,jj] i0,{jj,1,stages}],
s[j,0]==(frac[j,1]+frac[j,2])*s0,Table[r[nn,j,0]==(frac[j,1]+frac[j,2]) (r0/n),{nn,1,n}]},{j,1,ages}]]/.start,-1];
NDSolve[Flatten[{nstart/.pars,neqns/.pars}],nvars,{t,0,maxtime}]
]

all=Flatten[Table[{j,jj},{j,1,ages},{jj,1,stages}],1]

labels={"<70,Res","<70,Mutant","70+,Res","70+,Mutant"};

colours={RGBColor[0.4,0.4,0.8],ColorData["Pastel",1],RGBColor[0.4,0.6,0.4], ColorData["Pastel",3/4]}

coltab=Join[{{1->{0,colours[[1]]}}},Table[{i->{{i-1},colours[[i]]}},{i,2,Length[all]}]];
Table[{labels[[i]],colours[[i]]},{i,1,Length[all]}]//MatrixForm

partemp={TRYp,TRYf,TRYolder,TRYβ,TRYβV,TRYδ,TRYδO,TRYκ,TRYconvert};

solution[partemp,start];

Showing the % of the population infectious:

plot1=Show[
Plot[Evaluate@Table[ Sum[100 i[all[[c,1]],all[[c,2]],t]/.solution[partemp,start],{c,1,jcum}],{jcum,1,Length[all]}],{t,0,maxtime},PlotRange->{Automatic,{0,maxy}},
PlotStyle->None,Filling->{{1->{0,}},{2->{{1},}},{3->{{2},}},{4->{{3},}}}],
Plot[100TRYinf,{t,0,maxtime},PlotStyle->{Black,Dashed}],
Plot[100neweq,{t,0,maxtime},PlotStyle->{Black}],
ImageSize->250
]

Mutant fraction change between t=0 and t=50:

((Sum[i[j,1,t],{j,1,ages}])/(Sum[i[j,jj,t],{j,1,ages},{jj,1,stages}]))/.solution[partemp,start]/.t->50

Solve[Exp[s 50] ==( ((p50/(1-p50))/(p0/(1-p0)))/.p50->%/.p0->1-TRYf),s][[1]]

This is consistent with the predicted selection coefficient (updating the last part for the leaky model):

((TRYβV-TRYβ)/(TRYβ))TRYκ+Sum[R[j,t],{j,1+TRYn-TRYm,TRYn}]TRYβV/.startme/.I1[t]->TRYinf/.δ->TRYδ/.κ->TRYκ/.leak->TRYleak/.β1->TRYβV//N

New equilibrium:

neweq//N

% rise in the endemic equilibrium:

100*(((%-TRYinf))/TRYinf)

Clear[solution, i, s, r, ages, stages, n, neqns, nstart, nvars, didt, dsdt, drdt];

Alternative: SIR model (with leaky variants)
-> Consistent results (accounting for impact on resistant classes)

Now variant causes resistant class to be infectable (causing the recovered class to be leaky)

Figure S2A: more transmissible (equivalent to n=1 explored above)

maxtime=500;
maxy=8;

Parameters:

TRYp=TRYp; (*No heterogeneity, so p not relevant now.*)
TRYf=99/100;(*Starting fraction of resident*)
TRYolder=13/100; (*~ Age 70+ from StatsCan 5022509/38929902 *)
	(*https://www150.statcan.gc.ca/t1/tbl1/en/cv.action?pid=1710000501*)
TRYκ=1/5; (*Five day infectious period*)
TRYδ=1/125; (*Waning rate based on Menegale et al. (2023)*) 
TRYδO=1/125; (*No clear age difference noted in Menegale et al. (2023)*) 
TRYconvert=10/10; (*Only this fraction of cases seroconvert or boost immunity*)
TRYleak=0;

i->((δ (β-κ))/(β (δ+q κ)))

Solve[((δ (β-κ))/(β (δ+q κ)))==inf,β]

TRYinf=1/50; (*Assumed initial endemic frequency of infected individuals*)
TRYβ=((TRYδ TRYκ)/(TRYδ-TRYinf TRYδ-TRYconvert TRYinf TRYκ));

((TRYβ)/(TRYκ))//N

Variant has a higher transmissibility (set to have the same selection coefficient as for the immune evasive variant):

Flatten[Solve[((βV-TRYβ)/(TRYβ))TRYκ==0.0833,βV]]

TRYβV=βV/.%;

% increase in β:

100*((TRYβV-TRYβ)/(TRYβ))

and now has no immune evasive properties (neither the resident nor the variant can infect individuals in the RL class).

New equilibrium is the same as in the (Subscript[SIR, n]) model when there is no leakiness:

Limit[((leak β2-δ-κ+(Sqrt[(leak^2) ((β2)^2)+2 leak β2 (δ-κ)-4 leak δ κ+(((δ+κ))^2)]))/(2 leak β2))/.β2->TRYβV/.δ->TRYδ/.κ->TRYκ,leak->0]

(((TRYδ(TRYβV-TRYκ))/(TRYβV (TRYδ+TRYκ TRYconvert))))

Starting equilibrium would be (r0 includes all of the resistance classes):
{s->(κ/β),i->((δ (β-κ))/(β (δ+q κ)))}

start={s0->((TRYκ)/(TRYβ)),i0->((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))),r0->1-((TRYκ)/(TRYβ))-((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert)))}

[NOTE: For the ease of colouring, the first stage class is the variant.]

{ages=2,stages=2};

Clear[solution]
solution[{p_,f_,older_,β_,βV_,δY_,δO_,κ_,q_},start]:=solution[{p,f,older,β,βV,δY,δO,κ,q},start]=
Block[{ages=2,stages=2},
dsdt[j_,t_]=(1-q)*Sum[κ[j,jj]*i[j,jj,t],{jj,1,stages}]+δ[j]*r[j,t]-Sum[β[k,kk]*i[k,kk,t],{k,1,ages},{kk,1,stages}]*s[j,t] ;didt[j_,1,t_]=TRYleak Sum[β[k,1]*i[k,1,t],{k,1,ages}]*r[j,t]+
Sum[β[k,1]*i[k,1,t],{k,1,ages}]*s[j,t] -κ[j,1]*i[j,1,t];didt[j_,2,t_]=Sum[β[k,2]*i[k,2,t],{k,1,ages}]*s[j,t] -κ[j,2]*i[j,2,t];drdt[j_,t_]=q*Sum[κ[j,jj]*i[j,jj,t],{jj,1,stages}]-δ[j]*r[j,t]-TRYleak Sum[β[k,1]*i[k,1,t],{k,1,ages}]*r[j,t];pars={β[j_,1]->βV,β[j_,2]->β,δ[1]->δY,δ[2]->δO,κ[j_,jj_]->κ}; 
frac[1,1]=(1-older)(1-f);
frac[2,1]=older(1-f);
frac[1,2]=(1-older)f;
frac[2,2]=older f;
nvars=Drop[Flatten[Table[{s[j,t],Table[i[j,jj,t],{jj,1,stages}],r[j,t]},{j,1,ages}]],-1];
r[ages,t_]=1-Total[nvars];
neqns=Drop[Flatten[Table[{D[s[j,t],t]==dsdt[j,t],Table[D[i[j,jj,t],t]==didt[j,jj,t],{jj,1,stages}],D[r[j,t],t]==drdt[j,t]},{j,1,ages}]],-1];
nstart=Drop[Flatten[Table[{
Table[i[j,jj,0]==frac[j,jj] i0,{jj,1,stages}],
s[j,0]==(frac[j,1]+frac[j,2])*s0,r[j,0]==(frac[j,1]+frac[j,2]) r0},{j,1,ages}]]/.start,-1];
NDSolve[Flatten[{nstart/.pars,neqns/.pars}],nvars,{t,0,maxtime}]
]

all=Flatten[Table[{j,jj},{j,1,ages},{jj,1,stages}],1]

labels={"<70,Res","<70,Mutant","70+,Res","70+,Mutant"};

colours={RGBColor[0.4,0.4,0.8],ColorData["Pastel",1],RGBColor[0.4,0.6,0.4], ColorData["Pastel",3/4]}

coltab=Join[{{1->{0,colours[[1]]}}},Table[{i->{{i-1},colours[[i]]}},{i,2,Length[all]}]];
Table[{labels[[i]],colours[[i]]},{i,1,Length[all]}]//MatrixForm

partemp={TRYp,TRYf,TRYolder,TRYβ,TRYβV,TRYδ,TRYδO,TRYκ,TRYconvert};

solution[partemp,start];

Showing the % of the population infectious:

plot1 =Show[
Plot[Evaluate@Table[ Sum[100*i[all[[c,1]],all[[c,2]],t]/.solution[partemp,start],{c,1,jcum}],{jcum,1,Length[all]}],{t,0,maxtime},PlotRange->{Automatic,{0,maxy}},
PlotStyle->None,Filling->{{1->{0,}},{2->{{1},}},{3->{{2},}},{4->{{3},}}}],
Plot[100((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black,Dashed}],
Plot[100((TRYδ (TRYβV-TRYκ))/(TRYβV (TRYδ+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black}],
ImageSize->250
]

Mutant fraction change between t=0 and t=50:

((Sum[i[j,1,t],{j,1,ages}])/(Sum[i[j,jj,t],{j,1,ages},{jj,1,stages}]))/.solution[partemp,start]/.t->50

Solve[Exp[s 50] ==( ((p50/(1-p50))/(p0/(1-p0)))/.p50->%/.p0->1-TRYf),s][[1]]

This is consistent with the predicted selection coefficient:

((TRYβV-TRYβ)/(TRYβ))TRYκ+TRYleak(((TRYκ TRYconvert(TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))))TRYβV//N

New equilibrium:

(((TRYδ(TRYβV-TRYκ))/(TRYβV (TRYδ+TRYκ TRYconvert))))//N

% rise in the endemic equilibrium:

100*(((%-TRYinf))/TRYinf)

Clear[solution, i, s, r, ages, stages, n, neqns, nstart, nvars, didt, dsdt, drdt];

Figure S3A: more transmissible & immune evasive (variant causes leakiness: R class 10% as infective as S class)

maxtime=500;
maxy=12;

Parameters:

TRYp=TRYp; (*No heterogeneity, so p not relevant now.*)
TRYf=99/100;(*Starting fraction of resident*)
TRYolder=13/100; (*~ Age 70+ from StatsCan 5022509/38929902 *)
	(*https://www150.statcan.gc.ca/t1/tbl1/en/cv.action?pid=1710000501*)
TRYκ=1/5; (*Five day infectious period*)
TRYδ=1/125; (*Waning rate based on Menegale et al. (2023)*) 
TRYδO=1/125; (*No clear age difference noted in Menegale et al. (2023)*) 
TRYconvert=10/10; (*Only this fraction of cases seroconvert or boost immunity*)
TRYleak=1/10;

i->((δ (β-κ))/(β (δ+q κ)))

Solve[((δ (β-κ))/(β (δ+q κ)))==inf,β]

TRYinf=1/50; (*Assumed initial endemic frequency of infected individuals*)
TRYβ=((TRYδ TRYκ)/(TRYδ-TRYinf TRYδ-TRYconvert TRYinf TRYκ));

((TRYβ)/(TRYκ))//N

Variant has a higher transmissibility (set to have the same selection coefficient as for the immune evasive variant):

Flatten[Solve[TRYleak(1-((TRYκ)/(TRYβ))-((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))))βV+((βV-TRYβ)/(TRYβ))TRYκ==0.0833,βV]]

TRYβV=βV/.%;

% increase in β:

100*((TRYβV-TRYβ)/(TRYβ))

and now has some immune evasive properties (through leaky immunity)

New equilibrium:

neweq=((leak β2-δ-κ+(Sqrt[(leak^2) ((β2)^2)+2 leak β2 (δ-κ)-4 leak δ κ+(((δ+κ))^2)]))/(2 leak β2))/.β2->TRYβV/.δ->TRYδ/.κ->TRYκ/.leak->TRYleak

Starting equilibrium would be (r0 includes all of the resistance classes):
{s->(κ/β),i->((δ (β-κ))/(β (δ+q κ)))}

start={s0->((TRYκ)/(TRYβ)),i0->((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))),r0->1-((TRYκ)/(TRYβ))-((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert)))}

[NOTE: For the ease of colouring, the first stage class is the variant.]

{ages=2,stages=2};

Clear[solution]
solution[{p_,f_,older_,β_,βV_,δY_,δO_,κ_,q_},start]:=solution[{p,f,older,β,βV,δY,δO,κ,q},start]=
Block[{ages=2,stages=2},
dsdt[j_,t_]=(1-q)*Sum[κ[j,jj]*i[j,jj,t],{jj,1,stages}]+δ[j]*r[j,t]-Sum[β[k,kk]*i[k,kk,t],{k,1,ages},{kk,1,stages}]*s[j,t] ;didt[j_,1,t_]=TRYleak Sum[β[k,1]*i[k,1,t],{k,1,ages}]*r[j,t]+
Sum[β[k,1]*i[k,1,t],{k,1,ages}]*s[j,t] -κ[j,1]*i[j,1,t];didt[j_,2,t_]=Sum[β[k,2]*i[k,2,t],{k,1,ages}]*s[j,t] -κ[j,2]*i[j,2,t];drdt[j_,t_]=q*Sum[κ[j,jj]*i[j,jj,t],{jj,1,stages}]-δ[j]*r[j,t]-TRYleak Sum[β[k,1]*i[k,1,t],{k,1,ages}]*r[j,t];pars={β[j_,1]->βV,β[j_,2]->β,δ[1]->δY,δ[2]->δO,κ[j_,jj_]->κ}; 
frac[1,1]=(1-older)(1-f);
frac[2,1]=older(1-f);
frac[1,2]=(1-older)f;
frac[2,2]=older f;
nvars=Drop[Flatten[Table[{s[j,t],Table[i[j,jj,t],{jj,1,stages}],r[j,t]},{j,1,ages}]],-1];
r[ages,t_]=1-Total[nvars];
neqns=Drop[Flatten[Table[{D[s[j,t],t]==dsdt[j,t],Table[D[i[j,jj,t],t]==didt[j,jj,t],{jj,1,stages}],D[r[j,t],t]==drdt[j,t]},{j,1,ages}]],-1];
nstart=Drop[Flatten[Table[{
Table[i[j,jj,0]==frac[j,jj] i0,{jj,1,stages}],
s[j,0]==(frac[j,1]+frac[j,2])*s0,r[j,0]==(frac[j,1]+frac[j,2]) r0},{j,1,ages}]]/.start,-1];
NDSolve[Flatten[{nstart/.pars,neqns/.pars}],nvars,{t,0,maxtime}]
]

all=Flatten[Table[{j,jj},{j,1,ages},{jj,1,stages}],1]

labels={"<70,Res","<70,Mutant","70+,Res","70+,Mutant"};

colours={RGBColor[0.4,0.4,0.8],ColorData["Pastel",1],RGBColor[0.4,0.6,0.4], ColorData["Pastel",3/4]}

coltab=Join[{{1->{0,colours[[1]]}}},Table[{i->{{i-1},colours[[i]]}},{i,2,Length[all]}]];
Table[{labels[[i]],colours[[i]]},{i,1,Length[all]}]//MatrixForm

partemp={TRYp,TRYf,TRYolder,TRYβ,TRYβV,TRYδ,TRYδO,TRYκ,TRYconvert};

solution[partemp,start];

Showing the % of the population infectious:

plot1 =Show[
Plot[Evaluate@Table[ Sum[100*i[all[[c,1]],all[[c,2]],t]/.solution[partemp,start],{c,1,jcum}],{jcum,1,Length[all]}],{t,0,maxtime},PlotRange->{Automatic,{0,maxy}},
PlotStyle->None,Filling->{{1->{0,}},{2->{{1},}},{3->{{2},}},{4->{{3},}}}],
Plot[100((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black,Dashed}],
Plot[100neweq,{t,0,maxtime},PlotStyle->{Black}],
ImageSize->250
]

Mutant fraction change between t=0 and t=50:

((Sum[i[j,1,t],{j,1,ages}])/(Sum[i[j,jj,t],{j,1,ages},{jj,1,stages}]))/.solution[partemp,start]/.t->50

Solve[Exp[s 50] ==( ((p50/(1-p50))/(p0/(1-p0)))/.p50->%/.p0->1-TRYf),s][[1]]

This is consistent with the predicted selection coefficient:

((TRYβV-TRYβ)/(TRYβ))TRYκ+TRYleak(((TRYκ TRYconvert(TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))))TRYβV//N

% rise in the endemic equilibrium:

100*(((neweq-TRYinf))/TRYinf)

Clear[solution, i, s, r, ages, stages, n, neqns, nstart, nvars, didt, dsdt, drdt];

Figure S3B: more transmissible & immune evasive (variant causes leakiness: R class 20% as infective as S class)

maxtime=500;
maxy=12;

Parameters:

TRYp=TRYp; (*No heterogeneity, so p not relevant now.*)
TRYf=99/100;(*Starting fraction of resident*)
TRYolder=13/100; (*~ Age 70+ from StatsCan 5022509/38929902 *)
	(*https://www150.statcan.gc.ca/t1/tbl1/en/cv.action?pid=1710000501*)
TRYκ=1/5; (*Five day infectious period*)
TRYδ=1/125; (*Waning rate based on Menegale et al. (2023)*) 
TRYδO=1/125; (*No clear age difference noted in Menegale et al. (2023)*) 
TRYconvert=10/10; (*Only this fraction of cases seroconvert or boost immunity*)
TRYleak=1/5;

i->((δ (β-κ))/(β (δ+q κ)))

Solve[((δ (β-κ))/(β (δ+q κ)))==inf,β]

TRYinf=1/50; (*Assumed initial endemic frequency of infected individuals*)
TRYβ=((TRYδ TRYκ)/(TRYδ-TRYinf TRYδ-TRYconvert TRYinf TRYκ));

((TRYβ)/(TRYκ))//N

Variant has a higher transmissibility (set to have the same selection coefficient as for the immune evasive variant):

Flatten[Solve[TRYleak(1-((TRYκ)/(TRYβ))-((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))))βV+((βV-TRYβ)/(TRYβ))TRYκ==0.0833,βV]]

TRYβV=βV/.%;

% increase in β:

100*((TRYβV-TRYβ)/(TRYβ))

and now has some immune evasive properties (through leaky immunity)

New equilibrium:

neweq=((leak β2-δ-κ+(Sqrt[(leak^2) ((β2)^2)+2 leak β2 (δ-κ)-4 leak δ κ+(((δ+κ))^2)]))/(2 leak β2))/.β2->TRYβV/.δ->TRYδ/.κ->TRYκ/.leak->TRYleak

Starting equilibrium would be (r0 includes all of the resistance classes):
{s->(κ/β),i->((δ (β-κ))/(β (δ+q κ)))}

start={s0->((TRYκ)/(TRYβ)),i0->((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))),r0->1-((TRYκ)/(TRYβ))-((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert)))}

[NOTE: For the ease of colouring, the first stage class is the variant.]

{ages=2,stages=2};

Clear[solution]
solution[{p_,f_,older_,β_,βV_,δY_,δO_,κ_,q_},start]:=solution[{p,f,older,β,βV,δY,δO,κ,q},start]=
Block[{ages=2,stages=2},
dsdt[j_,t_]=(1-q)*Sum[κ[j,jj]*i[j,jj,t],{jj,1,stages}]+δ[j]*r[j,t]-Sum[β[k,kk]*i[k,kk,t],{k,1,ages},{kk,1,stages}]*s[j,t] ;didt[j_,1,t_]=TRYleak Sum[β[k,1]*i[k,1,t],{k,1,ages}]*r[j,t]+
Sum[β[k,1]*i[k,1,t],{k,1,ages}]*s[j,t] -κ[j,1]*i[j,1,t];didt[j_,2,t_]=Sum[β[k,2]*i[k,2,t],{k,1,ages}]*s[j,t] -κ[j,2]*i[j,2,t];drdt[j_,t_]=q*Sum[κ[j,jj]*i[j,jj,t],{jj,1,stages}]-δ[j]*r[j,t]-TRYleak Sum[β[k,1]*i[k,1,t],{k,1,ages}]*r[j,t];pars={β[j_,1]->βV,β[j_,2]->β,δ[1]->δY,δ[2]->δO,κ[j_,jj_]->κ}; 
frac[1,1]=(1-older)(1-f);
frac[2,1]=older(1-f);
frac[1,2]=(1-older)f;
frac[2,2]=older f;
nvars=Drop[Flatten[Table[{s[j,t],Table[i[j,jj,t],{jj,1,stages}],r[j,t]},{j,1,ages}]],-1];
r[ages,t_]=1-Total[nvars];
neqns=Drop[Flatten[Table[{D[s[j,t],t]==dsdt[j,t],Table[D[i[j,jj,t],t]==didt[j,jj,t],{jj,1,stages}],D[r[j,t],t]==drdt[j,t]},{j,1,ages}]],-1];
nstart=Drop[Flatten[Table[{
Table[i[j,jj,0]==frac[j,jj] i0,{jj,1,stages}],
s[j,0]==(frac[j,1]+frac[j,2])*s0,r[j,0]==(frac[j,1]+frac[j,2]) r0},{j,1,ages}]]/.start,-1];
NDSolve[Flatten[{nstart/.pars,neqns/.pars}],nvars,{t,0,maxtime}]
]

all=Flatten[Table[{j,jj},{j,1,ages},{jj,1,stages}],1]

labels={"<70,Res","<70,Mutant","70+,Res","70+,Mutant"};

colours={RGBColor[0.4,0.4,0.8],ColorData["Pastel",1],RGBColor[0.4,0.6,0.4], ColorData["Pastel",3/4]}

coltab=Join[{{1->{0,colours[[1]]}}},Table[{i->{{i-1},colours[[i]]}},{i,2,Length[all]}]];
Table[{labels[[i]],colours[[i]]},{i,1,Length[all]}]//MatrixForm

partemp={TRYp,TRYf,TRYolder,TRYβ,TRYβV,TRYδ,TRYδO,TRYκ,TRYconvert};

solution[partemp,start];

Showing the % of the population infectious:

plot1 =Show[
Plot[Evaluate@Table[ Sum[100*i[all[[c,1]],all[[c,2]],t]/.solution[partemp,start],{c,1,jcum}],{jcum,1,Length[all]}],{t,0,maxtime},PlotRange->{Automatic,{0,maxy}},
PlotStyle->None,Filling->{{1->{0,}},{2->{{1},}},{3->{{2},}},{4->{{3},}}}],
Plot[100((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black,Dashed}],
Plot[100neweq,{t,0,maxtime},PlotStyle->{Black}],
ImageSize->250
]

Mutant fraction change between t=0 and t=50:

((Sum[i[j,1,t],{j,1,ages}])/(Sum[i[j,jj,t],{j,1,ages},{jj,1,stages}]))/.solution[partemp,start]/.t->50

Solve[Exp[s 50] ==( ((p50/(1-p50))/(p0/(1-p0)))/.p50->%/.p0->1-TRYf),s][[1]]

This is consistent with the predicted selection coefficient:

((TRYβV-TRYβ)/(TRYβ))TRYκ+TRYleak(((TRYκ TRYconvert(TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))))TRYβV//N

% rise in the endemic equilibrium:

100*(((neweq-TRYinf))/TRYinf)

Clear[solution, i, s, r, ages, stages, n, neqns, nstart, nvars, didt, dsdt, drdt];

Figure S3C: persistently immune evasive (variant causes leakiness: R class 40% as infective as S class)

maxtime=500;
maxy=12;

Parameters:

TRYp=TRYp; (*No heterogeneity, so p not relevant now.*)
TRYf=99/100;(*Starting fraction of resident*)
TRYolder=13/100; (*~ Age 70+ from StatsCan 5022509/38929902 *)
	(*https://www150.statcan.gc.ca/t1/tbl1/en/cv.action?pid=1710000501*)
TRYκ=1/5; (*Five day infectious period*)
TRYδ=1/125; (*Waning rate based on Menegale et al. (2023)*) 
TRYδO=1/125; (*No clear age difference noted in Menegale et al. (2023)*) 
TRYconvert=10/10; (*Only this fraction of cases seroconvert or boost immunity*)
TRYleak=2/5;

i->((δ (β-κ))/(β (δ+q κ)))

Solve[((δ (β-κ))/(β (δ+q κ)))==inf,β]

TRYinf=1/50; (*Assumed initial endemic frequency of infected individuals*)
TRYβ=((TRYδ TRYκ)/(TRYδ-TRYinf TRYδ-TRYconvert TRYinf TRYκ));

((TRYβ)/(TRYκ))//N

Variant has the same transmissibility:

TRYβV=TRYβ;

but infects earlier during waning by allowing leaky immunity.

New equilibrium:

neweq=((leak β2-δ-κ+(Sqrt[(leak^2) ((β2)^2)+2 leak β2 (δ-κ)-4 leak δ κ+(((δ+κ))^2)]))/(2 leak β2))/.leak->TRYleak/.β2->TRYβV/.δ->TRYδ/.κ->TRYκ//N

Starting equilibrium would be (r0 includes all of the resistance classes):
{s->(κ/β),i->((δ (β-κ))/(β (δ+q κ)))}

start={s0->((TRYκ)/(TRYβ)),i0->((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))),r0->1-((TRYκ)/(TRYβ))-((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert)))}

[NOTE: For the ease of colouring, the first stage class is the variant.]

{ages=2,stages=2};

Clear[solution]
solution[{p_,f_,older_,β_,βV_,δY_,δO_,κ_,q_},start]:=solution[{p,f,older,β,βV,δY,δO,κ,q},start]=
Block[{ages=2,stages=2},
dsdt[j_,t_]=(1-q)*Sum[κ[j,jj]*i[j,jj,t],{jj,1,stages}]+δ[j]*r[j,t]-Sum[β[k,kk]*i[k,kk,t],{k,1,ages},{kk,1,stages}]*s[j,t] ;didt[j_,1,t_]=TRYleak Sum[β[k,1]*i[k,1,t],{k,1,ages}]*r[j,t]+
Sum[β[k,1]*i[k,1,t],{k,1,ages}]*s[j,t] -κ[j,1]*i[j,1,t];didt[j_,2,t_]=Sum[β[k,2]*i[k,2,t],{k,1,ages}]*s[j,t] -κ[j,2]*i[j,2,t];drdt[j_,t_]=q*Sum[κ[j,jj]*i[j,jj,t],{jj,1,stages}]-δ[j]*r[j,t]-TRYleak Sum[β[k,1]*i[k,1,t],{k,1,ages}]*r[j,t];pars={β[j_,1]->βV,β[j_,2]->β,δ[1]->δY,δ[2]->δO,κ[j_,jj_]->κ}; 
frac[1,1]=(1-older)(1-f);
frac[2,1]=older(1-f);
frac[1,2]=(1-older)f;
frac[2,2]=older f;
nvars=Drop[Flatten[Table[{s[j,t],Table[i[j,jj,t],{jj,1,stages}],r[j,t]},{j,1,ages}]],-1];
r[ages,t_]=1-Total[nvars];
neqns=Drop[Flatten[Table[{D[s[j,t],t]==dsdt[j,t],Table[D[i[j,jj,t],t]==didt[j,jj,t],{jj,1,stages}],D[r[j,t],t]==drdt[j,t]},{j,1,ages}]],-1];
nstart=Drop[Flatten[Table[{
Table[i[j,jj,0]==frac[j,jj] i0,{jj,1,stages}],
s[j,0]==(frac[j,1]+frac[j,2])*s0,r[j,0]==(frac[j,1]+frac[j,2]) r0},{j,1,ages}]]/.start,-1];
NDSolve[Flatten[{nstart/.pars,neqns/.pars}],nvars,{t,0,maxtime}]
]

all=Flatten[Table[{j,jj},{j,1,ages},{jj,1,stages}],1]

labels={"<70,Res","<70,Mutant","70+,Res","70+,Mutant"};

colours={RGBColor[0.4,0.4,0.8],ColorData["Pastel",1],RGBColor[0.4,0.6,0.4], ColorData["Pastel",3/4]}

coltab=Join[{{1->{0,colours[[1]]}}},Table[{i->{{i-1},colours[[i]]}},{i,2,Length[all]}]];
Table[{labels[[i]],colours[[i]]},{i,1,Length[all]}]//MatrixForm

partemp={TRYp,TRYf,TRYolder,TRYβ,TRYβV,TRYδ,TRYδO,TRYκ,TRYconvert};

solution[partemp,start];

Showing the % of the population infectious:

plot1 =Show[
Plot[Evaluate@Table[ Sum[100*i[all[[c,1]],all[[c,2]],t]/.solution[partemp,start],{c,1,jcum}],{jcum,1,Length[all]}],{t,0,maxtime},PlotRange->{Automatic,{0,maxy}},
PlotStyle->None,Filling->{{1->{0,}},{2->{{1},}},{3->{{2},}},{4->{{3},}}}],
Plot[100((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black,Dashed}],
Plot[100neweq,{t,0,maxtime},PlotStyle->{Black}],
ImageSize->250
]

Mutant fraction change between t=0 and t=50:

((Sum[i[j,1,t],{j,1,ages}])/(Sum[i[j,jj,t],{j,1,ages},{jj,1,stages}]))/.solution[partemp,start]/.t->50

Solve[Exp[s 50] ==( ((p50/(1-p50))/(p0/(1-p0)))/.p50->%/.p0->1-TRYf),s][[1]]

This is consistent with the predicted selection coefficient:

((TRYβV-TRYβ)/(TRYβ))TRYκ+TRYleak(((TRYκ TRYconvert(TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))))TRYβV//N

% rise in the endemic equilibrium:

100*(((neweq-TRYinf))/TRYinf)

The following illustrates the relatively minor effect caused by reducing seroconversion (to q = 80%), as long as δ/q is held constant:

partemp={TRYp,TRYf,TRYolder,TRYβ,TRYβV,(8/10)TRYδ,(8/10)TRYδO,TRYκ,(8/10)};

solution[partemp,start];

Showing the % of the population infectious:

plot2=Show[plot1,
Plot[Evaluate@Table[ Sum[100*i[all[[c,1]],all[[c,2]],t]/.solution[partemp,start],{c,1,jcum}],{jcum,1,Length[all]}],{t,0,maxtime},PlotRange->{Automatic,{0,maxy}},
PlotStyle->{,,,}],
Plot[100((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black,Dashed}],
Plot[100((((TRYδ)/(1-TRYx)) (TRYβV-TRYκ))/(TRYβV (((TRYδ)/(1-TRYx))+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black}],
ImageSize->250
]

Seroconversion affects the dynamics (curves compared to the previous shaded plot) but not the initial speed of spread of the variant (described by the selection coefficient) or the equilibrium.

Clear[solution, i, s, r, ages, stages, n, neqns, nstart, nvars, didt, dsdt, drdt];

Alt Figure S3A: more transmissible & immune evasive  (resident and variant leaky)

maxtime=500;
maxy=20;

Parameters:

TRYp=TRYp; (*No heterogeneity, so p not relevant now.*)
TRYf=99/100;(*Starting fraction of resident*)
TRYolder=13/100; (*~ Age 70+ from StatsCan 5022509/38929902 *)
	(*https://www150.statcan.gc.ca/t1/tbl1/en/cv.action?pid=1710000501*)
TRYκ=1/5; (*Five day infectious period*)
TRYδ=1/125; (*Waning rate based on Menegale et al. (2023)*) 
TRYδO=1/125; (*No clear age difference noted in Menegale et al. (2023)*) 
TRYconvert=10/10; (*Only this fraction of cases seroconvert or boost immunity*)
TRYleak=1/10;
TRYleakV=0.153;

TRYinf=1/50; (*Assumed initial endemic frequency of infected individuals*)

Now the equilibrium must be numerically calculated:

Simplify[Solve[{0==δ R- β1 S I1,
0==β1 S I1+ ξ1 β1 R I1-κ I1,
0==κ I1-δ R -ξ1 β1 R I1,
1==S+I1+R},{S,I1,R}],{β1>0,ξ1>0,δ>0,κ>0}]

startme=%[[3]];

Finding the value of β that gives the target incidence (TRYinf):

Solve[(S+I1+R/.I1->TRYinf/.startme/.δ->TRYδ/.κ->TRYκ/.ξ1->TRYleak)==1,β1]

TRYβ=β1/.%[[2]]

(*TRYβ=((TRYδ TRYκ)/(TRYδ-TRYinf TRYδ-TRYconvert TRYinf TRYκ));*)

((TRYβ)/(TRYκ))//N

Starting equilibrium would be:

start={s0->S,i0->I1,r0->R}/.startme/.δ->TRYδ/.κ->TRYκ/.ξ1->TRYleak/.β1->TRYβ;
N[start]

Variant has a higher transmissibility (set to have the same selection coefficient as for the immune evasive variant):

Flatten[Solve[(TRYleakV-TRYleak)(r0/.start)βV+((βV-TRYβ)/(TRYβ))TRYκ==0.0833,βV]]

TRYβV=βV/.%;

and infects earlier during waning by allowing leaky immunity.

New equilibrium:

neweq=((-δ-κ+β1 ξ1+(Sqrt[(δ^2)+(((κ-β1 ξ1))^2)+2 δ (κ+β1 ξ1-2 κ ξ1)]))/(2 β1 ξ1))/.ξ1->TRYleakV/.β1->TRYβV/.δ->TRYδ/.κ->TRYκ//N

[NOTE: For the ease of colouring, the first stage class is the variant.]

{ages=2,stages=2};

Clear[solution]
solution[{p_,f_,older_,β_,βV_,δY_,δO_,κ_,q_},start]:=solution[{p,f,older,β,βV,δY,δO,κ,q},start]=
Block[{ages=2,stages=2},
dsdt[j_,t_]=(1-q)*Sum[κ[j,jj]*i[j,jj,t],{jj,1,stages}]+δ[j]*r[j,t]-Sum[β[k,kk]*i[k,kk,t],{k,1,ages},{kk,1,stages}]*s[j,t] ;didt[j_,1,t_]=TRYleakV Sum[β[k,1]*i[k,1,t],{k,1,ages}]*r[j,t]+
Sum[β[k,1]*i[k,1,t],{k,1,ages}]*s[j,t] -κ[j,1]*i[j,1,t];didt[j_,2,t_]=TRYleak Sum[β[k,2]*i[k,2,t],{k,1,ages}]*r[j,t]+Sum[β[k,2]*i[k,2,t],{k,1,ages}]*s[j,t] -κ[j,2]*i[j,2,t];drdt[j_,t_]=q*Sum[κ[j,jj]*i[j,jj,t],{jj,1,stages}]-δ[j]*r[j,t]-TRYleakV Sum[β[k,1]*i[k,1,t],{k,1,ages}]*r[j,t]-TRYleak Sum[β[k,2]*i[k,2,t],{k,1,ages}]*r[j,t];pars={β[j_,1]->βV,β[j_,2]->β,δ[1]->δY,δ[2]->δO,κ[j_,jj_]->κ}; 
frac[1,1]=(1-older)(1-f);
frac[2,1]=older(1-f);
frac[1,2]=(1-older)f;
frac[2,2]=older f;
nvars=Drop[Flatten[Table[{s[j,t],Table[i[j,jj,t],{jj,1,stages}],r[j,t]},{j,1,ages}]],-1];
r[ages,t_]=1-Total[nvars];
neqns=Drop[Flatten[Table[{D[s[j,t],t]==dsdt[j,t],Table[D[i[j,jj,t],t]==didt[j,jj,t],{jj,1,stages}],D[r[j,t],t]==drdt[j,t]},{j,1,ages}]],-1];
nstart=Drop[Flatten[Table[{
Table[i[j,jj,0]==frac[j,jj] i0,{jj,1,stages}],
s[j,0]==(frac[j,1]+frac[j,2])*s0,r[j,0]==(frac[j,1]+frac[j,2]) r0},{j,1,ages}]]/.start,-1];
NDSolve[Flatten[{nstart/.pars,neqns/.pars}],nvars,{t,0,maxtime}]
]

all=Flatten[Table[{j,jj},{j,1,ages},{jj,1,stages}],1]

labels={"<70,Res","<70,Mutant","70+,Res","70+,Mutant"};

colours={RGBColor[0.4,0.4,0.8],ColorData["Pastel",1],RGBColor[0.4,0.6,0.4], ColorData["Pastel",3/4]}

coltab=Join[{{1->{0,colours[[1]]}}},Table[{i->{{i-1},colours[[i]]}},{i,2,Length[all]}]];
Table[{labels[[i]],colours[[i]]},{i,1,Length[all]}]//MatrixForm

partemp={TRYp,TRYf,TRYolder,TRYβ,TRYβV,TRYδ,TRYδO,TRYκ,TRYconvert};

solution[partemp,start];

Showing the % of the population infectious:

plot1 =Show[
Plot[Evaluate@Table[ Sum[100*i[all[[c,1]],all[[c,2]],t]/.solution[partemp,start],{c,1,jcum}],{jcum,1,Length[all]}],{t,0,maxtime},PlotRange->{Automatic,{0,maxy}},
PlotStyle->None,Filling->{{1->{0,}},{2->{{1},}},{3->{{2},}},{4->{{3},}}}],
Plot[100((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black,Dashed}],
Plot[100neweq,{t,0,maxtime},PlotStyle->{Black}],
ImageSize->250
]

Mutant fraction change between t=0 and t=50:

((Sum[i[j,1,t],{j,1,ages}])/(Sum[i[j,jj,t],{j,1,ages},{jj,1,stages}]))/.solution[partemp,start]/.t->50

Solve[Exp[s 50] ==( ((p50/(1-p50))/(p0/(1-p0)))/.p50->%/.p0->1-TRYf),s][[1]]

This is consistent with the predicted selection coefficient:

((TRYβV-TRYβ)/(TRYβ))TRYκ+(TRYleakV-TRYleak)(r0/.start)TRYβV//N

Adjusted leakiness of variant to get target selection coefficient in Alt Figure S3C):

Solve[(((TRYβV-TRYβ)/(TRYβ))TRYκ+(leakV-TRYleak)(r0/.start)TRYβV==0.0833),leakV]

{{leakV -> 0.613542}}

% rise in the endemic equilibrium:

100*(((neweq-TRYinf))/TRYinf)

Clear[solution, i, s, r, ages, stages, n, neqns, nstart, nvars, didt, dsdt, drdt];

Alt Figure S3B: more transmissible & immune evasive  (resident and variant leaky)

maxtime=500;
maxy=20;

Parameters:

TRYp=TRYp; (*No heterogeneity, so p not relevant now.*)
TRYf=99/100;(*Starting fraction of resident*)
TRYolder=13/100; (*~ Age 70+ from StatsCan 5022509/38929902 *)
	(*https://www150.statcan.gc.ca/t1/tbl1/en/cv.action?pid=1710000501*)
TRYκ=1/5; (*Five day infectious period*)
TRYδ=1/125; (*Waning rate based on Menegale et al. (2023)*) 
TRYδO=1/125; (*No clear age difference noted in Menegale et al. (2023)*) 
TRYconvert=10/10; (*Only this fraction of cases seroconvert or boost immunity*)
TRYleak=1/10;
TRYleakV=0.307;

TRYinf=1/50; (*Assumed initial endemic frequency of infected individuals*)

Now the equilibrium must be numerically calculated:

Simplify[Solve[{0==δ R- β1 S I1,
0==β1 S I1+ ξ1 β1 R I1-κ I1,
0==κ I1-δ R -ξ1 β1 R I1,
1==S+I1+R},{S,I1,R}],{β1>0,ξ1>0,δ>0,κ>0}]

startme=%[[3]];

Finding the value of β that gives the target incidence (TRYinf):

Solve[(S+I1+R/.I1->TRYinf/.startme/.δ->TRYδ/.κ->TRYκ/.ξ1->TRYleak)==1,β1]

TRYβ=β1/.%[[2]]

(*TRYβ=((TRYδ TRYκ)/(TRYδ-TRYinf TRYδ-TRYconvert TRYinf TRYκ));*)

((TRYβ)/(TRYκ))//N

Starting equilibrium would be:

start={s0->S,i0->I1,r0->R}/.startme/.δ->TRYδ/.κ->TRYκ/.ξ1->TRYleak/.β1->TRYβ;
N[start]

Variant has a higher transmissibility (set to have the same selection coefficient as for the immune evasive variant):

Flatten[Solve[(TRYleakV-TRYleak)(r0/.start)βV+((βV-TRYβ)/(TRYβ))TRYκ==0.0833,βV]]

TRYβV=βV/.%;

and infects earlier during waning by allowing leaky immunity.

New equilibrium:

neweq=((-δ-κ+β1 ξ1+(Sqrt[(δ^2)+(((κ-β1 ξ1))^2)+2 δ (κ+β1 ξ1-2 κ ξ1)]))/(2 β1 ξ1))/.ξ1->TRYleakV/.β1->TRYβV/.δ->TRYδ/.κ->TRYκ//N

[NOTE: For the ease of colouring, the first stage class is the variant.]

{ages=2,stages=2};

Clear[solution]
solution[{p_,f_,older_,β_,βV_,δY_,δO_,κ_,q_},start]:=solution[{p,f,older,β,βV,δY,δO,κ,q},start]=
Block[{ages=2,stages=2},
dsdt[j_,t_]=(1-q)*Sum[κ[j,jj]*i[j,jj,t],{jj,1,stages}]+δ[j]*r[j,t]-Sum[β[k,kk]*i[k,kk,t],{k,1,ages},{kk,1,stages}]*s[j,t] ;didt[j_,1,t_]=TRYleakV Sum[β[k,1]*i[k,1,t],{k,1,ages}]*r[j,t]+
Sum[β[k,1]*i[k,1,t],{k,1,ages}]*s[j,t] -κ[j,1]*i[j,1,t];didt[j_,2,t_]=TRYleak Sum[β[k,2]*i[k,2,t],{k,1,ages}]*r[j,t]+Sum[β[k,2]*i[k,2,t],{k,1,ages}]*s[j,t] -κ[j,2]*i[j,2,t];drdt[j_,t_]=q*Sum[κ[j,jj]*i[j,jj,t],{jj,1,stages}]-δ[j]*r[j,t]-TRYleakV Sum[β[k,1]*i[k,1,t],{k,1,ages}]*r[j,t]-TRYleak Sum[β[k,2]*i[k,2,t],{k,1,ages}]*r[j,t];pars={β[j_,1]->βV,β[j_,2]->β,δ[1]->δY,δ[2]->δO,κ[j_,jj_]->κ}; 
frac[1,1]=(1-older)(1-f);
frac[2,1]=older(1-f);
frac[1,2]=(1-older)f;
frac[2,2]=older f;
nvars=Drop[Flatten[Table[{s[j,t],Table[i[j,jj,t],{jj,1,stages}],r[j,t]},{j,1,ages}]],-1];
r[ages,t_]=1-Total[nvars];
neqns=Drop[Flatten[Table[{D[s[j,t],t]==dsdt[j,t],Table[D[i[j,jj,t],t]==didt[j,jj,t],{jj,1,stages}],D[r[j,t],t]==drdt[j,t]},{j,1,ages}]],-1];
nstart=Drop[Flatten[Table[{
Table[i[j,jj,0]==frac[j,jj] i0,{jj,1,stages}],
s[j,0]==(frac[j,1]+frac[j,2])*s0,r[j,0]==(frac[j,1]+frac[j,2]) r0},{j,1,ages}]]/.start,-1];
NDSolve[Flatten[{nstart/.pars,neqns/.pars}],nvars,{t,0,maxtime}]
]

all=Flatten[Table[{j,jj},{j,1,ages},{jj,1,stages}],1]

labels={"<70,Res","<70,Mutant","70+,Res","70+,Mutant"};

colours={RGBColor[0.4,0.4,0.8],ColorData["Pastel",1],RGBColor[0.4,0.6,0.4], ColorData["Pastel",3/4]}

coltab=Join[{{1->{0,colours[[1]]}}},Table[{i->{{i-1},colours[[i]]}},{i,2,Length[all]}]];
Table[{labels[[i]],colours[[i]]},{i,1,Length[all]}]//MatrixForm

partemp={TRYp,TRYf,TRYolder,TRYβ,TRYβV,TRYδ,TRYδO,TRYκ,TRYconvert};

solution[partemp,start];

Showing the % of the population infectious:

plot1 =Show[
Plot[Evaluate@Table[ Sum[100*i[all[[c,1]],all[[c,2]],t]/.solution[partemp,start],{c,1,jcum}],{jcum,1,Length[all]}],{t,0,maxtime},PlotRange->{Automatic,{0,maxy}},
PlotStyle->None,Filling->{{1->{0,}},{2->{{1},}},{3->{{2},}},{4->{{3},}}}],
Plot[100((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black,Dashed}],
Plot[100neweq,{t,0,maxtime},PlotStyle->{Black}],
ImageSize->250
]

Mutant fraction change between t=0 and t=50:

((Sum[i[j,1,t],{j,1,ages}])/(Sum[i[j,jj,t],{j,1,ages},{jj,1,stages}]))/.solution[partemp,start]/.t->50

Solve[Exp[s 50] ==( ((p50/(1-p50))/(p0/(1-p0)))/.p50->%/.p0->1-TRYf),s][[1]]

This is consistent with the predicted selection coefficient:

((TRYβV-TRYβ)/(TRYβ))TRYκ+(TRYleakV-TRYleak)(r0/.start)TRYβV//N

Adjusted leakiness of variant to get target selection coefficient:

Solve[(((TRYβV-TRYβ)/(TRYβ))TRYκ+(leakV-TRYleak)(r0/.start)TRYβV==0.0833),leakV]

{{leakV -> 0.613542}}

% rise in the endemic equilibrium:

100*(((neweq-TRYinf))/TRYinf)

Clear[solution, i, s, r, ages, stages, n, neqns, nstart, nvars, didt, dsdt, drdt];

Alt Figure S3C: persistently immune evasive (resident and variant leaky)

maxtime=500;
maxy=20;

Parameters:

TRYp=TRYp; (*No heterogeneity, so p not relevant now.*)
TRYf=99/100;(*Starting fraction of resident*)
TRYolder=13/100; (*~ Age 70+ from StatsCan 5022509/38929902 *)
	(*https://www150.statcan.gc.ca/t1/tbl1/en/cv.action?pid=1710000501*)
TRYκ=1/5; (*Five day infectious period*)
TRYδ=1/125; (*Waning rate based on Menegale et al. (2023)*) 
TRYδO=1/125; (*No clear age difference noted in Menegale et al. (2023)*) 
TRYconvert=10/10; (*Only this fraction of cases seroconvert or boost immunity*)
TRYleak=1/10;
TRYleakV=0.614;

TRYinf=1/50; (*Assumed initial endemic frequency of infected individuals*)

Now the equilibrium must be numerically calculated:

Simplify[Solve[{0==δ R- β1 S I1,
0==β1 S I1+ ξ1 β1 R I1-κ I1,
0==κ I1-δ R -ξ1 β1 R I1,
1==S+I1+R},{S,I1,R}],{β1>0,ξ1>0,δ>0,κ>0}]

startme=%[[3]];

Finding the value of β that gives the target incidence (TRYinf):

Solve[(S+I1+R/.I1->TRYinf/.startme/.δ->TRYδ/.κ->TRYκ/.ξ1->TRYleak)==1,β1]

TRYβ=β1/.%[[2]]

(*TRYβ=((TRYδ TRYκ)/(TRYδ-TRYinf TRYδ-TRYconvert TRYinf TRYκ));*)

((TRYβ)/(TRYκ))//N

Variant has the same transmissibility:

TRYβV=TRYβ;

but infects earlier during waning by allowing leaky immunity.

New equilibrium:

neweq=((-δ-κ+β1 ξ1+(Sqrt[(δ^2)+(((κ-β1 ξ1))^2)+2 δ (κ+β1 ξ1-2 κ ξ1)]))/(2 β1 ξ1))/.ξ1->TRYleakV/.β1->TRYβV/.δ->TRYδ/.κ->TRYκ//N

Starting equilibrium would be:

start={s0->S,i0->I1,r0->R}/.startme/.δ->TRYδ/.κ->TRYκ/.ξ1->TRYleak/.β1->TRYβ;
N[start]

[NOTE: For the ease of colouring, the first stage class is the variant.]

{ages=2,stages=2};

Clear[solution]
solution[{p_,f_,older_,β_,βV_,δY_,δO_,κ_,q_},start]:=solution[{p,f,older,β,βV,δY,δO,κ,q},start]=
Block[{ages=2,stages=2},
dsdt[j_,t_]=(1-q)*Sum[κ[j,jj]*i[j,jj,t],{jj,1,stages}]+δ[j]*r[j,t]-Sum[β[k,kk]*i[k,kk,t],{k,1,ages},{kk,1,stages}]*s[j,t] ;didt[j_,1,t_]=TRYleakV Sum[β[k,1]*i[k,1,t],{k,1,ages}]*r[j,t]+
Sum[β[k,1]*i[k,1,t],{k,1,ages}]*s[j,t] -κ[j,1]*i[j,1,t];didt[j_,2,t_]=TRYleak Sum[β[k,2]*i[k,2,t],{k,1,ages}]*r[j,t]+Sum[β[k,2]*i[k,2,t],{k,1,ages}]*s[j,t] -κ[j,2]*i[j,2,t];drdt[j_,t_]=q*Sum[κ[j,jj]*i[j,jj,t],{jj,1,stages}]-δ[j]*r[j,t]-TRYleakV Sum[β[k,1]*i[k,1,t],{k,1,ages}]*r[j,t]-TRYleak Sum[β[k,2]*i[k,2,t],{k,1,ages}]*r[j,t];pars={β[j_,1]->βV,β[j_,2]->β,δ[1]->δY,δ[2]->δO,κ[j_,jj_]->κ}; 
frac[1,1]=(1-older)(1-f);
frac[2,1]=older(1-f);
frac[1,2]=(1-older)f;
frac[2,2]=older f;
nvars=Drop[Flatten[Table[{s[j,t],Table[i[j,jj,t],{jj,1,stages}],r[j,t]},{j,1,ages}]],-1];
r[ages,t_]=1-Total[nvars];
neqns=Drop[Flatten[Table[{D[s[j,t],t]==dsdt[j,t],Table[D[i[j,jj,t],t]==didt[j,jj,t],{jj,1,stages}],D[r[j,t],t]==drdt[j,t]},{j,1,ages}]],-1];
nstart=Drop[Flatten[Table[{
Table[i[j,jj,0]==frac[j,jj] i0,{jj,1,stages}],
s[j,0]==(frac[j,1]+frac[j,2])*s0,r[j,0]==(frac[j,1]+frac[j,2]) r0},{j,1,ages}]]/.start,-1];
NDSolve[Flatten[{nstart/.pars,neqns/.pars}],nvars,{t,0,maxtime}]
]

all=Flatten[Table[{j,jj},{j,1,ages},{jj,1,stages}],1]

labels={"<70,Res","<70,Mutant","70+,Res","70+,Mutant"};

colours={RGBColor[0.4,0.4,0.8],ColorData["Pastel",1],RGBColor[0.4,0.6,0.4], ColorData["Pastel",3/4]}

coltab=Join[{{1->{0,colours[[1]]}}},Table[{i->{{i-1},colours[[i]]}},{i,2,Length[all]}]];
Table[{labels[[i]],colours[[i]]},{i,1,Length[all]}]//MatrixForm

partemp={TRYp,TRYf,TRYolder,TRYβ,TRYβV,TRYδ,TRYδO,TRYκ,TRYconvert};

solution[partemp,start];

Showing the % of the population infectious:

plot1 =Show[
Plot[Evaluate@Table[ Sum[100*i[all[[c,1]],all[[c,2]],t]/.solution[partemp,start],{c,1,jcum}],{jcum,1,Length[all]}],{t,0,maxtime},PlotRange->{Automatic,{0,maxy}},
PlotStyle->None,Filling->{{1->{0,}},{2->{{1},}},{3->{{2},}},{4->{{3},}}}],
Plot[100((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black,Dashed}],
Plot[100neweq,{t,0,maxtime},PlotStyle->{Black}],
ImageSize->250
]

Mutant fraction change between t=0 and t=50:

((Sum[i[j,1,t],{j,1,ages}])/(Sum[i[j,jj,t],{j,1,ages},{jj,1,stages}]))/.solution[partemp,start]/.t->50

Solve[Exp[s 50] ==( ((p50/(1-p50))/(p0/(1-p0)))/.p50->%/.p0->1-TRYf),s][[1]]

This is consistent with the predicted selection coefficient:

((TRYβV-TRYβ)/(TRYβ))TRYκ+(TRYleakV-TRYleak)(r0/.start)TRYβV//N

Adjusted leakiness of variant to get target selection coefficient:

Solve[(((TRYβV-TRYβ)/(TRYβ))TRYκ+(leakV-TRYleak)(r0/.start)TRYβV==0.0833),leakV]

{{leakV -> 0.613542}}

% rise in the endemic equilibrium:

100*(((neweq-TRYinf))/TRYinf)

Clear[solution, i, s, r, ages, stages, n, neqns, nstart, nvars, didt, dsdt, drdt];

Alternative: (Subscript[SIR, 2]) model
-> Less oscillatory

Model includes two waning classes, corresponding to high immunity early, decaying to a lower level of immunity that can be infected by a new variant.

We let the waning rate in the two classes be:

{(δ/(1-x)),(δ/x)};

so that the expected time to waning for the resident is:

1/%[[1]]+1/%[[2]]//Factor

x determines the fraction of the waning period in which more immune evasive variants can infect.

Figure S2D: more transmissible

maxtime=500;
maxy=8;

Parameters:

TRYp=TRYp; (*No heterogeneity, so p not relevant now.*)
TRYf=99/100;(*Starting fraction of resident*)
TRYolder=13/100; (*~ Age 70+ from StatsCan 5022509/38929902 *)
	(*https://www150.statcan.gc.ca/t1/tbl1/en/cv.action?pid=1710000501*)
TRYκ=1/5; (*Five day infectious period*)
TRYδ=1/125; (*Waning rate based on Menegale et al. (2023)*) 
TRYδO=1/125; (*No clear age difference noted in Menegale et al. (2023)*) 
TRYconvert=10/10; (*Only this fraction of cases seroconvert or boost immunity*)
TRYx=2/5;

i->((δ (β-κ))/(β (δ+q κ)))

Solve[((δ (β-κ))/(β (δ+q κ)))==inf,β]

TRYinf=1/50; (*Assumed initial endemic frequency of infected individuals*)
TRYβ=((TRYδ TRYκ)/(TRYδ-TRYinf TRYδ-TRYconvert TRYinf TRYκ));

((TRYβ)/(TRYκ))//N

Variant has a higher transmissibility (set to have the same selection coefficient as for the immune evasive variant):

Flatten[Solve[((βV-TRYβ)/(TRYβ))TRYκ==0.0833,βV]]

TRYβV=βV/.%;

% increase in β:

100*((TRYβV-TRYβ)/(TRYβ))

and now has no immune evasive properties (neither the resident nor the variant can infect individuals in the RL class).

Starting equilibrium would be (r0 includes all of the resistance classes):
{s->(κ/β),i->((δ (β-κ))/(β (δ+q κ)))}

start={s0->((TRYκ)/(TRYβ)),i0->((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))),r0->1-((TRYκ)/(TRYβ))-((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert)))}

[NOTE: For the ease of colouring, the first stage class is the variant.]

{ages=2,stages=2};

Clear[solution]
solution[{p_,f_,older_,β_,βV_,δY_,δO_,κ_,q_},start]:=solution[{p,f,older,β,βV,δY,δO,κ,q},start]=
Block[{ages=2,stages=2},
dsdt[j_,t_]=(1-q)*Sum[κ[j,jj]*i[j,jj,t],{jj,1,stages}]+((δ[j])/TRYx)*rL[j,t]-Sum[β[k,kk]*i[k,kk,t],{k,1,ages},{kk,1,stages}]*s[j,t] ;didt[j_,1,t_]=Sum[β[k,1]*i[k,1,t],{k,1,ages}]*s[j,t] -κ[j,1]*i[j,1,t];didt[j_,2,t_]=Sum[β[k,2]*i[k,2,t],{k,1,ages}]*s[j,t] -κ[j,2]*i[j,2,t];drHdt[j_,t_]=q*Sum[κ[j,jj]*i[j,jj,t],{jj,1,stages}]-((δ[j])/(1-TRYx))*rH[j,t];drLdt[j_,t_]=((δ[j])/(1-TRYx))*rH[j,t]-((δ[j])/TRYx)*rL[j,t];pars={β[j_,1]->βV,β[j_,2]->β,δ[1]->δY,δ[2]->δO,κ[j_,jj_]->κ}; 
frac[1,1]=(1-older)(1-f);
frac[2,1]=older(1-f);
frac[1,2]=(1-older)f;
frac[2,2]=older f;
nvars=Drop[Flatten[Table[{s[j,t],Table[i[j,jj,t],{jj,1,stages}],rH[j,t],rL[j,t]},{j,1,ages}]],-1];
rL[ages,t_]=1-Total[nvars];
neqns=Drop[Flatten[Table[{D[s[j,t],t]==dsdt[j,t],Table[D[i[j,jj,t],t]==didt[j,jj,t],{jj,1,stages}],D[rH[j,t],t]==drHdt[j,t],D[rL[j,t],t]==drLdt[j,t]},{j,1,ages}]],-1];
nstart=Drop[Flatten[Table[{
Table[i[j,jj,0]==frac[j,jj] i0,{jj,1,stages}],
s[j,0]==(frac[j,1]+frac[j,2])*s0,rH[j,0]==(frac[j,1]+frac[j,2]) (1-TRYx)r0,rL[j,0]==(frac[j,1]+frac[j,2]) TRYx r0},{j,1,ages}]]/.start,-1];
NDSolve[Flatten[{nstart/.pars,neqns/.pars}],nvars,{t,0,maxtime}]
]

all=Flatten[Table[{j,jj},{j,1,ages},{jj,1,stages}],1]

labels={"<70,Res","<70,Mutant","70+,Res","70+,Mutant"};

colours={RGBColor[0.4,0.4,0.8],ColorData["Pastel",1],RGBColor[0.4,0.6,0.4], ColorData["Pastel",3/4]}

coltab=Join[{{1->{0,colours[[1]]}}},Table[{i->{{i-1},colours[[i]]}},{i,2,Length[all]}]];
Table[{labels[[i]],colours[[i]]},{i,1,Length[all]}]//MatrixForm

partemp={TRYp,TRYf,TRYolder,TRYβ,TRYβV,TRYδ,TRYδO,TRYκ,TRYconvert};

solution[partemp,start];

Showing the % of the population infectious:

plot1 =Show[
Plot[Evaluate@Table[ Sum[100*i[all[[c,1]],all[[c,2]],t]/.solution[partemp,start],{c,1,jcum}],{jcum,1,Length[all]}],{t,0,maxtime},PlotRange->{Automatic,{0,maxy}},
PlotStyle->None,Filling->{{1->{0,}},{2->{{1},}},{3->{{2},}},{4->{{3},}}}],
Plot[100((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black,Dashed}],
Plot[100((TRYδ (TRYβV-TRYκ))/(TRYβV (TRYδ+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black}],
ImageSize->250
]

Mutant fraction change between t=0 and t=50:

((Sum[i[j,1,t],{j,1,ages}])/(Sum[i[j,jj,t],{j,1,ages},{jj,1,stages}]))/.solution[partemp,start]/.t->50

Solve[Exp[s 50] ==( ((p50/(1-p50))/(p0/(1-p0)))/.p50->%/.p0->1-TRYf),s][[1]]

This is consistent with the predicted selection coefficient:

((TRYβV-TRYβ)/(TRYβ))TRYκ+TRYx(((TRYκ TRYconvert(TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))))TRYβV//N

New equilibrium:

((((TRYδ+Δδ)(TRYβV-TRYκ))/(TRYβV ((TRYδ+Δδ)+TRYκ TRYconvert))))/.Δδ->(TRYx/(1-TRYx))TRYδ//N

% rise in the endemic equilibrium:

100*(((%-TRYinf))/TRYinf)

The following illustrates the relatively minor effect caused by reducing seroconversion (to q = 80%), as long as δ/q is held constant:

partemp={TRYp,TRYf,TRYolder,TRYβ,TRYβV,(8/10)TRYδ,(8/10)TRYδO,TRYκ,(8/10)};

solution[partemp,start];

Showing the % of the population infectious:

plot2=Show[plot1,
Plot[Evaluate@Table[ Sum[100*i[all[[c,1]],all[[c,2]],t]/.solution[partemp,start],{c,1,jcum}],{jcum,1,Length[all]}],{t,0,maxtime},PlotRange->{Automatic,{0,maxy}},
PlotStyle->{,,,}],
Plot[100((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black,Dashed}],
Plot[100((TRYδ (TRYβV-TRYκ))/(TRYβV (TRYδ+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black}],
ImageSize->250
]

Seroconversion affects the dynamics (curves compared to the previous shaded plot) but not the initial speed of spread of the variant (described by the selection coefficient) or the equilibrium.

Clear[solution, i, s, r, ages, stages, n, neqns, nstart, nvars, didt, dsdt, drdt];

Figure S2E: persistently immune evasive (increases mean waning rate by 67%)

maxtime=500;
maxy=8;

Parameters:

TRYp=TRYp; (*No heterogeneity, so p not relevant now.*)
TRYf=99/100;(*Starting fraction of resident*)
TRYolder=13/100; (*~ Age 70+ from StatsCan 5022509/38929902 *)
	(*https://www150.statcan.gc.ca/t1/tbl1/en/cv.action?pid=1710000501*)
TRYκ=1/5; (*Five day infectious period*)
TRYδ=1/125; (*Waning rate based on Menegale et al. (2023)*) 
TRYδO=1/125; (*No clear age difference noted in Menegale et al. (2023)*) 
TRYconvert=10/10; (*Only this fraction of cases seroconvert or boost immunity*)
TRYx=2/5;

i->((δ (β-κ))/(β (δ+q κ)))

Solve[((δ (β-κ))/(β (δ+q κ)))==inf,β]

TRYinf=1/50; (*Assumed initial endemic frequency of infected individuals*)
TRYβ=((TRYδ TRYκ)/(TRYδ-TRYinf TRYδ-TRYconvert TRYinf TRYκ));

((TRYβ)/(TRYκ))//N

Variant has the same transmissibility:

TRYβV=TRYβ;

but infects earlier during waning.

Specifically, the mean waning rate (once fixed) increases by 67%:

TRYδ/(1-x)-TRYδ//Factor

((TRYδ+Δδ)/(TRYδ))/.Δδ->(TRYx/(1-TRYx))TRYδ//N

while the mean waiting time decreases by 40%:

1/%

Starting equilibrium would be (r0 includes all of the resistance classes):
{s->(κ/β),i->((δ (β-κ))/(β (δ+q κ)))}

start={s0->((TRYκ)/(TRYβ)),i0->((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))),r0->1-((TRYκ)/(TRYβ))-((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert)))}

[NOTE: For the ease of colouring, the first stage class is the variant.]

{ages=2,stages=2};

Clear[solution]
solution[{p_,f_,older_,β_,βV_,δY_,δO_,κ_,q_},start]:=solution[{p,f,older,β,βV,δY,δO,κ,q},start]=
Block[{ages=2,stages=2},
dsdt[j_,t_]=(1-q)*Sum[κ[j,jj]*i[j,jj,t],{jj,1,stages}]+((δ[j])/TRYx)*rL[j,t]-Sum[β[k,kk]*i[k,kk,t],{k,1,ages},{kk,1,stages}]*s[j,t] ;didt[j_,1,t_]=Sum[β[k,1]*i[k,1,t],{k,1,ages}]*rL[j,t]+
Sum[β[k,1]*i[k,1,t],{k,1,ages}]*s[j,t] -κ[j,1]*i[j,1,t];didt[j_,2,t_]=Sum[β[k,2]*i[k,2,t],{k,1,ages}]*s[j,t] -κ[j,2]*i[j,2,t];drHdt[j_,t_]=q*Sum[κ[j,jj]*i[j,jj,t],{jj,1,stages}]-((δ[j])/(1-TRYx))*rH[j,t];drLdt[j_,t_]=((δ[j])/(1-TRYx))*rH[j,t]-((δ[j])/TRYx)*rL[j,t]-Sum[β[k,1]*i[k,1,t],{k,1,ages}]*rL[j,t];pars={β[j_,1]->βV,β[j_,2]->β,δ[1]->δY,δ[2]->δO,κ[j_,jj_]->κ}; 
frac[1,1]=(1-older)(1-f);
frac[2,1]=older(1-f);
frac[1,2]=(1-older)f;
frac[2,2]=older f;
nvars=Drop[Flatten[Table[{s[j,t],Table[i[j,jj,t],{jj,1,stages}],rH[j,t],rL[j,t]},{j,1,ages}]],-1];
rL[ages,t_]=1-Total[nvars];
neqns=Drop[Flatten[Table[{D[s[j,t],t]==dsdt[j,t],Table[D[i[j,jj,t],t]==didt[j,jj,t],{jj,1,stages}],D[rH[j,t],t]==drHdt[j,t],D[rL[j,t],t]==drLdt[j,t]},{j,1,ages}]],-1];
nstart=Drop[Flatten[Table[{
Table[i[j,jj,0]==frac[j,jj] i0,{jj,1,stages}],
s[j,0]==(frac[j,1]+frac[j,2])*s0,rH[j,0]==(frac[j,1]+frac[j,2]) (1-TRYx)r0,rL[j,0]==(frac[j,1]+frac[j,2]) TRYx r0},{j,1,ages}]]/.start,-1];
NDSolve[Flatten[{nstart/.pars,neqns/.pars}],nvars,{t,0,maxtime}]
]

all=Flatten[Table[{j,jj},{j,1,ages},{jj,1,stages}],1]

labels={"<70,Res","<70,Mutant","70+,Res","70+,Mutant"};

colours={RGBColor[0.4,0.4,0.8],ColorData["Pastel",1],RGBColor[0.4,0.6,0.4], ColorData["Pastel",3/4]}

coltab=Join[{{1->{0,colours[[1]]}}},Table[{i->{{i-1},colours[[i]]}},{i,2,Length[all]}]];
Table[{labels[[i]],colours[[i]]},{i,1,Length[all]}]//MatrixForm

partemp={TRYp,TRYf,TRYolder,TRYβ,TRYβV,TRYδ,TRYδO,TRYκ,TRYconvert};

solution[partemp,start];

Showing the % of the population infectious:

plot1 =Show[
Plot[Evaluate@Table[ Sum[100*i[all[[c,1]],all[[c,2]],t]/.solution[partemp,start],{c,1,jcum}],{jcum,1,Length[all]}],{t,0,maxtime},PlotRange->{Automatic,{0,maxy}},
PlotStyle->None,Filling->{{1->{0,}},{2->{{1},}},{3->{{2},}},{4->{{3},}}}],
Plot[100((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black,Dashed}],
Plot[100((((TRYδ)/(1-TRYx)) (TRYβV-TRYκ))/(TRYβV (((TRYδ)/(1-TRYx))+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black}],
ImageSize->250
]

Mutant fraction change between t=0 and t=50:

((Sum[i[j,1,t],{j,1,ages}])/(Sum[i[j,jj,t],{j,1,ages},{jj,1,stages}]))/.solution[partemp,start]/.t->50

Solve[Exp[s 50] ==( ((p50/(1-p50))/(p0/(1-p0)))/.p50->%/.p0->1-TRYf),s][[1]]

This is consistent with the predicted selection coefficient:

((TRYβV-TRYβ)/(TRYβ))TRYκ+TRYx(((TRYκ TRYconvert(TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))))TRYβV//N

New equilibrium:

((((TRYδ+Δδ)(TRYβV-TRYκ))/(TRYβV ((TRYδ+Δδ)+TRYκ TRYconvert))))/.Δδ->(TRYx/(1-TRYx))TRYδ//N

% rise in the endemic equilibrium:

100*(((%-TRYinf))/TRYinf)

The following illustrates the relatively minor effect caused by reducing seroconversion (to q = 80%), as long as δ/q is held constant:

partemp={TRYp,TRYf,TRYolder,TRYβ,TRYβV,(8/10)TRYδ,(8/10)TRYδO,TRYκ,(8/10)};

solution[partemp,start];

Showing the % of the population infectious:

plot2=Show[plot1,
Plot[Evaluate@Table[ Sum[100*i[all[[c,1]],all[[c,2]],t]/.solution[partemp,start],{c,1,jcum}],{jcum,1,Length[all]}],{t,0,maxtime},PlotRange->{Automatic,{0,maxy}},
PlotStyle->{,,,}],
Plot[100((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black,Dashed}],
Plot[100((((TRYδ)/(1-TRYx)) (TRYβV-TRYκ))/(TRYβV (((TRYδ)/(1-TRYx))+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black}],
ImageSize->250
]

Seroconversion affects the dynamics (curves compared to the previous shaded plot) but not the initial speed of spread of the variant (described by the selection coefficient) or the equilibrium.

Clear[solution, i, s, r, ages, stages, n, neqns, nstart, nvars, didt, dsdt, drdt];

Figure S2F: more transmissible and persistently immune evasive (increases mean waning rate by 25%)

maxtime=500;
maxy=8;

Parameters:

TRYp=TRYp; (*No heterogeneity, so p not relevant now.*)
TRYf=99/100;(*Starting fraction of resident*)
TRYolder=13/100; (*~ Age 70+ from StatsCan 5022509/38929902 *)
	(*https://www150.statcan.gc.ca/t1/tbl1/en/cv.action?pid=1710000501*)
TRYκ=1/5; (*Five day infectious period*)
TRYδ=1/125; (*Waning rate based on Menegale et al. (2023)*) 
TRYδO=1/125; (*No clear age difference noted in Menegale et al. (2023)*) 
TRYconvert=10/10; (*Only this fraction of cases seroconvert or boost immunity*)
TRYx=1/5;

i->((δ (β-κ))/(β (δ+q κ)))

Solve[((δ (β-κ))/(β (δ+q κ)))==inf,β]

TRYinf=1/50; (*Assumed initial endemic frequency of infected individuals*)
TRYβ=((TRYδ TRYκ)/(TRYδ-TRYinf TRYδ-TRYconvert TRYinf TRYκ));

((TRYβ)/(TRYκ))//N

Variant has a higher transmissibility (set to have the same selection coefficient as for the immune evasive variant):

Flatten[Solve[((βV-TRYβ)/(TRYβ))TRYκ+TRYx(((TRYκ TRYconvert(TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))))βV==0.0833,βV]]

TRYβV=βV/.%;

% increase in β:

100*((TRYβV-TRYβ)/(TRYβ))

and infects earlier during waning.

Specifically, the mean waning rate (once fixed) increases by 25%:

TRYδ/(1-x)-TRYδ//Factor

((TRYδ+Δδ)/(TRYδ))/.Δδ->(TRYx/(1-TRYx))TRYδ//N

while the mean waiting time decreases by 20%:

1/%

Starting equilibrium would be (r0 includes all of the resistance classes):
{s->(κ/β),i->((δ (β-κ))/(β (δ+q κ)))}

start={s0->((TRYκ)/(TRYβ)),i0->((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))),r0->1-((TRYκ)/(TRYβ))-((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert)))}

[NOTE: For the ease of colouring, the first stage class is the variant.]

{ages=2,stages=2};

Clear[solution]
solution[{p_,f_,older_,β_,βV_,δY_,δO_,κ_,q_},start]:=solution[{p,f,older,β,βV,δY,δO,κ,q},start]=
Block[{ages=2,stages=2},
dsdt[j_,t_]=(1-q)*Sum[κ[j,jj]*i[j,jj,t],{jj,1,stages}]+((δ[j])/TRYx)*rL[j,t]-Sum[β[k,kk]*i[k,kk,t],{k,1,ages},{kk,1,stages}]*s[j,t] ;didt[j_,1,t_]=Sum[β[k,1]*i[k,1,t],{k,1,ages}]*rL[j,t]+
Sum[β[k,1]*i[k,1,t],{k,1,ages}]*s[j,t] -κ[j,1]*i[j,1,t];didt[j_,2,t_]=Sum[β[k,2]*i[k,2,t],{k,1,ages}]*s[j,t] -κ[j,2]*i[j,2,t];drHdt[j_,t_]=q*Sum[κ[j,jj]*i[j,jj,t],{jj,1,stages}]-((δ[j])/(1-TRYx))*rH[j,t];drLdt[j_,t_]=((δ[j])/(1-TRYx))*rH[j,t]-((δ[j])/TRYx)*rL[j,t]-Sum[β[k,1]*i[k,1,t],{k,1,ages}]*rL[j,t];pars={β[j_,1]->βV,β[j_,2]->β,δ[1]->δY,δ[2]->δO,κ[j_,jj_]->κ}; 
frac[1,1]=(1-older)(1-f);
frac[2,1]=older(1-f);
frac[1,2]=(1-older)f;
frac[2,2]=older f;
nvars=Drop[Flatten[Table[{s[j,t],Table[i[j,jj,t],{jj,1,stages}],rH[j,t],rL[j,t]},{j,1,ages}]],-1];
rL[ages,t_]=1-Total[nvars];
neqns=Drop[Flatten[Table[{D[s[j,t],t]==dsdt[j,t],Table[D[i[j,jj,t],t]==didt[j,jj,t],{jj,1,stages}],D[rH[j,t],t]==drHdt[j,t],D[rL[j,t],t]==drLdt[j,t]},{j,1,ages}]],-1];
nstart=Drop[Flatten[Table[{
Table[i[j,jj,0]==frac[j,jj] i0,{jj,1,stages}],
s[j,0]==(frac[j,1]+frac[j,2])*s0,rH[j,0]==(frac[j,1]+frac[j,2]) (1-TRYx)r0,rL[j,0]==(frac[j,1]+frac[j,2]) TRYx r0},{j,1,ages}]]/.start,-1];
NDSolve[Flatten[{nstart/.pars,neqns/.pars}],nvars,{t,0,maxtime}]
]

all=Flatten[Table[{j,jj},{j,1,ages},{jj,1,stages}],1]

labels={"<70,Res","<70,Mutant","70+,Res","70+,Mutant"};

colours={RGBColor[0.4,0.4,0.8],ColorData["Pastel",1],RGBColor[0.4,0.6,0.4], ColorData["Pastel",3/4]}

coltab=Join[{{1->{0,colours[[1]]}}},Table[{i->{{i-1},colours[[i]]}},{i,2,Length[all]}]];
Table[{labels[[i]],colours[[i]]},{i,1,Length[all]}]//MatrixForm

partemp={TRYp,TRYf,TRYolder,TRYβ,TRYβV,TRYδ,TRYδO,TRYκ,TRYconvert};

solution[partemp,start];

Showing the % of the population infectious:

plot1 =Show[
Plot[Evaluate@Table[ Sum[100*i[all[[c,1]],all[[c,2]],t]/.solution[partemp,start],{c,1,jcum}],{jcum,1,Length[all]}],{t,0,maxtime},PlotRange->{Automatic,{0,maxy}},
PlotStyle->None,Filling->{{1->{0,}},{2->{{1},}},{3->{{2},}},{4->{{3},}}}],
Plot[100((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black,Dashed}],
Plot[100((((TRYδ)/(1-TRYx)) (TRYβV-TRYκ))/(TRYβV (((TRYδ)/(1-TRYx))+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black}],
ImageSize->250
]

Mutant fraction change between t=0 and t=50:

((Sum[i[j,1,t],{j,1,ages}])/(Sum[i[j,jj,t],{j,1,ages},{jj,1,stages}]))/.solution[partemp,start]/.t->50

Solve[Exp[s 50] ==( ((p50/(1-p50))/(p0/(1-p0)))/.p50->%/.p0->1-TRYf),s][[1]]

This is consistent with the predicted selection coefficient:

((TRYβV-TRYβ)/(TRYβ))TRYκ+TRYx(((TRYκ TRYconvert(TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))))TRYβV//N

New equilibrium:

((((TRYδ+Δδ)(TRYβV-TRYκ))/(TRYβV ((TRYδ+Δδ)+TRYκ TRYconvert))))/.Δδ->(TRYx/(1-TRYx))TRYδ//N

% rise in the endemic equilibrium:

100*(((%-TRYinf))/TRYinf)

The following illustrates the relatively minor effect caused by reducing seroconversion (to q = 80%), as long as δ/q is held constant:

partemp={TRYp,TRYf,TRYolder,TRYβ,TRYβV,(8/10)TRYδ,(8/10)TRYδO,TRYκ,(8/10)};

solution[partemp,start];

Showing the % of the population infectious:

plot2=Show[plot1,
Plot[Evaluate@Table[ Sum[100*i[all[[c,1]],all[[c,2]],t]/.solution[partemp,start],{c,1,jcum}],{jcum,1,Length[all]}],{t,0,maxtime},PlotRange->{Automatic,{0,maxy}},
PlotStyle->{,,,}],
Plot[100((TRYδ (TRYβ-TRYκ))/(TRYβ (TRYδ+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black,Dashed}],
Plot[100((((TRYδ)/(1-TRYx)) (TRYβV-TRYκ))/(TRYβV (((TRYδ)/(1-TRYx))+TRYκ TRYconvert))),{t,0,maxtime},PlotStyle->{Black}],
ImageSize->250
]

Seroconversion affects the dynamics (curves compared to the previous shaded plot) but not the initial speed of spread of the variant (described by the selection coefficient) or the equilibrium.

Clear[solution, i, s, r, ages, stages, n, neqns, nstart, nvars, didt, dsdt, drdt];