Supplementary Material: 
“The potential for viral evolution during the COVID-19 pandemic”

Troy Day, Sylvain Gandon, Sébastien Lion, & Sarah P. Otto
May 16, 2020

Analysis of model with genetic variation (wildtype and mutant subtypes)

Based on Figure 2:

We develop ordinary differential equations for both the wildtype (w) and the mutant (m):

dSdt=-((βI*Iw+βA*Aw+βP*Pw)+(βIm*Im+βAm*Am+βP*Pm))*S;
dEwdt=(βI*Iw+βA*Aw+βP*Pw)*S-κE*Ew;
dPwdt=(1-f)*κE*Ew-κP*Pw;
dIwdt=κP*Pw-α Iw - κI Iw;
dAwdt=f*κE*Ew-κA*Aw;

dEmdt=(βIm*Im+βAm*Am+βPm*Pm)*S-κEm*Em;
dPmdt=(1-fm)*κEm*Em-κPm*Pm;
dImdt=κPm*Pm-αm Im - κIm Im;
dAmdt=fm*κEm*Em-κAm*Am;

We assume that the disease is fairly rare and study its growth by analysing the properties of the local stability matrix:

start={Ew->0,Aw->0,Pw->0,Iw->0,Em->0,Am->0,Pm->0,Im->0};

stabmat={{D[dSdt,S],D[dSdt,Ew],D[dSdt,Pw],D[dSdt,Iw],D[dSdt,Aw],D[dSdt,Em],D[dSdt,Pm],D[dSdt,Im],D[dSdt,Am]},
{D[dEwdt,S],D[dEwdt,Ew],D[dEwdt,Pw],D[dEwdt,Iw],D[dEwdt,Aw],D[dEwdt,Em],D[dEwdt,Pm],D[dEwdt,Im],D[dEwdt,Am]},
{D[dPwdt,S],D[dPwdt,Ew],D[dPwdt,Pw],D[dPwdt,Iw],D[dPwdt,Aw],D[dPwdt,Em],D[dPwdt,Pm],D[dPwdt,Im],D[dPwdt,Am]},
{D[dIwdt,S],D[dIwdt,Ew],D[dIwdt,Pw],D[dIwdt,Iw],D[dIwdt,Aw],D[dIwdt,Em],D[dIwdt,Pm],D[dIwdt,Im],D[dIwdt,Am]},{D[dAwdt,S],D[dAwdt,Ew],D[dAwdt,Pw],D[dAwdt,Iw],D[dAwdt,Aw],D[dAwdt,Em],D[dAwdt,Pm],D[dAwdt,Im],D[dAwdt,Am]},
{D[dEmdt,S],D[dEmdt,Ew],D[dEmdt,Pw],D[dEmdt,Iw],D[dEmdt,Aw],D[dEmdt,Em],D[dEmdt,Pm],D[dEmdt,Im],D[dEmdt,Am]},
{D[dPmdt,S],D[dPmdt,Ew],D[dPmdt,Pw],D[dPmdt,Iw],D[dPmdt,Aw],D[dPmdt,Em],D[dPmdt,Pm],D[dPmdt,Im],D[dPmdt,Am]},
{D[dImdt,S],D[dImdt,Ew],D[dImdt,Pw],D[dImdt,Iw],D[dImdt,Aw],D[dImdt,Em],D[dImdt,Pm],D[dImdt,Im],D[dImdt,Am]},
{D[dAmdt,S],D[dAmdt,Ew],D[dAmdt,Pw],D[dAmdt,Iw],D[dAmdt,Aw],D[dAmdt,Em],D[dAmdt,Pm],D[dAmdt,Im],D[dAmdt,Am]}}/.start;

In the special case where the exposed class is highly transitory and most infections occur during the pre-symptomatic phase, we have the simpler equations:

	dPwdt=(1-f)*βP*Pw-κP*Pw
	dPmdt = (1-fm)*βPm*Pm-κPm*Pm;
	
The change in allele frequency, p=(Pm/(Pm+Pw)), is thus given by (dp/dt)=((d((Pm/(Pm+Pw))))/dt), which is:

D[((Pm[t])/(Pw[t]+Pm[t])),t]/.Pm'[t]->(1-fm)*βPm*S Pm-κPm*Pm/.Pw'[t]->(1-f)*βP*S*Pw-κP*Pw/.Pm[t]->Pm/.Pw[t]->Pw/.fm->f+Δf/.βPm->βP+Δβ/.κPm->κP+Δκ/.{Pm->p(Pm+Pw),Pw->(1-p)(Pm+Pw)}//Factor

Assuming that Δf Δβ is small, relative to the other terms (weak selection), gives equations (1) and (2).

Returning to the generic model:

stabmat//MatrixForm

The eigenvalues, λ, are determined by the roots of the following characteristic polynomial:

charpoly=Det[λ IdentityMatrix[9]-stabmat]//Factor

This separates into a characteristic polynomial describing the growth of the resident:

poly=Collect[ charpoly[[2]],λ,Factor]

poly/.βI->0/.βA->0//Factor

And a characteristic polynomial describing the growth of the mutant:

polym=Collect[ charpoly[[3]],λ,Factor]/.λ->λm

As expected, when most individuals are susceptible, mutations spread independently and the one with the largest eigenvalue will spread fastest.  To determine the difference in the leading eigenvalue, we assume that the mutational effects are small (all of the “m” terms above are near the wildtype rate, with the difference of order ϵ). We can then calculate Δλm=λm-λ, by taking the Taylor series and solving the characteristic polynomial for Δλm.

smallmutant={αm->α+ϵ*Δα,fm->f+ϵ*Δf,κIm->κI+ϵ*ΔκI,κEm->κE+ϵ*ΔκE,κAm->κA+ϵ*ΔκA,κPm->κP+ϵ*ΔκP,βAm->βA+ϵ*ΔβA,βPm->βP+ϵ*ΔβP,βIm->βI+ϵ*ΔβI};
smallterms={Δα,Δf,ΔκI,ΔκE,ΔκA,ΔκP,ΔβA,ΔβP,ΔβI};

As expected, polym is zero to O(1) (reduces to poly for the wildtype, which is zero by definition of the eigenvalues):

Factor[Normal[Series[polym/.λm->λ+ϵ*Δλm/.smallmutant,{ϵ,0,0}]]]

Factor[%-poly]

To order ϵ, we can find effect of the mutation on the rate of growth, Δλm:

Factor[SeriesCoefficient[Series[polym/.λm->λ+ϵ*Δλm/.smallmutant,{ϵ,0,1}],1]];
Factor[Flatten[Solve[%==0,Δλm]]]

The denominator is shared and equals the slope of the characteristic polynomial at the leading eigenvalue (positive given the form of the cubic):

denom=D[poly,λ]//Factor

selection=Collect[denom Δλm/slope/.%%,smallterms,Factor]

We can also understand selection on mutants by using the method in Box 10.2 of Otto and Day, calculating selection as ((dλ)/dz)=(((v^T) ((dM)/dz)u)/((v^T) u)), where z represents some underlying trait that the mutation affects (e.g., replication rate) and (v^T) and u are the left and right eigenvectors of M, the matrix describing movement between the infected classes for the resident:

submat=stabmat[[2;;5,2;;5]];
%//MatrixForm

and for the mutant:

submatm=stabmat[[6;;9,6;;9]];
%//MatrixForm

Calculating the relative reproductive values (the left eigenvector):

left={vE,vP,vI,vA};

left.submat-λ left

leftvector=left/.Flatten[Solve[(%[[2;;4]])==0,Table[left[[i]],{i,2,4}]]]//Simplify

In addition, the fact that the characteristic polynomial is zero guarantees that the first term in v M - v λ  is also zero:

((left.submat-λ left)[[1]]/.Table[left[[i]]->leftvector[[i]],{i,1,4}]//Factor)/.Factor[poly]->0

Calculating the proportion in each class (the right eigenvector):

right={uE,uP,uI,uA};

submat.right-λ right

rightvector=right/.Flatten[Solve[(%[[2;;4]])==0,Table[right[[i]],{i,2,4}]]]//Simplify

In addition, the fact that the characteristic polynomial is zero guarantees that the first term in v M - v λ  is also zero:

((submat.right-λ right)[[1]]/.Table[right[[i]]->rightvector[[i]],{i,1,4}]//Factor)/.Factor[poly]->0

Note that the elements of both left and right eigenvectors are positive in a growing population (λ>0).  Although the length of these eigenvectors is arbitrary, we constrain them so that (v^T) u in the denominator of ((dλ)/dz)=(((v^T) ((dM)/dz)u)/((v^T) u)) is one by rescaling the right eigenvector (choosing this one keeps the left eigenvector easier to read):

rightvector=(rightvector/(leftvector.rightvector));

((dM)/dz) represents the mutational effect on the transition matrix for a viral genotype:

mutmat= Factor[Normal[Series[submatm-submat/.smallmutant,{ϵ,0,1}]]/.ϵ->1];
MatrixForm[%]

This gives us an equivalent way to write the selection coefficient on the mutation, using the fact that the characteristic polynomial is zero:

Factor[leftvector.mutmat.rightvector-(selection/.slope->denom)]/.Factor[poly]->0

Selection is thus proportional to:

Collect[Table[left[[i]],{i,1,4}].mutmat.Table[right[[i]],{i,1,4}],smallterms,Factor]

S vE (uA ΔβA+ uI ΔβI + uP ΔβP)-uI vI Δα+uE ((1-f)vP+f vA-vE) ΔκE+uE (vA-vP) Δf κE-uI vI ΔκI-uP (vP-vI) ΔκP-uA vA ΔκA

This way of calculating selection makes it clearer what selection acts upon when a disease is spreading (λ>0):

	* Always favors increased transmission [ΔβA, ΔβP, ΔβI]
	* Always favors lower virulence [Δα] and slower recovery [ΔκI] 
	* Always favors faster progression from E -> A/P [ΔκE] because (1-f)vP+f vA>vE [see Note 1]
	* Favors slower progression from P -> I [ΔκP] if vP > vI [always true if βI ≤ βP (1+1/(Dλ) - see Note 2]
	* Favors more asymptomatics individuals [Δf] only if vA>vP [true only with high enough βA - see Note 3]
	* Never favors faster progression from A -> R [ΔκAm]

[Note 1] In a spreading disease, (1-f)vP+f vA is always greater than vE because the next stage individuals (P&A)  pass on the disease earlier than exposed individuals, contributing earlier to growth.

Factor[(1-f)leftvector[[2]]+(f)leftvector[[4]]-leftvector[[1]]/.Solve[poly==0,βP]]

[Note 2] vP is greater than vI as long as βI ≤ βP (1+(α+γ)/(λ)) because infected individuals have less time to pass on the disease than pre-asymptomatics.

Factor[leftvector[[3]]-leftvector[[2]]]

Solve[%==0,βI]//Simplify

[Note 3] vP is greater than vA as long as presymptomatics are expected to give rise to enough more cases than asymptomatics that ((βP)/(κP+λ))≥((βA)/(κA+λ)) (symptomatic infectious individuals make this condition even easier to satisfy).

Collect[leftvector[[4]]-leftvector[[2]],{βP,βI,βA},Factor]

Figure 3 - Numerical analyses

Code

Common settings

Clear["Global`*"]

Off[General::spell1]
Off[General::spell]
Off[NDSolve::nlnum]

tfinal=500;

tfig=500;

label={"Time (days)","Frequency"};
label=None;

parvec={β, βm,α, αm,κyP,κyPm,f,fm,κyE, κyI,κyA, c1,c2,c3,mut};

Clear[finalfreq,NSolution]
finalfreq[vec_]:=finalfreq[vec]=Block[{β=vec[[1]], βm=vec[[2]],α=vec[[3]], αm=vec[[4]],κyP=vec[[5]],κyPm=vec[[6]],f=vec[[7]],fm=vec[[8]],κyE=vec[[9]], κyI=vec[[10]],κyA=vec[[11]],c1=vec[[12]],c2=vec[[13]],c3=vec[[14]],mut=vec[[15]]},

(*REMAINING PARAMETERS*)
βyA=β/10;βyP=β;βyI=β/3;
κyAm=κyA;κyIm=κyI;
βyAm=βm/10;βyPm=βm;βyIm=βm/3;

(* E: Exposed (non-infectious) cases *)
(* A: asymptomatic cases *)
(* P: pre-symptomatic cases *)
(* S: symptomatic cases *)
(* f: proportion of fsymptomatic hosts *)
(* c: amount of distancing applied *)

(*INOCULATION*)
start=(10^(-5));

t=.;
c=.;

c[t_]:=Evaluate[If[t<t1min,0,If[t<t1max,c1,If[t<t2min,0,If[t<t2max,c2,If[t<t3min,0,If[t<t3max,c3,0]]]]]]];
h[t]=(1-c[t])(βyA yA[t]+βyP yP[t]+βyI yI[t]);
hm[t]=(1-c[t])(βyAm yAm[t]+βyPm yPm[t]+βyIm yIm[t]);

sys={S'[t]==-(h[t]+hm[t]) S[t],
yE'[t]==h[t] S[t]-κyE yE[t],
yA'[t]==f κyE yE[t]-κyA yA[t],
yP'[t]==(1-f) κyE yE[t]-κyP yP[t],
yI'[t]==κyP yP[t]-(κyI+α) yI[t],
yR'[t]==κyA yA[t]+κyI yI[t],
yD'[t]==α yI[t]};  (*yD serves as a counter of the cumulative number of deaths*)

sysm={
yEm'[t]==hm[t] S[t]-κyE yEm[t],
yAm'[t]==fm κyE yEm[t]-κyAm yAm[t],
yPm'[t]==(1-fm) κyE yEm[t]-κyPm yPm[t],
yIm'[t]==κyPm yPm[t]-(κyIm+αm) yIm[t],
yRm'[t]==κyAm yAm[t]+κyIm yIm[t],
yDm'[t]==αm yIm[t]};  (*yDm serves as a counter of the cumulative number of deaths*)

(*FIRST BOUT OF CONTROL*)
t1min=75;
t1max=150;
(*SECOND BOUT OF CONTROL*)
t2min=225;
t2max=300;
(*THIRD BOUT OF CONTROL*)
t3min=375;
t3max=450;

init={S[0]==1,yE[0]==start(1-mut),yA[0]==0,yP[0]==0,yI[0]==0,yR[0]==0,yD[0]==0};
initm={yEm[0]==start(mut),yAm[0]==0,yPm[0]==0,yIm[0]==0,yRm[0]==0,yDm[0]==0};

var={S,yE,yA,yP,yI,yR,yD};
varm={yEm,yAm,yPm,yIm,yRm,yDm};

NSolution[vec]=
NDSolve[
Flatten[{sys,sysm,init,initm}],Flatten[{var,varm}],{t,0,tfinal},AccuracyGoal->100];

(*R0*)
Print["R0 of the wild type = ",((βyA f)/(κyA))+(((1-f) (βyP (α+κyI)+βyI κyP))/((α+κyI) κyP))];

(*Rm*)
Print["R0 of the mutant type = ",((βyAm f)/(κyAm))+(((1-f) (βyPm (αm+κyIm)+βyIm κyPm))/((αm+κyIm) κyPm))];

(*CASE MORTALITY*)
Print["Case mortality = ",(1-f)(α/(κyI+α))];
Print["Case mortality (mutant) = ",(1-fm)((αm)/(κyI+αm))];

(*CUMULATIVE MORTALITY*)
Print["Cumulative mortality (total) = ",yD[tfinal]+yDm[tfinal]/.Flatten[NSolution[vec]]];
Print["Cumulative mortality (mutant) = ",yDm[tfinal]/.Flatten[NSolution[vec]]];

]

Figure 3a: Evolution of transmission β, no control measures

pars={β->1, βm->1.2, (*TRANSMISSION OF THE MUTANT*)
α->0.005, αm->0.005, (*VIRULENCE OF THE MUTANT*)
κyP->1,κyPm->1, (* 1/TIME IN PRESYMPTOMATIC CLASS FOR THE MUTANT *)

f->0.2, fm->0.2, (* PROPORTION OF ASYMPTOMATIC CASES *)
κyE->0.25, (* 1/TIME IN EXPOSED CLASS *)
κyI->0.2, (* 1/TIME IN SYMPTOMATIC CLASS *)
κyA->0.11, (* 1/TIME IN ASYMPTOMATIC CLASS *)

(*Strength of control measures*)
c1->0,c2->0,c3->0
};

The doublng time of the disease when rare, as calculated from the eigenvalue (“poly”) of the stability matrix:

Sort[Solve[(κE (-S α βP κA+f S α βP κA-S βP κA κI+f S βP κA κI-f S α βA κP+α κA κP-S βI κA κP+f S βI κA κP-f S βA κI κP+κA κI κP)+(-f S α βA κE-S α βP κE+f S α βP κE+α κA κE-S βP κA κE+f S βP κA κE-f S βA κE κI-S βP κE κI+f S βP κE κI+κA κE κI+α κA κP+α κE κP-f S βA κE κP-S βI κE κP+f S βI κE κP+κA κE κP+κA κI κP+κE κI κP) λ+(α κA+α κE-f S βA κE-S βP κE+f S βP κE+κA κE+κA κI+κE κI+α κP+κA κP+κE κP+κI κP) (λ^2)+(α+κA+κE+κI+κP) (λ^3)+(λ^4)/.S->1/.βA->β/10/.βP->β/.βI->β/3/.κA->κyA/.κE->κyE/.κP->κyP/.κI->κyI/.pars)==0,λ]]
Solve[(Exp[λ t]==2)/.Last[%],t]

parset=parvec/.pars/.mut->0;
finalfreq[parset]

tNOevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]

fnoevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tNOevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red,Dashed}},AspectRatio->0.75];

fSDashed=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005],Dashed, Blue},AspectRatio->0.75];

mortalitynoevol=yD[tNOevolNC]+yDm[tNOevolNC]/.Flatten[NSolution[parset]];

parset=parvec/.pars/.mut->0.01;
finalfreq[parset]

tevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]
fevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red}},AspectRatio->0.75];

freqERADIC=Evaluate[(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))/.NSolution[parset]]/.t->tevolNC

FIGfreq=Plot[Evaluate[{(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{0,1}},PlotStyle->{{Thickness[0.008],Black}},AspectRatio->0.75];

fS=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005], Blue},AspectRatio->0.75];

fcontrol=Plot[2*c[t]/.pars,{t,0,tfig},Filling->Axis,FillingStyle->LightGray,PlotStyle->None,Frame->True,PlotRange->{{0,tfig},{0,1}},AspectRatio->0.75,Frame->True];


line1=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC-0.02}],Flatten[{tevolNC+5,freqERADIC+0.02}]}]}];
line2=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC+0.02}],Flatten[{tevolNC+5,freqERADIC-0.02}]}]}];

mortalityevol=yD[tevolNC]+yDm[tevolNC]/.Flatten[NSolution[parset]];
(*mortality=Graphics[{Text[Style[StringJoin[ToString[N[Round[100 mortalitynoevol,(10^(-1))] ]],"% (",ToString[N[Round[100 (mortalityevol),(10^(-1))]]],"%)"],12,FontFamily->"Helvetica"],{(t3min+t3max)/2,0.3}]}];*)

DeathRates=Graphics[{Pink,Rectangle[{t3max-tfig/20,0.7},{t3max,0.7+10*mortalitynoevol}],Red,Rectangle[{t3max,0.7},{t3max+tfig/20,0.7+10*mortalityevol}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.2/20,0.705}],Black,Rectangle[{t3max-tfig 1.4/20,0.8},{t3max-tfig 1.2/20,0.805}],Black,Rectangle[{t3max-tfig 1.4/20,0.9},{t3max-tfig 1.2/20,0.905}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.42/20,0.905}]}];

FigTOP1=Show[fS,FIGfreq,fnoevol,fevol,fSDashed,line1,line2,DeathRates,Frame->True,
TicksStyle->Directive[Black,FontFamily->"Helvetica",11],
LabelStyle->Directive[12,FontFamily->"Helvetica",Black],
FrameTicks->{{N[{0.2,0.4,0.6,0.8}],N[{0.2,0.4,0.6,0.8}]},{{0,100,200,300,400,500},{0,100,200,300,400,500}}},FrameTicksStyle->{{Directive[FontSize->12],Directive[FontOpacity->0,FontSize->0]},{Directive[FontSize->12],Directive[FontOpacity->0,FontSize->0]}}
]

Figure 3b: Evolution of asymptomatic fraction f,  no control measures

pars={β->1, βm->1, (*TRANSMISSION OF THE MUTANT*)
α->0.005, αm->0.005, (*VIRULENCE OF THE MUTANT*)
κyP->1,κyPm->1, (* 1/TIME IN PRESYMPTOMATIC CLASS FOR THE MUTANT *)

f->0.2, fm->0.1, (* PROPORTION OF ASYMPTOMATIC CASES *)
κyE->0.25, (* 1/TIME IN EXPOSED CLASS *)
κyI->0.2, (* 1/TIME IN SYMPTOMATIC CLASS *)
κyA->0.11, (* 1/TIME IN ASYMPTOMATIC CLASS *)

(*Strength of control measures*)
c1->0,c2->0,c3->0
};

parset=parvec/.pars/.mut->0;
finalfreq[parset]

tNOevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]

fnoevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tNOevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red,Dashed}},AspectRatio->0.75];

fSDashed=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005],Dashed, Blue},AspectRatio->0.75];

mortalitynoevol=yD[tNOevolNC]+yDm[tNOevolNC]/.Flatten[NSolution[parset]];

parset=parvec/.pars/.mut->0.01;
finalfreq[parset]

tevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]
fevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red}},AspectRatio->0.75];

freqERADIC=Evaluate[(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))/.NSolution[parset]]/.t->tevolNC

FIGfreq=Plot[Evaluate[{(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{0,1}},PlotStyle->{{Thickness[0.008],Black}},AspectRatio->0.75];

fS=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005], Blue},AspectRatio->0.75];

fcontrol=Plot[2*c[t]/.pars,{t,0,tfig},Filling->Axis,FillingStyle->LightGray,PlotStyle->None,Frame->True,PlotRange->{{0,tfig},{0,1}},AspectRatio->0.75,Frame->True];


line1=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC-0.02}],Flatten[{tevolNC+5,freqERADIC+0.02}]}]}];
line2=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC+0.02}],Flatten[{tevolNC+5,freqERADIC-0.02}]}]}];


mortalityevol=yD[tevolNC]+yDm[tevolNC]/.Flatten[NSolution[parset]];
(*mortality=Graphics[{Text[Style[StringJoin[ToString[N[Round[100 mortalitynoevol,(10^(-1))] ]],"% (",ToString[N[Round[100 (mortalityevol),(10^(-1))]]],"%)"],12,FontFamily->"Helvetica"],{(t3min+t3max)/2,0.3}]}];*)

DeathRates=Graphics[{Pink,Rectangle[{t3max-tfig/20,0.7},{t3max,0.7+10*mortalitynoevol}],Red,Rectangle[{t3max,0.7},{t3max+tfig/20,0.7+10*mortalityevol}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.2/20,0.705}],Black,Rectangle[{t3max-tfig 1.4/20,0.8},{t3max-tfig 1.2/20,0.805}],Black,Rectangle[{t3max-tfig 1.4/20,0.9},{t3max-tfig 1.2/20,0.905}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.42/20,0.905}]}];

FigTOP2=Show[fS,FIGfreq,fnoevol,fevol,fSDashed,line1,line2,DeathRates,Frame->True,
TicksStyle->Directive[Black,FontFamily->"Helvetica",11],
LabelStyle->Directive[12,FontFamily->"Helvetica",Black],FrameTicks->{{N[{0.2,0.4,0.6,0.8}],N[{0.2,0.4,0.6,0.8}]},{{0,100,200,300,400,500},{0,100,200,300,400,500}}},FrameTicksStyle->{{Directive[FontOpacity->0,FontSize->0],Directive[FontOpacity->0,FontSize->0]},{Directive[FontSize->12],Directive[FontOpacity->0,FontSize->0]}}
]

Figure 3c: Evolution of pre-symptomatic phase κ, no control measures

pars={β->1, βm->1, (*TRANSMISSION OF THE MUTANT*)
α->0.005, αm->0.005, (*VIRULENCE OF THE MUTANT*)
κyP->1,κyPm->2/3, (* 1/TIME IN PRESYMPTOMATIC CLASS FOR THE MUTANT *)

f->0.2, fm->0.2, (* PROPORTION OF ASYMPTOMATIC CASES *)
κyE->0.25, (* 1/TIME IN EXPOSED CLASS *)
κyI->0.2, (* 1/TIME IN SYMPTOMATIC CLASS *)
κyA->0.11, (* 1/TIME IN ASYMPTOMATIC CLASS *)

(*Strength of control measures*)
c1->0,c2->0,c3->0
};

parset=parvec/.pars/.mut->0;
finalfreq[parset]

tNOevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]

fnoevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tNOevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red,Dashed}},AspectRatio->0.75];

fSDashed=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005],Dashed, Blue},AspectRatio->0.75];

mortalitynoevol=yD[tNOevolNC]+yDm[tNOevolNC]/.Flatten[NSolution[parset]];

parset=parvec/.pars/.mut->0.01;
finalfreq[parset]

tevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]
fevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red}},AspectRatio->0.75];

freqERADIC=Evaluate[(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))/.NSolution[parset]]/.t->tevolNC

FIGfreq=Plot[Evaluate[{(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{0,1}},PlotStyle->{{Thickness[0.008],Black}},AspectRatio->0.75];

fS=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005], Blue},AspectRatio->0.75];

fcontrol=Plot[2*c[t]/.pars,{t,0,tfig},Filling->Axis,FillingStyle->LightGray,PlotStyle->None,Frame->True,PlotRange->{{0,tfig},{0,1}},AspectRatio->0.75,Frame->True];


line1=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC-0.02}],Flatten[{tevolNC+5,freqERADIC+0.02}]}]}];
line2=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC+0.02}],Flatten[{tevolNC+5,freqERADIC-0.02}]}]}];


mortalityevol=yD[tevolNC]+yDm[tevolNC]/.Flatten[NSolution[parset]];
(*mortality=Graphics[{Text[Style[StringJoin[ToString[N[Round[100 mortalitynoevol,(10^(-1))] ]],"% (",ToString[N[Round[100 (mortalityevol),(10^(-1))]]],"%)"],12,FontFamily->"Helvetica"],{(t3min+t3max)/2,0.3}]}];*)

DeathRates=Graphics[{Pink,Rectangle[{t3max-tfig/20,0.7},{t3max,0.7+10*mortalitynoevol}],Red,Rectangle[{t3max,0.7},{t3max+tfig/20,0.7+10*mortalityevol}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.2/20,0.705}],Black,Rectangle[{t3max-tfig 1.4/20,0.8},{t3max-tfig 1.2/20,0.805}],Black,Rectangle[{t3max-tfig 1.4/20,0.9},{t3max-tfig 1.2/20,0.905}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.42/20,0.905}]}];

FigTOP3=Show[fS,FIGfreq,fnoevol,fevol,fSDashed,line1,line2,DeathRates,Frame->True,
TicksStyle->Directive[Black,FontFamily->"Helvetica",11],
LabelStyle->Directive[12,FontFamily->"Helvetica",Black],FrameTicks->{{N[{0.2,0.4,0.6,0.8}],N[{0.2,0.4,0.6,0.8}]},{{0,100,200,300,400,500},{0,100,200,300,400,500}}},FrameTicksStyle->{{Directive[FontOpacity->0,FontSize->0],Directive[FontOpacity->0,FontSize->0]},{Directive[FontSize->12],Directive[FontOpacity->0,FontSize->0]}}
]

Figure 3d: Evolution of virulence α, no control measures

pars={β->1, βm->1, (*TRANSMISSION OF THE MUTANT*)
α->0.005, αm->0.0, (*VIRULENCE OF THE MUTANT*)
κyP->1,κyPm->1, (* 1/TIME IN PRESYMPTOMATIC CLASS FOR THE MUTANT *)

f->0.2, fm->0.2, (* PROPORTION OF ASYMPTOMATIC CASES *)
κyE->0.25, (* 1/TIME IN EXPOSED CLASS *)
κyI->0.2, (* 1/TIME IN SYMPTOMATIC CLASS *)
κyA->0.11, (* 1/TIME IN ASYMPTOMATIC CLASS *)

(*Strength of control measures*)
c1->0,c2->0,c3->0
};

parset=parvec/.pars/.mut->0;
finalfreq[parset]

tNOevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]

fnoevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tNOevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red,Dashed}},AspectRatio->0.75];

fSDashed=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005],Dashed, Blue},AspectRatio->0.75];

mortalitynoevol=yD[tNOevolNC]+yDm[tNOevolNC]/.Flatten[NSolution[parset]];

parset=parvec/.pars/.mut->0.01;
finalfreq[parset]

tevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]
fevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red}},AspectRatio->0.75];

freqERADIC=Evaluate[(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))/.NSolution[parset]]/.t->tevolNC

FIGfreq=Plot[Evaluate[{(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{0,1}},PlotStyle->{{Thickness[0.008],Black}},AspectRatio->0.75];

fS=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005], Blue},AspectRatio->0.75];

fcontrol=Plot[2*c[t]/.pars,{t,0,tfig},Filling->Axis,FillingStyle->LightGray,PlotStyle->None,Frame->True,PlotRange->{{0,tfig},{0,1}},AspectRatio->0.75,Frame->True];


line1=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC-0.02}],Flatten[{tevolNC+5,freqERADIC+0.02}]}]}];
line2=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC+0.02}],Flatten[{tevolNC+5,freqERADIC-0.02}]}]}];


mortalityevol=yD[tevolNC]+yDm[tevolNC]/.Flatten[NSolution[parset]];
(*mortality=Graphics[{Text[Style[StringJoin[ToString[N[Round[100 mortalitynoevol,(10^(-1))] ]],"% (",ToString[N[Round[100 (mortalityevol),(10^(-1))]]],"%)"],12,FontFamily->"Helvetica"],{(t3min+t3max)/2,0.3}]}];*)

DeathRates=Graphics[{Pink,Rectangle[{t3max-tfig/20,0.7},{t3max,0.7+10*mortalitynoevol}],Red,Rectangle[{t3max,0.7},{t3max+tfig/20,0.7+10*mortalityevol}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.2/20,0.705}],Black,Rectangle[{t3max-tfig 1.4/20,0.8},{t3max-tfig 1.2/20,0.805}],Black,Rectangle[{t3max-tfig 1.4/20,0.9},{t3max-tfig 1.2/20,0.905}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.42/20,0.905}]}];

FigTOP4=Show[fS,FIGfreq,fnoevol,fevol,fSDashed,line1,line2,DeathRates,Frame->True,
TicksStyle->Directive[Black,FontFamily->"Helvetica",11],
LabelStyle->Directive[12,FontFamily->"Helvetica",Black],FrameTicks->{{N[{0.2,0.4,0.6,0.8}],N[{0.2,0.4,0.6,0.8}]},{{0,100,200,300,400,500},{0,100,200,300,400,500}}},FrameTicksStyle->{{Directive[FontOpacity->0,FontSize->0],Directive[FontOpacity->0,FontSize->0]},{Directive[FontSize->12],Directive[FontOpacity->0,FontSize->0]}}
]

Figure 3e: Evolution of transmission β, control measures

pars={β->1, βm->1.2, (*TRANSMISSION OF THE MUTANT*)
α->0.005, αm->0.005, (*VIRULENCE OF THE MUTANT*)
κyP->1,κyPm->1, (* 1/TIME IN PRESYMPTOMATIC CLASS FOR THE MUTANT *)

f->0.2, fm->0.2, (* PROPORTION OF ASYMPTOMATIC CASES *)
κyE->0.25, (* 1/TIME IN EXPOSED CLASS *)
κyI->0.2, (* 1/TIME IN SYMPTOMATIC CLASS *)
κyA->0.11, (* 1/TIME IN ASYMPTOMATIC CLASS *)

(*Strength of control measures*)
c1->0.6,c2->0.6,c3->0.6
};

parset=parvec/.pars/.mut->0;
finalfreq[parset]

tNOevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]

fnoevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tNOevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red,Dashed}},AspectRatio->0.75];

fSDashed=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005],Dashed, Blue},AspectRatio->0.75];

mortalitynoevol=yD[tNOevolNC]+yDm[tNOevolNC]/.Flatten[NSolution[parset]];

parset=parvec/.pars/.mut->0.01;
finalfreq[parset]

tevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]
fevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red}},AspectRatio->0.75];

freqERADIC=Evaluate[(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))/.NSolution[parset]]/.t->tevolNC

FIGfreq=Plot[Evaluate[{(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{0,1}},PlotStyle->{{Thickness[0.008],Black}},AspectRatio->0.75];

fS=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005], Blue},AspectRatio->0.75];

fcontrol=Plot[2*c[t]/.pars,{t,0,tfig},Filling->Axis,FillingStyle->LightGray,PlotStyle->None,Frame->True,PlotRange->{{0,tfig},{0,1}},AspectRatio->0.75,Frame->True];

line1=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC-0.02}],Flatten[{tevolNC+5,freqERADIC+0.02}]}]}];
line2=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC+0.02}],Flatten[{tevolNC+5,freqERADIC-0.02}]}]}];

mortalityevol=yD[tevolNC]+yDm[tevolNC]/.Flatten[NSolution[parset]];
(*mortality=Graphics[{Text[Style[StringJoin[ToString[N[Round[100 mortalitynoevol,(10^(-1))] ]],"% (",ToString[N[Round[100 (mortalityevol),(10^(-1))]]],"%)"],12,FontFamily->"Helvetica"],{(t3min+t3max)/2,0.3}]}];*)

DeathRates=Graphics[{Pink,Rectangle[{t3max-tfig/20,0.7},{t3max,0.7+10*mortalitynoevol}],Red,Rectangle[{t3max,0.7},{t3max+tfig/20,0.7+10*mortalityevol}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.2/20,0.705}],Black,Rectangle[{t3max-tfig 1.4/20,0.8},{t3max-tfig 1.2/20,0.805}],Black,Rectangle[{t3max-tfig 1.4/20,0.9},{t3max-tfig 1.2/20,0.905}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.42/20,0.905}]}];

FigBOT1=Show[fcontrol,fS,fSDashed, FIGfreq,fnoevol,fevol,line1,line2,DeathRates,Frame->True,
TicksStyle->Directive[Black,FontFamily->"Helvetica",11],
LabelStyle->Directive[12,FontFamily->"Helvetica",Black],FrameTicks->{{N[{0.2,0.4,0.6,0.8}],N[{0.2,0.4,0.6,0.8}]},{{0,100,200,300,400,500},{0,100,200,300,400,500}}},FrameTicksStyle->{{Directive[FontSize->12],Directive[FontOpacity->0,FontSize->0]},{Directive[FontSize->12],Directive[FontOpacity->0,FontSize->0]}}
]

Figure 3f: Evolution of asymptomatic fraction f, control measures

pars={β->1, βm->1, (*TRANSMISSION OF THE MUTANT*)
α->0.005, αm->0.005, (*VIRULENCE OF THE MUTANT*)
κyP->1,κyPm->1, (* 1/TIME IN PRESYMPTOMATIC CLASS FOR THE MUTANT *)

f->0.2, fm->0.1, (* PROPORTION OF ASYMPTOMATIC CASES *)
κyE->0.25, (* 1/TIME IN EXPOSED CLASS *)
κyI->0.2, (* 1/TIME IN SYMPTOMATIC CLASS *)
κyA->0.11, (* 1/TIME IN ASYMPTOMATIC CLASS *)

(*Strength of control measures*)
c1->0.6,c2->0.6,c3->0.6
};

parset=parvec/.pars/.mut->0;
finalfreq[parset]

tNOevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]

fnoevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tNOevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red,Dashed}},AspectRatio->0.75];

fSDashed=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005],Dashed, Blue},AspectRatio->0.75];

mortalitynoevol=yD[tNOevolNC]+yDm[tNOevolNC]/.Flatten[NSolution[parset]];

parset=parvec/.pars/.mut->0.01;
finalfreq[parset]

tevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]
fevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red}},AspectRatio->0.75];

freqERADIC=Evaluate[(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))/.NSolution[parset]]/.t->tevolNC

FIGfreq=Plot[Evaluate[{(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{0,1}},PlotStyle->{{Thickness[0.008],Black}},AspectRatio->0.75];

fS=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005], Blue},AspectRatio->0.75];

fcontrol=Plot[2*c[t]/.pars,{t,0,tfig},Filling->Axis,FillingStyle->LightGray,PlotStyle->None,Frame->True,PlotRange->{{0,tfig},{0,1}},AspectRatio->0.75,Frame->True];

line1=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC-0.02}],Flatten[{tevolNC+5,freqERADIC+0.02}]}]}];
line2=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC+0.02}],Flatten[{tevolNC+5,freqERADIC-0.02}]}]}];

mortalityevol=yD[tevolNC]+yDm[tevolNC]/.Flatten[NSolution[parset]];
(*mortality=Graphics[{Text[Style[StringJoin[ToString[N[Round[100 mortalitynoevol,(10^(-1))] ]],"% (",ToString[N[Round[100 (mortalityevol),(10^(-1))]]],"%)"],12,FontFamily->"Helvetica"],{(t3min+t3max)/2,0.3}]}];*)

DeathRates=Graphics[{Pink,Rectangle[{t3max-tfig/20,0.7},{t3max,0.7+10*mortalitynoevol}],Red,Rectangle[{t3max,0.7},{t3max+tfig/20,0.7+10*mortalityevol}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.2/20,0.705}],Black,Rectangle[{t3max-tfig 1.4/20,0.8},{t3max-tfig 1.2/20,0.805}],Black,Rectangle[{t3max-tfig 1.4/20,0.9},{t3max-tfig 1.2/20,0.905}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.42/20,0.905}]}];

FigBOT2=Show[fcontrol,fS,FIGfreq,fnoevol,fevol,fSDashed,line1,line2,DeathRates,Frame->True,
TicksStyle->Directive[Black,FontFamily->"Helvetica",11],
LabelStyle->Directive[12,FontFamily->"Helvetica",Black],FrameTicks->{{N[{0.2,0.4,0.6,0.8}],N[{0.2,0.4,0.6,0.8}]},{{0,100,200,300,400,500},{0,100,200,300,400,500}}},FrameTicksStyle->{{Directive[FontOpacity->0,FontSize->0],Directive[FontOpacity->0,FontSize->0]},{Directive[FontSize->12],Directive[FontOpacity->0,FontSize->0]}}
]

Figure 3g: Evolution of pre-symptomatic phase κ, control measures

pars={β->1, βm->1, (*TRANSMISSION OF THE MUTANT*)
α->0.005, αm->0.005, (*VIRULENCE OF THE MUTANT*)
κyP->1,κyPm->2/3, (* 1/TIME IN PRESYMPTOMATIC CLASS FOR THE MUTANT *)

f->0.2, fm->0.2, (* PROPORTION OF ASYMPTOMATIC CASES *)
κyE->0.25, (* 1/TIME IN EXPOSED CLASS *)
κyI->0.2, (* 1/TIME IN SYMPTOMATIC CLASS *)
κyA->0.11, (* 1/TIME IN ASYMPTOMATIC CLASS *)

(*Strength of control measures*)
c1->0.6,c2->0.6,c3->0.6
};

parset=parvec/.pars/.mut->0;
finalfreq[parset]

tNOevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]

fnoevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tNOevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red,Dashed}},AspectRatio->0.75];

fSDashed=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005],Dashed, Blue},AspectRatio->0.75];

mortalitynoevol=yD[tNOevolNC]+yDm[tNOevolNC]/.Flatten[NSolution[parset]];

parset=parvec/.pars/.mut->0.01;
finalfreq[parset]

tevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]
fevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red}},AspectRatio->0.75];

freqERADIC=Evaluate[(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))/.NSolution[parset]]/.t->tevolNC

FIGfreq=Plot[Evaluate[{(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{0,1}},PlotStyle->{{Thickness[0.008],Black}},AspectRatio->0.75];

fS=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005], Blue},AspectRatio->0.75];

fcontrol=Plot[2*c[t]/.pars,{t,0,tfig},Filling->Axis,FillingStyle->LightGray,PlotStyle->None,Frame->True,PlotRange->{{0,tfig},{0,1}},AspectRatio->0.75,Frame->True];

line1=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC-0.02}],Flatten[{tevolNC+5,freqERADIC+0.02}]}]}];
line2=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC+0.02}],Flatten[{tevolNC+5,freqERADIC-0.02}]}]}];

mortalityevol=yD[tevolNC]+yDm[tevolNC]/.Flatten[NSolution[parset]];
(*mortality=Graphics[{Text[Style[StringJoin[ToString[N[Round[100 mortalitynoevol,(10^(-1))] ]],"% (",ToString[N[Round[100 (mortalityevol),(10^(-1))]]],"%)"],12,FontFamily->"Helvetica"],{(t3min+t3max)/2,0.3}]}];*)

DeathRates=Graphics[{Pink,Rectangle[{t3max-tfig/20,0.7},{t3max,0.7+10*mortalitynoevol}],Red,Rectangle[{t3max,0.7},{t3max+tfig/20,0.7+10*mortalityevol}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.2/20,0.705}],Black,Rectangle[{t3max-tfig 1.4/20,0.8},{t3max-tfig 1.2/20,0.805}],Black,Rectangle[{t3max-tfig 1.4/20,0.9},{t3max-tfig 1.2/20,0.905}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.42/20,0.905}]}];

FigBOT3=Show[fcontrol,fS,FIGfreq,fnoevol,fevol,fSDashed,line1,line2,DeathRates,Frame->True,
TicksStyle->Directive[Black,FontFamily->"Helvetica",11],
LabelStyle->Directive[12,FontFamily->"Helvetica",Black],FrameTicks->{{N[{0.2,0.4,0.6,0.8}],N[{0.2,0.4,0.6,0.8}]},{{0,100,200,300,400,500},{0,100,200,300,400,500}}},FrameTicksStyle->{{Directive[FontOpacity->0,FontSize->0],Directive[FontOpacity->0,FontSize->0]},{Directive[FontSize->12],Directive[FontOpacity->0,FontSize->0]}}
]

Figure 3h: Evolution of virulence α, control measures

pars={β->1, βm->1, (*TRANSMISSION OF THE MUTANT*)
α->0.005, αm->0.00, (*VIRULENCE OF THE MUTANT*)
κyP->1,κyPm->1, (* 1/TIME IN PRESYMPTOMATIC CLASS FOR THE MUTANT *)

f->0.2, fm->0.2, (* PROPORTION OF ASYMPTOMATIC CASES *)
κyE->0.25, (* 1/TIME IN EXPOSED CLASS *)
κyI->0.2, (* 1/TIME IN SYMPTOMATIC CLASS *)
κyA->0.11, (* 1/TIME IN ASYMPTOMATIC CLASS *)

(*Strength of control measures*)
c1->0.6,c2->0.6,c3->0.6
};

parset=parvec/.pars/.mut->0;
finalfreq[parset]

tNOevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]

fnoevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tNOevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red,Dashed}},AspectRatio->0.75];

fSDashed=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005],Dashed, Blue},AspectRatio->0.75];

mortalitynoevol=yD[tNOevolNC]+yDm[tNOevolNC]/.Flatten[NSolution[parset]];

parset=parvec/.pars/.mut->0.01;
finalfreq[parset]

tevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]
fevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red}},AspectRatio->0.75];

freqERADIC=Evaluate[(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))/.NSolution[parset]]/.t->tevolNC

FIGfreq=Plot[Evaluate[{(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{0,1}},PlotStyle->{{Thickness[0.008],Black}},AspectRatio->0.75];

fS=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005], Blue},AspectRatio->0.75];

fcontrol=Plot[2*c[t]/.pars,{t,0,tfig},Filling->Axis,FillingStyle->LightGray,PlotStyle->None,Frame->True,PlotRange->{{0,tfig},{0,1}},AspectRatio->0.75,Frame->True];

line1=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC-0.02}],Flatten[{tevolNC+5,freqERADIC+0.02}]}]}];
line2=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC+0.02}],Flatten[{tevolNC+5,freqERADIC-0.02}]}]}];

mortalityevol=yD[tevolNC]+yDm[tevolNC]/.Flatten[NSolution[parset]];
(*mortality=Graphics[{Text[Style[StringJoin[ToString[N[Round[100 mortalitynoevol,(10^(-1))] ]],"% (",ToString[N[Round[100 (mortalityevol),(10^(-1))]]],"%)"],12,FontFamily->"Helvetica"],{(t3min+t3max)/2,0.3}]}];*)

DeathRates=Graphics[{Pink,Rectangle[{t3max-tfig/20,0.7},{t3max,0.7+10*mortalitynoevol}],Red,Rectangle[{t3max,0.7},{t3max+tfig/20,0.7+10*mortalityevol}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.2/20,0.705}],Black,Rectangle[{t3max-tfig 1.4/20,0.8},{t3max-tfig 1.2/20,0.805}],Black,Rectangle[{t3max-tfig 1.4/20,0.9},{t3max-tfig 1.2/20,0.905}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.42/20,0.905}]}];

FigBOT4=Show[fcontrol,fS,FIGfreq,fnoevol,fevol,fSDashed,line1,line2,DeathRates,Frame->True,
TicksStyle->Directive[Black,FontFamily->"Helvetica",11],
LabelStyle->Directive[12,FontFamily->"Helvetica",Black],FrameTicks->{{N[{0.2,0.4,0.6,0.8}],N[{0.2,0.4,0.6,0.8}]},{{0,100,200,300,400,500},{0,100,200,300,400,500}}},FrameTicksStyle->{{Directive[FontOpacity->0,FontSize->0],Directive[FontOpacity->0,FontSize->0]},{Directive[FontSize->12],Directive[FontOpacity->0,FontSize->0]}}
]

Altogether

Show[
GraphicsGrid[{{FigTOP1,FigTOP2,FigTOP3,FigTOP4},{FigBOT1,FigBOT2,FigBOT3,FigBOT4}},FrameLabel->label,Spacings->{-10,10}],
Graphics[Text[Style["(a)",16,FontFamily->"Helvetica"],{50,-33}]],
Graphics[Text[Style["(b)",16,FontFamily->"Helvetica"],{384,-33}]],
Graphics[Text[Style["(c)",16,FontFamily->"Helvetica"],{734,-33}]],
Graphics[Text[Style["(d)",16,FontFamily->"Helvetica"],{1084,-33}]],
Graphics[Text[Style["(e)",16,FontFamily->"Helvetica"],{50,-315}]],
Graphics[Text[Style["(f)",16,FontFamily->"Helvetica"],{384,-315}]],
Graphics[Text[Style["(g)",16,FontFamily->"Helvetica"],{734,-315}]],
Graphics[Text[Style["(h)",16,FontFamily->"Helvetica"],{1084,-315}]],
Graphics[Text[Style["Transmission",16,FontFamily->"Helvetica"],{180,10}]],
Graphics[Text[Style["Asymptomatic fraction",16,FontFamily->"Helvetica"],{ 520,10}]],
Graphics[Text[Style["Pre-symptomatic phase",16,FontFamily->"Helvetica"],{870,10}]],
Graphics[Text[Style["Virulence",16,FontFamily->"Helvetica"],{1220,10}]],
Graphics[Text[Style["Time (days)",16,FontFamily->"Helvetica"],{700,-580}]],
Graphics[Rotate[Text[Style["Frequency",16,FontFamily->"Helvetica"],{-20,-270}],Pi/2]]
]

Export["Fig3.pdf",%];

Figure 4 - Numerical analyses

Code

Common settings

Clear["Global`*"]

Off[General::spell1]
Off[General::spell]
Off[NDSolve::nlnum]

tfinal=500;

tfig=500;

label={"Time (days)","Frequency"};
label=None;

parvec={β, βm,α, αm,κyP,κyPm,f,fm,κyE, κyI,κyA, c1,c2,c3,mut};

Clear[finalfreq,NSolution]
finalfreq[vec_]:=finalfreq[vec]=Block[{β=vec[[1]], βm=vec[[2]],α=vec[[3]], αm=vec[[4]],κyP=vec[[5]],κyPm=vec[[6]],f=vec[[7]],fm=vec[[8]],κyE=vec[[9]], κyI=vec[[10]],κyA=vec[[11]],c1=vec[[12]],c2=vec[[13]],c3=vec[[14]],mut=vec[[15]]},

(*REMAINING PARAMETERS*)
βyA=β/10;βyP=β;βyI=β/3;
κyAm=κyA;κyIm=κyI;
βyAm=βm/10;βyPm=βm;βyIm=βm/3;

(* E: Exposed (non-infectious) cases *)
(* A: asymptomatic cases *)
(* P: pre-symptomatic cases *)
(* S: symptomatic cases *)
(* f: proportion of fsymptomatic hosts *)
(* c: amount of distancing applied *)

(*INOCULATION*)
start=(10^(-5));

t=.;
c=.;

c[t_]:=Evaluate[If[t<t1min,0,If[t<t1max,c1,If[t<t2min,0,If[t<t2max,c2,If[t<t3min,0,If[t<t3max,c3,0]]]]]]];
h[t]=(1-c[t])(βyA yA[t]+βyP yP[t]+βyI yI[t]);
hm[t]=(1-c[t])(βyAm yAm[t]+βyPm yPm[t]+βyIm yIm[t]);

sys={S'[t]==-(h[t]+hm[t]) S[t],
yE'[t]==h[t] S[t]-κyE yE[t],
yA'[t]==f κyE yE[t]-κyA yA[t],
yP'[t]==(1-f) κyE yE[t]-κyP yP[t],
yI'[t]==κyP yP[t]-(κyI+α) yI[t],
yR'[t]==κyA yA[t]+κyI yI[t],
yD'[t]==α yI[t]};  (*yD serves as a counter of the cumulative number of deaths*)

sysm={
yEm'[t]==hm[t] S[t]-κyE yEm[t],
yAm'[t]==fm κyE yEm[t]-κyAm yAm[t],
yPm'[t]==(1-fm) κyE yEm[t]-κyPm yPm[t],
yIm'[t]==κyPm yPm[t]-(κyIm+αm) yIm[t],
yRm'[t]==κyAm yAm[t]+κyIm yIm[t],
yDm'[t]==αm yIm[t]};  (*yDm serves as a counter of the cumulative number of deaths*)

(*FIRST BOUT OF CONTROL*)
t1min=75;
t1max=150;
(*SECOND BOUT OF CONTROL*)
t2min=225;
t2max=300;
(*THIRD BOUT OF CONTROL*)
t3min=375;
t3max=450;

init={S[0]==1,yE[0]==start(1-mut),yA[0]==0,yP[0]==0,yI[0]==0,yR[0]==0,yD[0]==0};
initm={yEm[0]==start(mut),yAm[0]==0,yPm[0]==0,yIm[0]==0,yRm[0]==0,yDm[0]==0};

var={S,yE,yA,yP,yI,yR,yD};
varm={yEm,yAm,yPm,yIm,yRm,yDm};

NSolution[vec]=
NDSolve[
Flatten[{sys,sysm,init,initm}],Flatten[{var,varm}],{t,0,tfinal},AccuracyGoal->100];

(*R0*)
Print["R0 of the wild type = ",((βyA f)/(κyA))+(((1-f) (βyP (α+κyI)+βyI κyP))/((α+κyI) κyP))];

(*Rm*)
Print["R0 of the mutant type = ",((βyAm f)/(κyAm))+(((1-f) (βyPm (αm+κyIm)+βyIm κyPm))/((αm+κyIm) κyPm))];

(*CASE MORTALITY*)
Print["Case mortality = ",(1-f)(α/(κyI+α))];
Print["Case mortality (mutant) = ",(1-fm)((αm)/(κyI+αm))];

(*CUMULATIVE MORTALITY*)
Print["Cumulative mortality (total) = ",yD[tfinal]+yDm[tfinal]/.Flatten[NSolution[vec]]];
Print["Cumulative mortality (mutant) = ",yDm[tfinal]/.Flatten[NSolution[vec]]];

]

Figure 4a: Evolution with positive pleiotropy between transmission β and virulence α, no control measures

pars={β->1, βm->1.2, (*TRANSMISSION OF THE MUTANT*)
α->0.005, αm->0.01, (*VIRULENCE OF THE MUTANT*)
κyP->1,κyPm->1, (* 1/TIME IN PRESYMPTOMATIC CLASS FOR THE MUTANT *)

f->0.2, fm->0.2, (* PROPORTION OF ASYMPTOMATIC CASES *)
κyE->0.25, (* 1/TIME IN EXPOSED CLASS *)
κyI->0.2, (* 1/TIME IN SYMPTOMATIC CLASS *)
κyA->0.11, (* 1/TIME IN ASYMPTOMATIC CLASS *)

(*Strength of control measures*)
c1->0,c2->0,c3->0
};
Print["Case mortality = ",(α/(κyI+α))/.pars];
Print["Case mortality (mutant) = ",((αm)/(κyI+αm))/.pars];

parset=parvec/.pars/.mut->0;
finalfreq[parset]

tNOevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]

fnoevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tNOevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red,Dashed}},AspectRatio->0.75];

fSDashed=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005],Dashed, Blue},AspectRatio->0.75];

mortalitynoevol=yD[tNOevolNC]+yDm[tNOevolNC]/.Flatten[NSolution[parset]];

parset=parvec/.pars/.mut->0.01;
finalfreq[parset]

tevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]
fevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red}},AspectRatio->0.75];

freqERADIC=Evaluate[(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))/.NSolution[parset]]/.t->tevolNC

FIGfreq=Plot[Evaluate[{(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{0,1}},PlotStyle->{{Thickness[0.008],Black}},AspectRatio->0.75];

fS=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005], Blue},AspectRatio->0.75];

fcontrol=Plot[2*c[t]/.pars,{t,0,tfig},Filling->Axis,FillingStyle->LightGray,PlotStyle->None,Frame->True,PlotRange->{{0,tfig},{0,1}},AspectRatio->0.75,Frame->True];


line1=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC-0.02}],Flatten[{tevolNC+5,freqERADIC+0.02}]}]}];
line2=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC+0.02}],Flatten[{tevolNC+5,freqERADIC-0.02}]}]}];

mortalityevol=yD[tevolNC]+yDm[tevolNC]/.Flatten[NSolution[parset]];
(*mortality=Graphics[{Text[Style[StringJoin[ToString[N[Round[100 mortalitynoevol,(10^(-1))] ]],"% (",ToString[N[Round[100 (mortalityevol),(10^(-1))]]],"%)"],12,FontFamily->"Helvetica"],{(t3min+t3max)/2,0.3}]}];*)

DeathRates=Graphics[{Pink,Rectangle[{t3max-tfig/20,0.7},{t3max,0.7+10*mortalitynoevol}],Red,Rectangle[{t3max,0.7},{t3max+tfig/20,0.7+10*mortalityevol}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.2/20,0.705}],Black,Rectangle[{t3max-tfig 1.4/20,0.8},{t3max-tfig 1.2/20,0.805}],Black,Rectangle[{t3max-tfig 1.4/20,0.9},{t3max-tfig 1.2/20,0.905}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.42/20,0.905}]}];

FigTOP1=Show[fS,FIGfreq,fnoevol,fevol,fSDashed,line1,line2,DeathRates,Frame->True,
TicksStyle->Directive[Black,FontFamily->"Helvetica",11],
LabelStyle->Directive[12,FontFamily->"Helvetica",Black],FrameTicks->{{N[{0,0.2,0.4,0.6,0.8,1}],None},{{0,100,200,300,400,500},None}},FrameLabel->{None,None}]

Figure 4b: Evolution with negative pleiotropy between transmission β and virulence α, no control measures

pars={β->1, βm->1.2, (*TRANSMISSION OF THE MUTANT*)
α->0.005, αm->0.000, (*VIRULENCE OF THE MUTANT*)
κyP->1,κyPm->1, (* 1/TIME IN PRESYMPTOMATIC CLASS FOR THE MUTANT *)

f->0.2, fm->0.2, (* PROPORTION OF ASYMPTOMATIC CASES *)
κyE->0.25, (* 1/TIME IN EXPOSED CLASS *)
κyI->0.2, (* 1/TIME IN SYMPTOMATIC CLASS *)
κyA->0.11, (* 1/TIME IN ASYMPTOMATIC CLASS *)

(*Strength of control measures*)
c1->0,c2->0,c3->0
};
Print["Case mortality = ",(α/(κyI+α))/.pars];
Print["Case mortality (mutant) = ",((αm)/(κyI+αm))/.pars];

parset=parvec/.pars/.mut->0;
finalfreq[parset]

tNOevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]

fnoevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tNOevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red,Dashed}},AspectRatio->0.75];

fSDashed=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005],Dashed, Blue},AspectRatio->0.75];

mortalitynoevol=yD[tNOevolNC]+yDm[tNOevolNC]/.Flatten[NSolution[parset]];

parset=parvec/.pars/.mut->0.01;
finalfreq[parset]

tevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]
fevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red}},AspectRatio->0.75];

freqERADIC=Evaluate[(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))/.NSolution[parset]]/.t->tevolNC

FIGfreq=Plot[Evaluate[{(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{0,1}},PlotStyle->{{Thickness[0.008],Black}},AspectRatio->0.75];

fS=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005], Blue},AspectRatio->0.75];

fcontrol=Plot[2*c[t]/.pars,{t,0,tfig},Filling->Axis,FillingStyle->LightGray,PlotStyle->None,Frame->True,PlotRange->{{0,tfig},{0,1}},AspectRatio->0.75,Frame->True];


line1=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC-0.02}],Flatten[{tevolNC+5,freqERADIC+0.02}]}]}];
line2=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC+0.02}],Flatten[{tevolNC+5,freqERADIC-0.02}]}]}];

mortalityevol=yD[tevolNC]+yDm[tevolNC]/.Flatten[NSolution[parset]];
(*mortality=Graphics[{Text[Style[StringJoin[ToString[N[Round[100 mortalitynoevol,(10^(-1))] ]],"% (",ToString[N[Round[100 (mortalityevol),(10^(-1))]]],"%)"],12,FontFamily->"Helvetica"],{(t3min+t3max)/2,0.3}]}];*)

DeathRates=Graphics[{Pink,Rectangle[{t3max-tfig/20,0.7},{t3max,0.7+10*mortalitynoevol}],Red,Rectangle[{t3max,0.7},{t3max+tfig/20,0.7+10*mortalityevol}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.2/20,0.705}],Black,Rectangle[{t3max-tfig 1.4/20,0.8},{t3max-tfig 1.2/20,0.805}],Black,Rectangle[{t3max-tfig 1.4/20,0.9},{t3max-tfig 1.2/20,0.905}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.42/20,0.905}]}];

FigTOP2=Show[fS,FIGfreq,fnoevol,fevol,fSDashed,line1,line2,DeathRates,Frame->True,
TicksStyle->Directive[Black,FontFamily->"Helvetica",11],
LabelStyle->Directive[12,FontFamily->"Helvetica",Black],FrameTicks->{{N[{0,0.2,0.4,0.6,0.8,1}],None},{{0,100,200,300,400,500},None}},FrameLabel->{None,None}]

Figure 4c: Evolution with positive pleiotropy between transmission β and virulence α, control measures

pars={β->1, βm->1.2, (*TRANSMISSION OF THE MUTANT*)
α->0.005, αm->0.01, (*VIRULENCE OF THE MUTANT*)
κyP->1,κyPm->1, (* 1/TIME IN PRESYMPTOMATIC CLASS FOR THE MUTANT *)

f->0.2, fm->0.2, (* PROPORTION OF ASYMPTOMATIC CASES *)
κyE->0.25, (* 1/TIME IN EXPOSED CLASS *)
κyI->0.2, (* 1/TIME IN SYMPTOMATIC CLASS *)
κyA->0.11, (* 1/TIME IN ASYMPTOMATIC CLASS *)

(*Strength of control measures*)
c1->0.6,c2->0.6,c3->0.6
};
Print["Case mortality = ",(α/(κyI+α))/.pars];
Print["Case mortality (mutant) = ",((αm)/(κyI+αm))/.pars];

parset=parvec/.pars/.mut->0;
finalfreq[parset]

tNOevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]

fnoevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tNOevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red,Dashed}},AspectRatio->0.75];

fSDashed=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005],Dashed, Blue},AspectRatio->0.75];

mortalitynoevol=yD[tNOevolNC]+yDm[tNOevolNC]/.Flatten[NSolution[parset]];

parset=parvec/.pars/.mut->0.01;
finalfreq[parset]

tevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]
fevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red}},AspectRatio->0.75];

freqERADIC=Evaluate[(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))/.NSolution[parset]]/.t->tevolNC

FIGfreq=Plot[Evaluate[{(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{0,1}},PlotStyle->{{Thickness[0.008],Black}},AspectRatio->0.75];

fS=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005], Blue},AspectRatio->0.75];

fcontrol=Plot[2*c[t]/.pars,{t,0,tfig},Filling->Axis,FillingStyle->LightGray,PlotStyle->None,Frame->True,PlotRange->{{0,tfig},{0,1}},AspectRatio->0.75,Frame->True];

line1=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC-0.02}],Flatten[{tevolNC+5,freqERADIC+0.02}]}]}];
line2=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC+0.02}],Flatten[{tevolNC+5,freqERADIC-0.02}]}]}];

mortalityevol=yD[tevolNC]+yDm[tevolNC]/.Flatten[NSolution[parset]];
(*mortality=Graphics[{Text[Style[StringJoin[ToString[N[Round[100 mortalitynoevol,(10^(-1))] ]],"% (",ToString[N[Round[100 (mortalityevol),(10^(-1))]]],"%)"],12,FontFamily->"Helvetica"],{(t3min+t3max)/2,0.3}]}];*)

DeathRates=Graphics[{Pink,Rectangle[{t3max-tfig/20,0.7},{t3max,0.7+10*mortalitynoevol}],Red,Rectangle[{t3max,0.7},{t3max+tfig/20,0.7+10*mortalityevol}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.2/20,0.705}],Black,Rectangle[{t3max-tfig 1.4/20,0.8},{t3max-tfig 1.2/20,0.805}],Black,Rectangle[{t3max-tfig 1.4/20,0.9},{t3max-tfig 1.2/20,0.905}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.42/20,0.905}]}];

FigBOT1=Show[fcontrol,fS,FIGfreq,fnoevol,fevol,fSDashed,line1,line2, DeathRates,Frame->True,
TicksStyle->Directive[Black,FontFamily->"Helvetica",11],
LabelStyle->Directive[12,FontFamily->"Helvetica",Black],FrameTicks->{{N[{0,0.2,0.4,0.6,0.8,1}],None},{{0,100,200,300,400,500},None}},FrameLabel->{None,None}]

If we increase the pleiotropic effect on the length of the pre-symptomatic phase, however, the mutant allele no longer rises in frequency if κyPm>

tab=Table[Flatten[{κ,finalfreq[parvec/.κyPm->κ/.pars/.mut->0.01];Evaluate[(yEm[tfig]+yAm[tfig]+yPm[tfig]+yIm[tfig])/(yE[tfig]+yA[tfig]+yP[tfig]+yI[tfig]+yEm[tfig]+yAm[tfig]+yPm[tfig]+yIm[tfig])/.Flatten[NSolution[parvec/.κyPm->κ/.pars/.mut->0.01]]]}],{κ,1,2,0.05}];

tab

plot1=ListPlot[tab,Joined->True,PlotStyle->Gray]

Figure 4d: Evolution with negative pleiotropy between transmission β and virulence α, control measures

pars={β->1, βm->1.2, (*TRANSMISSION OF THE MUTANT*)
α->0.005, αm->0.000, (*VIRULENCE OF THE MUTANT*)
κyP->1,κyPm->1, (* 1/TIME IN PRESYMPTOMATIC CLASS FOR THE MUTANT *)

f->0.2, fm->0.2, (* PROPORTION OF ASYMPTOMATIC CASES *)
κyE->0.25, (* 1/TIME IN EXPOSED CLASS *)
κyI->0.2, (* 1/TIME IN SYMPTOMATIC CLASS *)
κyA->0.11, (* 1/TIME IN ASYMPTOMATIC CLASS *)

(*Strength of control measures*)
c1->0.6,c2->0.6,c3->0.6
};
Print["Case mortality = ",(α/(κyI+α))/.pars];
Print["Case mortality (mutant) = ",((αm)/(κyI+αm))/.pars];

parset=parvec/.pars/.mut->0;
finalfreq[parset]

tNOevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]

fnoevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tNOevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red,Dashed}},AspectRatio->0.75];

fSDashed=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005],Dashed, Blue},AspectRatio->0.75];

mortalitynoevol=yD[tNOevolNC]+yDm[tNOevolNC]/.Flatten[NSolution[parset]];

parset=parvec/.pars/.mut->0.01;
finalfreq[parset]

tevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]
fevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red}},AspectRatio->0.75];

freqERADIC=Evaluate[(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))/.NSolution[parset]]/.t->tevolNC

FIGfreq=Plot[Evaluate[{(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{0,1}},PlotStyle->{{Thickness[0.008],Black}},AspectRatio->0.75];

fS=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005], Blue},AspectRatio->0.75];

fcontrol=Plot[2*c[t]/.pars,{t,0,tfig},Filling->Axis,FillingStyle->LightGray,PlotStyle->None,Frame->True,PlotRange->{{0,tfig},{0,1}},AspectRatio->0.75,Frame->True];

line1=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC-0.02}],Flatten[{tevolNC+5,freqERADIC+0.02}]}]}];
line2=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC+0.02}],Flatten[{tevolNC+5,freqERADIC-0.02}]}]}];

mortalityevol=yD[tevolNC]+yDm[tevolNC]/.Flatten[NSolution[parset]];
(*mortality=Graphics[{Text[Style[StringJoin[ToString[N[Round[100 mortalitynoevol,(10^(-1))] ]],"% (",ToString[N[Round[100 (mortalityevol),(10^(-1))]]],"%)"],12,FontFamily->"Helvetica"],{(t3min+t3max)/2,0.3}]}];*)

DeathRates=Graphics[{Pink,Rectangle[{t3max-tfig/20,0.7},{t3max,0.7+10*mortalitynoevol}],Red,Rectangle[{t3max,0.7},{t3max+tfig/20,0.7+10*mortalityevol}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.2/20,0.705}],Black,Rectangle[{t3max-tfig 1.4/20,0.8},{t3max-tfig 1.2/20,0.805}],Black,Rectangle[{t3max-tfig 1.4/20,0.9},{t3max-tfig 1.2/20,0.905}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.42/20,0.905}]}];

FigBOT2=Show[fcontrol,fS,FIGfreq,fnoevol,fevol,fSDashed,line1,line2, DeathRates,Frame->True,
TicksStyle->Directive[Black,FontFamily->"Helvetica",11],
LabelStyle->Directive[12,FontFamily->"Helvetica",Black],FrameTicks->{{N[{0,0.2,0.4,0.6,0.8,1}],None},{{0,100,200,300,400,500},None}},FrameLabel->{None,None}]

Altogether

Show[
GraphicsGrid[{{FigTOP1,FigTOP2},{FigBOT1,FigBOT2}},FrameLabel->label],
Graphics[Text[Style["(a)",16,FontFamily->"Helvetica"],{55,-43}]],
Graphics[Text[Style["(b)",16,FontFamily->"Helvetica"],{443,-43}]],
Graphics[Text[Style["(c)",16,FontFamily->"Helvetica"],{50,-295}]],
Graphics[Text[Style["(d)",16,FontFamily->"Helvetica"],{443,-295}]],
Graphics[Text[Style["Positive trade-off",18,FontFamily->"Helvetica"],{200,0}]],
Graphics[Text[Style["Negative trade-off",18,FontFamily->"Helvetica"],{ 580,0}]],
Graphics[Text[Style["Time (days)",20,FontFamily->"Helvetica"],{400,-520}]],
Graphics[Rotate[Text[Style["Frequency",20,FontFamily->"Helvetica"],{0,-250}],Pi/2]]
]

Export["Fig4.pdf",%];

Figure 3alt - Main change is a longer presymptomatic phase (κyP=0.4)

β adjusted to keep R0 and doubling time of the disease within range observed.

Code

Common settings

Clear["Global`*"]

Off[General::spell1]
Off[General::spell]
Off[NDSolve::nlnum]

tfinal=500;

tfig=500;

label={"Time (days)","Frequency"};
label=None;

parvec={β, βm,α, αm,κyP,κyPm,f,fm,κyE, κyI,κyA, c1,c2,c3,mut};

Clear[finalfreq,NSolution]
finalfreq[vec_]:=finalfreq[vec]=Block[{β=vec[[1]], βm=vec[[2]],α=vec[[3]], αm=vec[[4]],κyP=vec[[5]],κyPm=vec[[6]],f=vec[[7]],fm=vec[[8]],κyE=vec[[9]], κyI=vec[[10]],κyA=vec[[11]],c1=vec[[12]],c2=vec[[13]],c3=vec[[14]],mut=vec[[15]]},

(*REMAINING PARAMETERS*)
βyA=β/10;βyP=β;βyI=β/3;
κyAm=κyA;κyIm=κyI;
βyAm=βm/10;βyPm=βm;βyIm=βm/3;

(* E: Exposed (non-infectious) cases *)
(* A: asymptomatic cases *)
(* P: pre-symptomatic cases *)
(* S: symptomatic cases *)
(* f: proportion of fsymptomatic hosts *)
(* c: amount of distancing applied *)

(*INOCULATION*)
start=(10^(-5));

t=.;
c=.;

c[t_]:=Evaluate[If[t<t1min,0,If[t<t1max,c1,If[t<t2min,0,If[t<t2max,c2,If[t<t3min,0,If[t<t3max,c3,0]]]]]]];
h[t]=(1-c[t])(βyA yA[t]+βyP yP[t]+βyI yI[t]);
hm[t]=(1-c[t])(βyAm yAm[t]+βyPm yPm[t]+βyIm yIm[t]);

sys={S'[t]==-(h[t]+hm[t]) S[t],
yE'[t]==h[t] S[t]-κyE yE[t],
yA'[t]==f κyE yE[t]-κyA yA[t],
yP'[t]==(1-f) κyE yE[t]-κyP yP[t],
yI'[t]==κyP yP[t]-(κyI+α) yI[t],
yR'[t]==κyA yA[t]+κyI yI[t],
yD'[t]==α yI[t]};  (*yD serves as a counter of the cumulative number of deaths*)

sysm={
yEm'[t]==hm[t] S[t]-κyE yEm[t],
yAm'[t]==fm κyE yEm[t]-κyAm yAm[t],
yPm'[t]==(1-fm) κyE yEm[t]-κyPm yPm[t],
yIm'[t]==κyPm yPm[t]-(κyIm+αm) yIm[t],
yRm'[t]==κyAm yAm[t]+κyIm yIm[t],
yDm'[t]==αm yIm[t]};  (*yDm serves as a counter of the cumulative number of deaths*)

(*FIRST BOUT OF CONTROL*)
t1min=75;
t1max=150;
(*SECOND BOUT OF CONTROL*)
t2min=225;
t2max=300;
(*THIRD BOUT OF CONTROL*)
t3min=375;
t3max=450;

init={S[0]==1,yE[0]==start(1-mut),yA[0]==0,yP[0]==0,yI[0]==0,yR[0]==0,yD[0]==0};
initm={yEm[0]==start(mut),yAm[0]==0,yPm[0]==0,yIm[0]==0,yRm[0]==0,yDm[0]==0};

var={S,yE,yA,yP,yI,yR,yD};
varm={yEm,yAm,yPm,yIm,yRm,yDm};

NSolution[vec]=
NDSolve[
Flatten[{sys,sysm,init,initm}],Flatten[{var,varm}],{t,0,tfinal},AccuracyGoal->100];

(*R0*)
Print["R0 of the wild type = ",((βyA f)/(κyA))+(((1-f) (βyP (α+κyI)+βyI κyP))/((α+κyI) κyP))];

(*Rm*)
Print["R0 of the mutant type = ",((βyAm f)/(κyAm))+(((1-f) (βyPm (αm+κyIm)+βyIm κyPm))/((αm+κyIm) κyPm))];

(*CASE MORTALITY*)
Print["Case mortality = ",(1-f)(α/(κyI+α))];
Print["Case mortality (mutant) = ",(1-fm)((αm)/(κyI+αm))];

(*CUMULATIVE MORTALITY*)
Print["Cumulative mortality (total) = ",yD[tfinal]+yDm[tfinal]/.Flatten[NSolution[vec]]];
Print["Cumulative mortality (mutant) = ",yDm[tfinal]/.Flatten[NSolution[vec]]];

]

Figure 3a: Evolution of transmission β, no control measures

pars={β->0.7, βm->0.7 1.2, (*TRANSMISSION OF THE MUTANT*)
α->0.005, αm->0.005, (*VIRULENCE OF THE MUTANT*)
κyP->1/2.5,κyPm->1/2.5, (* 1/TIME IN PRESYMPTOMATIC CLASS FOR THE MUTANT *)

f->0.2, fm->0.2, (* PROPORTION OF ASYMPTOMATIC CASES *)
κyE->0.25, (* 1/TIME IN EXPOSED CLASS *)
κyI->0.2, (* 1/TIME IN SYMPTOMATIC CLASS *)
κyA->0.11, (* 1/TIME IN ASYMPTOMATIC CLASS *)

(*Strength of control measures*)
c1->0,c2->0,c3->0
};

The doublng time of the disease when rare, as calculated from the eigenvalue (“poly”) of the stability matrix:

Sort[Solve[(κE (-S α βP κA+f S α βP κA-S βP κA κI+f S βP κA κI-f S α βA κP+α κA κP-S βI κA κP+f S βI κA κP-f S βA κI κP+κA κI κP)+(-f S α βA κE-S α βP κE+f S α βP κE+α κA κE-S βP κA κE+f S βP κA κE-f S βA κE κI-S βP κE κI+f S βP κE κI+κA κE κI+α κA κP+α κE κP-f S βA κE κP-S βI κE κP+f S βI κE κP+κA κE κP+κA κI κP+κE κI κP) λ+(α κA+α κE-f S βA κE-S βP κE+f S βP κE+κA κE+κA κI+κE κI+α κP+κA κP+κE κP+κI κP) (λ^2)+(α+κA+κE+κI+κP) (λ^3)+(λ^4)/.S->1/.βA->β/10/.βP->β/.βI->β/3/.κA->κyA/.κE->κyE/.κP->κyP/.κI->κyI/.pars)==0,λ]]
Solve[(Exp[λ t]==2)/.Last[%],t]

parset=parvec/.pars/.mut->0;
finalfreq[parset]

tNOevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]

fnoevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tNOevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red,Dashed}},AspectRatio->0.75];

fSDashed=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005],Dashed, Blue},AspectRatio->0.75];

mortalitynoevol=yD[tNOevolNC]+yDm[tNOevolNC]/.Flatten[NSolution[parset]];

parset=parvec/.pars/.mut->0.01;
finalfreq[parset]

tevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]
fevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red}},AspectRatio->0.75];

freqERADIC=Evaluate[(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))/.NSolution[parset]]/.t->tevolNC

FIGfreq=Plot[Evaluate[{(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{0,1}},PlotStyle->{{Thickness[0.008],Black}},AspectRatio->0.75];

fS=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005], Blue},AspectRatio->0.75];

fcontrol=Plot[2*c[t]/.pars,{t,0,tfig},Filling->Axis,FillingStyle->LightGray,PlotStyle->None,Frame->True,PlotRange->{{0,tfig},{0,1}},AspectRatio->0.75,Frame->True];


line1=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC-0.02}],Flatten[{tevolNC+5,freqERADIC+0.02}]}]}];
line2=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC+0.02}],Flatten[{tevolNC+5,freqERADIC-0.02}]}]}];

mortalityevol=yD[tevolNC]+yDm[tevolNC]/.Flatten[NSolution[parset]];
(*mortality=Graphics[{Text[Style[StringJoin[ToString[N[Round[100 mortalitynoevol,(10^(-1))] ]],"% (",ToString[N[Round[100 (mortalityevol),(10^(-1))]]],"%)"],12,FontFamily->"Helvetica"],{(t3min+t3max)/2,0.3}]}];*)

DeathRates=Graphics[{Pink,Rectangle[{t3max-tfig/20,0.7},{t3max,0.7+10*mortalitynoevol}],Red,Rectangle[{t3max,0.7},{t3max+tfig/20,0.7+10*mortalityevol}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.2/20,0.705}],Black,Rectangle[{t3max-tfig 1.4/20,0.8},{t3max-tfig 1.2/20,0.805}],Black,Rectangle[{t3max-tfig 1.4/20,0.9},{t3max-tfig 1.2/20,0.905}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.42/20,0.905}]}];

FigTOP1=Show[fS,FIGfreq,fnoevol,fevol,fSDashed,line1,line2,DeathRates,Frame->True,
TicksStyle->Directive[Black,FontFamily->"Helvetica",11],
LabelStyle->Directive[12,FontFamily->"Helvetica",Black],
FrameTicks->{{N[{0.2,0.4,0.6,0.8}],N[{0.2,0.4,0.6,0.8}]},{{0,100,200,300,400,500},{0,100,200,300,400,500}}},FrameTicksStyle->{{Directive[FontSize->12],Directive[FontOpacity->0,FontSize->0]},{Directive[FontSize->12],Directive[FontOpacity->0,FontSize->0]}}
]

Figure 3b: Evolution of asymptomatic fraction f,  no control measures

pars={β->0.7, βm->0.7 , (*TRANSMISSION OF THE MUTANT*)
α->0.005, αm->0.005, (*VIRULENCE OF THE MUTANT*)
κyP->1/2.5,κyPm->1/2.5, (* 1/TIME IN PRESYMPTOMATIC CLASS FOR THE MUTANT *)

f->0.2, fm->0.1, (* PROPORTION OF ASYMPTOMATIC CASES *)
κyE->0.25, (* 1/TIME IN EXPOSED CLASS *)
κyI->0.2, (* 1/TIME IN SYMPTOMATIC CLASS *)
κyA->0.11, (* 1/TIME IN ASYMPTOMATIC CLASS *)

(*Strength of control measures*)
c1->0,c2->0,c3->0
};

parset=parvec/.pars/.mut->0;
finalfreq[parset]

tNOevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]

fnoevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tNOevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red,Dashed}},AspectRatio->0.75];

fSDashed=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005],Dashed, Blue},AspectRatio->0.75];

mortalitynoevol=yD[tNOevolNC]+yDm[tNOevolNC]/.Flatten[NSolution[parset]];

parset=parvec/.pars/.mut->0.01;
finalfreq[parset]

tevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]
fevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red}},AspectRatio->0.75];

freqERADIC=Evaluate[(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))/.NSolution[parset]]/.t->tevolNC

FIGfreq=Plot[Evaluate[{(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{0,1}},PlotStyle->{{Thickness[0.008],Black}},AspectRatio->0.75];

fS=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005], Blue},AspectRatio->0.75];

fcontrol=Plot[2*c[t]/.pars,{t,0,tfig},Filling->Axis,FillingStyle->LightGray,PlotStyle->None,Frame->True,PlotRange->{{0,tfig},{0,1}},AspectRatio->0.75,Frame->True];


line1=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC-0.02}],Flatten[{tevolNC+5,freqERADIC+0.02}]}]}];
line2=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC+0.02}],Flatten[{tevolNC+5,freqERADIC-0.02}]}]}];


mortalityevol=yD[tevolNC]+yDm[tevolNC]/.Flatten[NSolution[parset]];
(*mortality=Graphics[{Text[Style[StringJoin[ToString[N[Round[100 mortalitynoevol,(10^(-1))] ]],"% (",ToString[N[Round[100 (mortalityevol),(10^(-1))]]],"%)"],12,FontFamily->"Helvetica"],{(t3min+t3max)/2,0.3}]}];*)

DeathRates=Graphics[{Pink,Rectangle[{t3max-tfig/20,0.7},{t3max,0.7+10*mortalitynoevol}],Red,Rectangle[{t3max,0.7},{t3max+tfig/20,0.7+10*mortalityevol}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.2/20,0.705}],Black,Rectangle[{t3max-tfig 1.4/20,0.8},{t3max-tfig 1.2/20,0.805}],Black,Rectangle[{t3max-tfig 1.4/20,0.9},{t3max-tfig 1.2/20,0.905}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.42/20,0.905}]}];

FigTOP2=Show[fS,FIGfreq,fnoevol,fevol,fSDashed,line1,line2,DeathRates,Frame->True,
TicksStyle->Directive[Black,FontFamily->"Helvetica",11],
LabelStyle->Directive[12,FontFamily->"Helvetica",Black],FrameTicks->{{N[{0.2,0.4,0.6,0.8}],N[{0.2,0.4,0.6,0.8}]},{{0,100,200,300,400,500},{0,100,200,300,400,500}}},FrameTicksStyle->{{Directive[FontOpacity->0,FontSize->0],Directive[FontOpacity->0,FontSize->0]},{Directive[FontSize->12],Directive[FontOpacity->0,FontSize->0]}}
]

Figure 3c: Evolution of pre-symptomatic phase κ, no control measures

pars={β->0.7, βm->0.7, (*TRANSMISSION OF THE MUTANT*)
α->0.005, αm->0.005, (*VIRULENCE OF THE MUTANT*)
κyP->1/2.5,κyPm->1/4, (* 1/TIME IN PRESYMPTOMATIC CLASS FOR THE MUTANT *)

f->0.2, fm->0.2, (* PROPORTION OF ASYMPTOMATIC CASES *)
κyE->0.25, (* 1/TIME IN EXPOSED CLASS *)
κyI->0.2, (* 1/TIME IN SYMPTOMATIC CLASS *)
κyA->0.11, (* 1/TIME IN ASYMPTOMATIC CLASS *)

(*Strength of control measures*)
c1->0,c2->0,c3->0
};

parset=parvec/.pars/.mut->0;
finalfreq[parset]

tNOevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]

fnoevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tNOevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red,Dashed}},AspectRatio->0.75];

fSDashed=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005],Dashed, Blue},AspectRatio->0.75];

mortalitynoevol=yD[tNOevolNC]+yDm[tNOevolNC]/.Flatten[NSolution[parset]];

parset=parvec/.pars/.mut->0.01;
finalfreq[parset]

tevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]
fevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red}},AspectRatio->0.75];

freqERADIC=Evaluate[(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))/.NSolution[parset]]/.t->tevolNC

FIGfreq=Plot[Evaluate[{(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{0,1}},PlotStyle->{{Thickness[0.008],Black}},AspectRatio->0.75];

fS=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005], Blue},AspectRatio->0.75];

fcontrol=Plot[2*c[t]/.pars,{t,0,tfig},Filling->Axis,FillingStyle->LightGray,PlotStyle->None,Frame->True,PlotRange->{{0,tfig},{0,1}},AspectRatio->0.75,Frame->True];


line1=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC-0.02}],Flatten[{tevolNC+5,freqERADIC+0.02}]}]}];
line2=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC+0.02}],Flatten[{tevolNC+5,freqERADIC-0.02}]}]}];


mortalityevol=yD[tevolNC]+yDm[tevolNC]/.Flatten[NSolution[parset]];
(*mortality=Graphics[{Text[Style[StringJoin[ToString[N[Round[100 mortalitynoevol,(10^(-1))] ]],"% (",ToString[N[Round[100 (mortalityevol),(10^(-1))]]],"%)"],12,FontFamily->"Helvetica"],{(t3min+t3max)/2,0.3}]}];*)

DeathRates=Graphics[{Pink,Rectangle[{t3max-tfig/20,0.7},{t3max,0.7+10*mortalitynoevol}],Red,Rectangle[{t3max,0.7},{t3max+tfig/20,0.7+10*mortalityevol}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.2/20,0.705}],Black,Rectangle[{t3max-tfig 1.4/20,0.8},{t3max-tfig 1.2/20,0.805}],Black,Rectangle[{t3max-tfig 1.4/20,0.9},{t3max-tfig 1.2/20,0.905}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.42/20,0.905}]}];

FigTOP3=Show[fS,FIGfreq,fnoevol,fevol,fSDashed,line1,line2,DeathRates,Frame->True,
TicksStyle->Directive[Black,FontFamily->"Helvetica",11],
LabelStyle->Directive[12,FontFamily->"Helvetica",Black],FrameTicks->{{N[{0.2,0.4,0.6,0.8}],N[{0.2,0.4,0.6,0.8}]},{{0,100,200,300,400,500},{0,100,200,300,400,500}}},FrameTicksStyle->{{Directive[FontOpacity->0,FontSize->0],Directive[FontOpacity->0,FontSize->0]},{Directive[FontSize->12],Directive[FontOpacity->0,FontSize->0]}}
]

Figure 3d: Evolution of virulence α, no control measures

pars={β->0.7, βm->0.7, (*TRANSMISSION OF THE MUTANT*)
α->0.005, αm->0, (*VIRULENCE OF THE MUTANT*)
κyP->1/2.5,κyPm->1/2.5, (* 1/TIME IN PRESYMPTOMATIC CLASS FOR THE MUTANT *)

f->0.2, fm->0.2, (* PROPORTION OF ASYMPTOMATIC CASES *)
κyE->0.25, (* 1/TIME IN EXPOSED CLASS *)
κyI->0.2, (* 1/TIME IN SYMPTOMATIC CLASS *)
κyA->0.11, (* 1/TIME IN ASYMPTOMATIC CLASS *)

(*Strength of control measures*)
c1->0,c2->0,c3->0
};

parset=parvec/.pars/.mut->0;
finalfreq[parset]

tNOevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]

fnoevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tNOevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red,Dashed}},AspectRatio->0.75];

fSDashed=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005],Dashed, Blue},AspectRatio->0.75];

mortalitynoevol=yD[tNOevolNC]+yDm[tNOevolNC]/.Flatten[NSolution[parset]];

parset=parvec/.pars/.mut->0.01;
finalfreq[parset]

tevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]
fevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red}},AspectRatio->0.75];

freqERADIC=Evaluate[(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))/.NSolution[parset]]/.t->tevolNC

FIGfreq=Plot[Evaluate[{(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{0,1}},PlotStyle->{{Thickness[0.008],Black}},AspectRatio->0.75];

fS=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005], Blue},AspectRatio->0.75];

fcontrol=Plot[2*c[t]/.pars,{t,0,tfig},Filling->Axis,FillingStyle->LightGray,PlotStyle->None,Frame->True,PlotRange->{{0,tfig},{0,1}},AspectRatio->0.75,Frame->True];


line1=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC-0.02}],Flatten[{tevolNC+5,freqERADIC+0.02}]}]}];
line2=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC+0.02}],Flatten[{tevolNC+5,freqERADIC-0.02}]}]}];


mortalityevol=yD[tevolNC]+yDm[tevolNC]/.Flatten[NSolution[parset]];
(*mortality=Graphics[{Text[Style[StringJoin[ToString[N[Round[100 mortalitynoevol,(10^(-1))] ]],"% (",ToString[N[Round[100 (mortalityevol),(10^(-1))]]],"%)"],12,FontFamily->"Helvetica"],{(t3min+t3max)/2,0.3}]}];*)

DeathRates=Graphics[{Pink,Rectangle[{t3max-tfig/20,0.7},{t3max,0.7+10*mortalitynoevol}],Red,Rectangle[{t3max,0.7},{t3max+tfig/20,0.7+10*mortalityevol}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.2/20,0.705}],Black,Rectangle[{t3max-tfig 1.4/20,0.8},{t3max-tfig 1.2/20,0.805}],Black,Rectangle[{t3max-tfig 1.4/20,0.9},{t3max-tfig 1.2/20,0.905}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.42/20,0.905}]}];

FigTOP4=Show[fS,FIGfreq,fnoevol,fevol,fSDashed,line1,line2,DeathRates,Frame->True,
TicksStyle->Directive[Black,FontFamily->"Helvetica",11],
LabelStyle->Directive[12,FontFamily->"Helvetica",Black],FrameTicks->{{N[{0.2,0.4,0.6,0.8}],N[{0.2,0.4,0.6,0.8}]},{{0,100,200,300,400,500},{0,100,200,300,400,500}}},FrameTicksStyle->{{Directive[FontOpacity->0,FontSize->0],Directive[FontOpacity->0,FontSize->0]},{Directive[FontSize->12],Directive[FontOpacity->0,FontSize->0]}}
]

Figure 3e: Evolution of transmission β, control measures

pars={β->0.7, βm->0.7 1.2, (*TRANSMISSION OF THE MUTANT*)
α->0.005, αm->0.005, (*VIRULENCE OF THE MUTANT*)
κyP->1/2.5,κyPm->1/2.5, (* 1/TIME IN PRESYMPTOMATIC CLASS FOR THE MUTANT *)

f->0.2, fm->0.2, (* PROPORTION OF ASYMPTOMATIC CASES *)
κyE->0.25, (* 1/TIME IN EXPOSED CLASS *)
κyI->0.2, (* 1/TIME IN SYMPTOMATIC CLASS *)
κyA->0.11, (* 1/TIME IN ASYMPTOMATIC CLASS *)

(*Strength of control measures*)
c1->0.6,c2->0.6,c3->0.6
};

parset=parvec/.pars/.mut->0;
finalfreq[parset]

tNOevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]

fnoevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tNOevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red,Dashed}},AspectRatio->0.75];

fSDashed=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005],Dashed, Blue},AspectRatio->0.75];

mortalitynoevol=yD[tNOevolNC]+yDm[tNOevolNC]/.Flatten[NSolution[parset]];

parset=parvec/.pars/.mut->0.01;
finalfreq[parset]

tevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]
fevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red}},AspectRatio->0.75];

freqERADIC=Evaluate[(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))/.NSolution[parset]]/.t->tevolNC

FIGfreq=Plot[Evaluate[{(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{0,1}},PlotStyle->{{Thickness[0.008],Black}},AspectRatio->0.75];

fS=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005], Blue},AspectRatio->0.75];

fcontrol=Plot[2*c[t]/.pars,{t,0,tfig},Filling->Axis,FillingStyle->LightGray,PlotStyle->None,Frame->True,PlotRange->{{0,tfig},{0,1}},AspectRatio->0.75,Frame->True];

line1=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC-0.02}],Flatten[{tevolNC+5,freqERADIC+0.02}]}]}];
line2=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC+0.02}],Flatten[{tevolNC+5,freqERADIC-0.02}]}]}];

mortalityevol=yD[tevolNC]+yDm[tevolNC]/.Flatten[NSolution[parset]];
(*mortality=Graphics[{Text[Style[StringJoin[ToString[N[Round[100 mortalitynoevol,(10^(-1))] ]],"% (",ToString[N[Round[100 (mortalityevol),(10^(-1))]]],"%)"],12,FontFamily->"Helvetica"],{(t3min+t3max)/2,0.3}]}];*)

DeathRates=Graphics[{Pink,Rectangle[{t3max-tfig/20,0.7},{t3max,0.7+10*mortalitynoevol}],Red,Rectangle[{t3max,0.7},{t3max+tfig/20,0.7+10*mortalityevol}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.2/20,0.705}],Black,Rectangle[{t3max-tfig 1.4/20,0.8},{t3max-tfig 1.2/20,0.805}],Black,Rectangle[{t3max-tfig 1.4/20,0.9},{t3max-tfig 1.2/20,0.905}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.42/20,0.905}]}];

FigBOT1=Show[fcontrol,fS,fSDashed, FIGfreq,fnoevol,fevol,line1,line2,DeathRates,Frame->True,
TicksStyle->Directive[Black,FontFamily->"Helvetica",11],
LabelStyle->Directive[12,FontFamily->"Helvetica",Black],FrameTicks->{{N[{0.2,0.4,0.6,0.8}],N[{0.2,0.4,0.6,0.8}]},{{0,100,200,300,400,500},{0,100,200,300,400,500}}},FrameTicksStyle->{{Directive[FontSize->12],Directive[FontOpacity->0,FontSize->0]},{Directive[FontSize->12],Directive[FontOpacity->0,FontSize->0]}}
]

Figure 3f: Evolution of asymptomatic fraction f, control measures

pars={β->0.7, βm->0.7, (*TRANSMISSION OF THE MUTANT*)
α->0.005, αm->0.005, (*VIRULENCE OF THE MUTANT*)
κyP->1/2.5,κyPm->1/2.5, (* 1/TIME IN PRESYMPTOMATIC CLASS FOR THE MUTANT *)

f->0.2, fm->0.1, (* PROPORTION OF ASYMPTOMATIC CASES *)
κyE->0.25, (* 1/TIME IN EXPOSED CLASS *)
κyI->0.2, (* 1/TIME IN SYMPTOMATIC CLASS *)
κyA->0.11, (* 1/TIME IN ASYMPTOMATIC CLASS *)

(*Strength of control measures*)
c1->0.6,c2->0.6,c3->0.6
};

parset=parvec/.pars/.mut->0;
finalfreq[parset]

tNOevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]

fnoevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tNOevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red,Dashed}},AspectRatio->0.75];

fSDashed=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005],Dashed, Blue},AspectRatio->0.75];

mortalitynoevol=yD[tNOevolNC]+yDm[tNOevolNC]/.Flatten[NSolution[parset]];

parset=parvec/.pars/.mut->0.01;
finalfreq[parset]

tevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]
fevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red}},AspectRatio->0.75];

freqERADIC=Evaluate[(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))/.NSolution[parset]]/.t->tevolNC

FIGfreq=Plot[Evaluate[{(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{0,1}},PlotStyle->{{Thickness[0.008],Black}},AspectRatio->0.75];

fS=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005], Blue},AspectRatio->0.75];

fcontrol=Plot[2*c[t]/.pars,{t,0,tfig},Filling->Axis,FillingStyle->LightGray,PlotStyle->None,Frame->True,PlotRange->{{0,tfig},{0,1}},AspectRatio->0.75,Frame->True];

line1=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC-0.02}],Flatten[{tevolNC+5,freqERADIC+0.02}]}]}];
line2=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC+0.02}],Flatten[{tevolNC+5,freqERADIC-0.02}]}]}];

mortalityevol=yD[tevolNC]+yDm[tevolNC]/.Flatten[NSolution[parset]];
(*mortality=Graphics[{Text[Style[StringJoin[ToString[N[Round[100 mortalitynoevol,(10^(-1))] ]],"% (",ToString[N[Round[100 (mortalityevol),(10^(-1))]]],"%)"],12,FontFamily->"Helvetica"],{(t3min+t3max)/2,0.3}]}];*)

DeathRates=Graphics[{Pink,Rectangle[{t3max-tfig/20,0.7},{t3max,0.7+10*mortalitynoevol}],Red,Rectangle[{t3max,0.7},{t3max+tfig/20,0.7+10*mortalityevol}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.2/20,0.705}],Black,Rectangle[{t3max-tfig 1.4/20,0.8},{t3max-tfig 1.2/20,0.805}],Black,Rectangle[{t3max-tfig 1.4/20,0.9},{t3max-tfig 1.2/20,0.905}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.42/20,0.905}]}];

FigBOT2=Show[fcontrol,fS,FIGfreq,fnoevol,fevol,fSDashed,line1,line2,DeathRates,Frame->True,
TicksStyle->Directive[Black,FontFamily->"Helvetica",11],
LabelStyle->Directive[12,FontFamily->"Helvetica",Black],FrameTicks->{{N[{0.2,0.4,0.6,0.8}],N[{0.2,0.4,0.6,0.8}]},{{0,100,200,300,400,500},{0,100,200,300,400,500}}},FrameTicksStyle->{{Directive[FontOpacity->0,FontSize->0],Directive[FontOpacity->0,FontSize->0]},{Directive[FontSize->12],Directive[FontOpacity->0,FontSize->0]}}
]

Figure 3g: Evolution of pre-symptomatic phase κ, control measures

pars={β->0.7, βm->0.7, (*TRANSMISSION OF THE MUTANT*)
α->0.005, αm->0.005, (*VIRULENCE OF THE MUTANT*)
κyP->1/2.5,κyPm->1/4, (* 1/TIME IN PRESYMPTOMATIC CLASS FOR THE MUTANT *)

f->0.2, fm->0.2, (* PROPORTION OF ASYMPTOMATIC CASES *)
κyE->0.25, (* 1/TIME IN EXPOSED CLASS *)
κyI->0.2, (* 1/TIME IN SYMPTOMATIC CLASS *)
κyA->0.11, (* 1/TIME IN ASYMPTOMATIC CLASS *)

(*Strength of control measures*)
c1->0.6,c2->0.6,c3->0.6
};

parset=parvec/.pars/.mut->0;
finalfreq[parset]

tNOevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]

fnoevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tNOevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red,Dashed}},AspectRatio->0.75];

fSDashed=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005],Dashed, Blue},AspectRatio->0.75];

mortalitynoevol=yD[tNOevolNC]+yDm[tNOevolNC]/.Flatten[NSolution[parset]];

parset=parvec/.pars/.mut->0.01;
finalfreq[parset]

tevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]
fevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red}},AspectRatio->0.75];

freqERADIC=Evaluate[(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))/.NSolution[parset]]/.t->tevolNC

FIGfreq=Plot[Evaluate[{(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{0,1}},PlotStyle->{{Thickness[0.008],Black}},AspectRatio->0.75];

fS=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005], Blue},AspectRatio->0.75];

fcontrol=Plot[2*c[t]/.pars,{t,0,tfig},Filling->Axis,FillingStyle->LightGray,PlotStyle->None,Frame->True,PlotRange->{{0,tfig},{0,1}},AspectRatio->0.75,Frame->True];

line1=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC-0.02}],Flatten[{tevolNC+5,freqERADIC+0.02}]}]}];
line2=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC+0.02}],Flatten[{tevolNC+5,freqERADIC-0.02}]}]}];

mortalityevol=yD[tevolNC]+yDm[tevolNC]/.Flatten[NSolution[parset]];
(*mortality=Graphics[{Text[Style[StringJoin[ToString[N[Round[100 mortalitynoevol,(10^(-1))] ]],"% (",ToString[N[Round[100 (mortalityevol),(10^(-1))]]],"%)"],12,FontFamily->"Helvetica"],{(t3min+t3max)/2,0.3}]}];*)

DeathRates=Graphics[{Pink,Rectangle[{t3max-tfig/20,0.7},{t3max,0.7+10*mortalitynoevol}],Red,Rectangle[{t3max,0.7},{t3max+tfig/20,0.7+10*mortalityevol}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.2/20,0.705}],Black,Rectangle[{t3max-tfig 1.4/20,0.8},{t3max-tfig 1.2/20,0.805}],Black,Rectangle[{t3max-tfig 1.4/20,0.9},{t3max-tfig 1.2/20,0.905}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.42/20,0.905}]}];

FigBOT3=Show[fcontrol,fS,FIGfreq,fnoevol,fevol,fSDashed,line1,line2,DeathRates,Frame->True,
TicksStyle->Directive[Black,FontFamily->"Helvetica",11],
LabelStyle->Directive[12,FontFamily->"Helvetica",Black],FrameTicks->{{N[{0.2,0.4,0.6,0.8}],N[{0.2,0.4,0.6,0.8}]},{{0,100,200,300,400,500},{0,100,200,300,400,500}}},FrameTicksStyle->{{Directive[FontOpacity->0,FontSize->0],Directive[FontOpacity->0,FontSize->0]},{Directive[FontSize->12],Directive[FontOpacity->0,FontSize->0]}}
]

Figure 3h: Evolution of virulence α, control measures

pars={β->0.7, βm->0.7, (*TRANSMISSION OF THE MUTANT*)
α->0.005, αm->0, (*VIRULENCE OF THE MUTANT*)
κyP->1/2.5,κyPm->1/2.5, (* 1/TIME IN PRESYMPTOMATIC CLASS FOR THE MUTANT *)

f->0.2, fm->0.2, (* PROPORTION OF ASYMPTOMATIC CASES *)
κyE->0.25, (* 1/TIME IN EXPOSED CLASS *)
κyI->0.2, (* 1/TIME IN SYMPTOMATIC CLASS *)
κyA->0.11, (* 1/TIME IN ASYMPTOMATIC CLASS *)

(*Strength of control measures*)
c1->0.6,c2->0.6,c3->0.6
};

parset=parvec/.pars/.mut->0;
finalfreq[parset]

tNOevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]

fnoevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tNOevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red,Dashed}},AspectRatio->0.75];

fSDashed=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005],Dashed, Blue},AspectRatio->0.75];

mortalitynoevol=yD[tNOevolNC]+yDm[tNOevolNC]/.Flatten[NSolution[parset]];

parset=parvec/.pars/.mut->0.01;
finalfreq[parset]

tevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]
fevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red}},AspectRatio->0.75];

freqERADIC=Evaluate[(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))/.NSolution[parset]]/.t->tevolNC

FIGfreq=Plot[Evaluate[{(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{0,1}},PlotStyle->{{Thickness[0.008],Black}},AspectRatio->0.75];

fS=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005], Blue},AspectRatio->0.75];

fcontrol=Plot[2*c[t]/.pars,{t,0,tfig},Filling->Axis,FillingStyle->LightGray,PlotStyle->None,Frame->True,PlotRange->{{0,tfig},{0,1}},AspectRatio->0.75,Frame->True];

line1=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC-0.02}],Flatten[{tevolNC+5,freqERADIC+0.02}]}]}];
line2=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC+0.02}],Flatten[{tevolNC+5,freqERADIC-0.02}]}]}];

mortalityevol=yD[tevolNC]+yDm[tevolNC]/.Flatten[NSolution[parset]];
(*mortality=Graphics[{Text[Style[StringJoin[ToString[N[Round[100 mortalitynoevol,(10^(-1))] ]],"% (",ToString[N[Round[100 (mortalityevol),(10^(-1))]]],"%)"],12,FontFamily->"Helvetica"],{(t3min+t3max)/2,0.3}]}];*)

DeathRates=Graphics[{Pink,Rectangle[{t3max-tfig/20,0.7},{t3max,0.7+10*mortalitynoevol}],Red,Rectangle[{t3max,0.7},{t3max+tfig/20,0.7+10*mortalityevol}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.2/20,0.705}],Black,Rectangle[{t3max-tfig 1.4/20,0.8},{t3max-tfig 1.2/20,0.805}],Black,Rectangle[{t3max-tfig 1.4/20,0.9},{t3max-tfig 1.2/20,0.905}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.42/20,0.905}]}];

FigBOT4=Show[fcontrol,fS,FIGfreq,fnoevol,fevol,fSDashed,line1,line2,DeathRates,Frame->True,
TicksStyle->Directive[Black,FontFamily->"Helvetica",11],
LabelStyle->Directive[12,FontFamily->"Helvetica",Black],FrameTicks->{{N[{0.2,0.4,0.6,0.8}],N[{0.2,0.4,0.6,0.8}]},{{0,100,200,300,400,500},{0,100,200,300,400,500}}},FrameTicksStyle->{{Directive[FontOpacity->0,FontSize->0],Directive[FontOpacity->0,FontSize->0]},{Directive[FontSize->12],Directive[FontOpacity->0,FontSize->0]}}
]

Altogether

Show[
GraphicsGrid[{{FigTOP1,FigTOP2,FigTOP3,FigTOP4},{FigBOT1,FigBOT2,FigBOT3,FigBOT4}},FrameLabel->label,Spacings->{-10,10}],
Graphics[Text[Style["(a)",16,FontFamily->"Helvetica"],{50,-33}]],
Graphics[Text[Style["(b)",16,FontFamily->"Helvetica"],{384,-33}]],
Graphics[Text[Style["(c)",16,FontFamily->"Helvetica"],{734,-33}]],
Graphics[Text[Style["(d)",16,FontFamily->"Helvetica"],{1084,-33}]],
Graphics[Text[Style["(e)",16,FontFamily->"Helvetica"],{50,-315}]],
Graphics[Text[Style["(f)",16,FontFamily->"Helvetica"],{384,-315}]],
Graphics[Text[Style["(g)",16,FontFamily->"Helvetica"],{734,-315}]],
Graphics[Text[Style["(h)",16,FontFamily->"Helvetica"],{1084,-315}]],
Graphics[Text[Style["Transmission",16,FontFamily->"Helvetica"],{180,10}]],
Graphics[Text[Style["Asymptomatic fraction",16,FontFamily->"Helvetica"],{ 520,10}]],
Graphics[Text[Style["Pre-symptomatic phase",16,FontFamily->"Helvetica"],{870,10}]],
Graphics[Text[Style["Virulence",16,FontFamily->"Helvetica"],{1220,10}]],
Graphics[Text[Style["Time (days)",16,FontFamily->"Helvetica"],{700,-580}]],
Graphics[Rotate[Text[Style["Frequency",16,FontFamily->"Helvetica"],{-20,-270}],Pi/2]]
]

Figure 3alt - Main change is more asymptomatics (f=0.4)

β adjusted to keep R0 and doubling time of the disease within range observed.

Code

Common settings

Clear["Global`*"]

Off[General::spell1]
Off[General::spell]
Off[NDSolve::nlnum]

tfinal=500;

tfig=500;

label={"Time (days)","Frequency"};
label=None;

parvec={β, βm,α, αm,κyP,κyPm,f,fm,κyE, κyI,κyA, c1,c2,c3,mut};

Clear[finalfreq,NSolution]
finalfreq[vec_]:=finalfreq[vec]=Block[{β=vec[[1]], βm=vec[[2]],α=vec[[3]], αm=vec[[4]],κyP=vec[[5]],κyPm=vec[[6]],f=vec[[7]],fm=vec[[8]],κyE=vec[[9]], κyI=vec[[10]],κyA=vec[[11]],c1=vec[[12]],c2=vec[[13]],c3=vec[[14]],mut=vec[[15]]},

(*REMAINING PARAMETERS*)
βyA=β/10;βyP=β;βyI=β/3;
κyAm=κyA;κyIm=κyI;
βyAm=βm/10;βyPm=βm;βyIm=βm/3;

(* E: Exposed (non-infectious) cases *)
(* A: asymptomatic cases *)
(* P: pre-symptomatic cases *)
(* S: symptomatic cases *)
(* f: proportion of fsymptomatic hosts *)
(* c: amount of distancing applied *)

(*INOCULATION*)
start=(10^(-5));

t=.;
c=.;

c[t_]:=Evaluate[If[t<t1min,0,If[t<t1max,c1,If[t<t2min,0,If[t<t2max,c2,If[t<t3min,0,If[t<t3max,c3,0]]]]]]];
h[t]=(1-c[t])(βyA yA[t]+βyP yP[t]+βyI yI[t]);
hm[t]=(1-c[t])(βyAm yAm[t]+βyPm yPm[t]+βyIm yIm[t]);

sys={S'[t]==-(h[t]+hm[t]) S[t],
yE'[t]==h[t] S[t]-κyE yE[t],
yA'[t]==f κyE yE[t]-κyA yA[t],
yP'[t]==(1-f) κyE yE[t]-κyP yP[t],
yI'[t]==κyP yP[t]-(κyI+α) yI[t],
yR'[t]==κyA yA[t]+κyI yI[t],
yD'[t]==α yI[t]};  (*yD serves as a counter of the cumulative number of deaths*)

sysm={
yEm'[t]==hm[t] S[t]-κyE yEm[t],
yAm'[t]==fm κyE yEm[t]-κyAm yAm[t],
yPm'[t]==(1-fm) κyE yEm[t]-κyPm yPm[t],
yIm'[t]==κyPm yPm[t]-(κyIm+αm) yIm[t],
yRm'[t]==κyAm yAm[t]+κyIm yIm[t],
yDm'[t]==αm yIm[t]};  (*yDm serves as a counter of the cumulative number of deaths*)

(*FIRST BOUT OF CONTROL*)
t1min=75;
t1max=150;
(*SECOND BOUT OF CONTROL*)
t2min=225;
t2max=300;
(*THIRD BOUT OF CONTROL*)
t3min=375;
t3max=450;

init={S[0]==1,yE[0]==start(1-mut),yA[0]==0,yP[0]==0,yI[0]==0,yR[0]==0,yD[0]==0};
initm={yEm[0]==start(mut),yAm[0]==0,yPm[0]==0,yIm[0]==0,yRm[0]==0,yDm[0]==0};

var={S,yE,yA,yP,yI,yR,yD};
varm={yEm,yAm,yPm,yIm,yRm,yDm};

NSolution[vec]=
NDSolve[
Flatten[{sys,sysm,init,initm}],Flatten[{var,varm}],{t,0,tfinal},AccuracyGoal->100];

(*R0*)
Print["R0 of the wild type = ",((βyA f)/(κyA))+(((1-f) (βyP (α+κyI)+βyI κyP))/((α+κyI) κyP))];

(*Rm*)
Print["R0 of the mutant type = ",((βyAm f)/(κyAm))+(((1-f) (βyPm (αm+κyIm)+βyIm κyPm))/((αm+κyIm) κyPm))];

(*CASE MORTALITY*)
Print["Case mortality = ",(1-f)(α/(κyI+α))];
Print["Case mortality (mutant) = ",(1-fm)((αm)/(κyI+αm))];

(*CUMULATIVE MORTALITY*)
Print["Cumulative mortality (total) = ",yD[tfinal]+yDm[tfinal]/.Flatten[NSolution[vec]]];
Print["Cumulative mortality (mutant) = ",yDm[tfinal]/.Flatten[NSolution[vec]]];

]

Figure 3a: Evolution of transmission β, no control measures

pars={β->1.3, βm->1.2 1.3, (*TRANSMISSION OF THE MUTANT*)
α->0.005, αm->0.005, (*VIRULENCE OF THE MUTANT*)
κyP->1,κyPm->1, (* 1/TIME IN PRESYMPTOMATIC CLASS FOR THE MUTANT *)

f->0.4, fm->0.4, (* PROPORTION OF ASYMPTOMATIC CASES *)
κyE->0.25, (* 1/TIME IN EXPOSED CLASS *)
κyI->0.2, (* 1/TIME IN SYMPTOMATIC CLASS *)
κyA->0.11, (* 1/TIME IN ASYMPTOMATIC CLASS *)

(*Strength of control measures*)
c1->0,c2->0,c3->0
};

The doublng time of the disease when rare, as calculated from the eigenvalue (“poly”) of the stability matrix:

Sort[Solve[(κE (-S α βP κA+f S α βP κA-S βP κA κI+f S βP κA κI-f S α βA κP+α κA κP-S βI κA κP+f S βI κA κP-f S βA κI κP+κA κI κP)+(-f S α βA κE-S α βP κE+f S α βP κE+α κA κE-S βP κA κE+f S βP κA κE-f S βA κE κI-S βP κE κI+f S βP κE κI+κA κE κI+α κA κP+α κE κP-f S βA κE κP-S βI κE κP+f S βI κE κP+κA κE κP+κA κI κP+κE κI κP) λ+(α κA+α κE-f S βA κE-S βP κE+f S βP κE+κA κE+κA κI+κE κI+α κP+κA κP+κE κP+κI κP) (λ^2)+(α+κA+κE+κI+κP) (λ^3)+(λ^4)/.S->1/.βA->β/10/.βP->β/.βI->β/3/.κA->κyA/.κE->κyE/.κP->κyP/.κI->κyI/.pars)==0,λ]]
Solve[(Exp[λ t]==2)/.Last[%],t]

parset=parvec/.pars/.mut->0;
finalfreq[parset]

tNOevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]

fnoevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tNOevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red,Dashed}},AspectRatio->0.75];

fSDashed=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005],Dashed, Blue},AspectRatio->0.75];

mortalitynoevol=yD[tNOevolNC]+yDm[tNOevolNC]/.Flatten[NSolution[parset]];

parset=parvec/.pars/.mut->0.01;
finalfreq[parset]

tevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]
fevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red}},AspectRatio->0.75];

freqERADIC=Evaluate[(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))/.NSolution[parset]]/.t->tevolNC

FIGfreq=Plot[Evaluate[{(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{0,1}},PlotStyle->{{Thickness[0.008],Black}},AspectRatio->0.75];

fS=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005], Blue},AspectRatio->0.75];

fcontrol=Plot[2*c[t]/.pars,{t,0,tfig},Filling->Axis,FillingStyle->LightGray,PlotStyle->None,Frame->True,PlotRange->{{0,tfig},{0,1}},AspectRatio->0.75,Frame->True];


line1=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC-0.02}],Flatten[{tevolNC+5,freqERADIC+0.02}]}]}];
line2=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC+0.02}],Flatten[{tevolNC+5,freqERADIC-0.02}]}]}];

mortalityevol=yD[tevolNC]+yDm[tevolNC]/.Flatten[NSolution[parset]];
(*mortality=Graphics[{Text[Style[StringJoin[ToString[N[Round[100 mortalitynoevol,(10^(-1))] ]],"% (",ToString[N[Round[100 (mortalityevol),(10^(-1))]]],"%)"],12,FontFamily->"Helvetica"],{(t3min+t3max)/2,0.3}]}];*)

DeathRates=Graphics[{Pink,Rectangle[{t3max-tfig/20,0.7},{t3max,0.7+10*mortalitynoevol}],Red,Rectangle[{t3max,0.7},{t3max+tfig/20,0.7+10*mortalityevol}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.2/20,0.705}],Black,Rectangle[{t3max-tfig 1.4/20,0.8},{t3max-tfig 1.2/20,0.805}],Black,Rectangle[{t3max-tfig 1.4/20,0.9},{t3max-tfig 1.2/20,0.905}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.42/20,0.905}]}];

FigTOP1=Show[fS,FIGfreq,fnoevol,fevol,fSDashed,line1,line2,DeathRates,Frame->True,
TicksStyle->Directive[Black,FontFamily->"Helvetica",11],
LabelStyle->Directive[12,FontFamily->"Helvetica",Black],
FrameTicks->{{N[{0.2,0.4,0.6,0.8}],N[{0.2,0.4,0.6,0.8}]},{{0,100,200,300,400,500},{0,100,200,300,400,500}}},FrameTicksStyle->{{Directive[FontSize->12],Directive[FontOpacity->0,FontSize->0]},{Directive[FontSize->12],Directive[FontOpacity->0,FontSize->0]}}
]

Figure 3b: Evolution of asymptomatic fraction f,  no control measures

pars={β->1.3, βm->1.3, (*TRANSMISSION OF THE MUTANT*)
α->0.005, αm->0.005, (*VIRULENCE OF THE MUTANT*)
κyP->1,κyPm->1, (* 1/TIME IN PRESYMPTOMATIC CLASS FOR THE MUTANT *)

f->0.4, fm->0.3, (* PROPORTION OF ASYMPTOMATIC CASES *)
κyE->0.25, (* 1/TIME IN EXPOSED CLASS *)
κyI->0.2, (* 1/TIME IN SYMPTOMATIC CLASS *)
κyA->0.11, (* 1/TIME IN ASYMPTOMATIC CLASS *)

(*Strength of control measures*)
c1->0,c2->0,c3->0
};

parset=parvec/.pars/.mut->0;
finalfreq[parset]

tNOevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]

fnoevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tNOevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red,Dashed}},AspectRatio->0.75];

fSDashed=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005],Dashed, Blue},AspectRatio->0.75];

mortalitynoevol=yD[tNOevolNC]+yDm[tNOevolNC]/.Flatten[NSolution[parset]];

parset=parvec/.pars/.mut->0.01;
finalfreq[parset]

tevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]
fevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red}},AspectRatio->0.75];

freqERADIC=Evaluate[(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))/.NSolution[parset]]/.t->tevolNC

FIGfreq=Plot[Evaluate[{(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{0,1}},PlotStyle->{{Thickness[0.008],Black}},AspectRatio->0.75];

fS=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005], Blue},AspectRatio->0.75];

fcontrol=Plot[2*c[t]/.pars,{t,0,tfig},Filling->Axis,FillingStyle->LightGray,PlotStyle->None,Frame->True,PlotRange->{{0,tfig},{0,1}},AspectRatio->0.75,Frame->True];


line1=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC-0.02}],Flatten[{tevolNC+5,freqERADIC+0.02}]}]}];
line2=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC+0.02}],Flatten[{tevolNC+5,freqERADIC-0.02}]}]}];


mortalityevol=yD[tevolNC]+yDm[tevolNC]/.Flatten[NSolution[parset]];
(*mortality=Graphics[{Text[Style[StringJoin[ToString[N[Round[100 mortalitynoevol,(10^(-1))] ]],"% (",ToString[N[Round[100 (mortalityevol),(10^(-1))]]],"%)"],12,FontFamily->"Helvetica"],{(t3min+t3max)/2,0.3}]}];*)

DeathRates=Graphics[{Pink,Rectangle[{t3max-tfig/20,0.7},{t3max,0.7+10*mortalitynoevol}],Red,Rectangle[{t3max,0.7},{t3max+tfig/20,0.7+10*mortalityevol}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.2/20,0.705}],Black,Rectangle[{t3max-tfig 1.4/20,0.8},{t3max-tfig 1.2/20,0.805}],Black,Rectangle[{t3max-tfig 1.4/20,0.9},{t3max-tfig 1.2/20,0.905}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.42/20,0.905}]}];

FigTOP2=Show[fS,FIGfreq,fnoevol,fevol,fSDashed,line1,line2,DeathRates,Frame->True,
TicksStyle->Directive[Black,FontFamily->"Helvetica",11],
LabelStyle->Directive[12,FontFamily->"Helvetica",Black],FrameTicks->{{N[{0.2,0.4,0.6,0.8}],N[{0.2,0.4,0.6,0.8}]},{{0,100,200,300,400,500},{0,100,200,300,400,500}}},FrameTicksStyle->{{Directive[FontOpacity->0,FontSize->0],Directive[FontOpacity->0,FontSize->0]},{Directive[FontSize->12],Directive[FontOpacity->0,FontSize->0]}}
]

Figure 3c: Evolution of pre-symptomatic phase κ, no control measures

pars={β->1.3, βm->1.3, (*TRANSMISSION OF THE MUTANT*)
α->0.005, αm->0.005, (*VIRULENCE OF THE MUTANT*)
κyP->1,κyPm->2/3, (* 1/TIME IN PRESYMPTOMATIC CLASS FOR THE MUTANT *)

f->0.4, fm->0.4, (* PROPORTION OF ASYMPTOMATIC CASES *)
κyE->0.25, (* 1/TIME IN EXPOSED CLASS *)
κyI->0.2, (* 1/TIME IN SYMPTOMATIC CLASS *)
κyA->0.11, (* 1/TIME IN ASYMPTOMATIC CLASS *)

(*Strength of control measures*)
c1->0,c2->0,c3->0
};

parset=parvec/.pars/.mut->0;
finalfreq[parset]

tNOevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]

fnoevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tNOevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red,Dashed}},AspectRatio->0.75];

fSDashed=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005],Dashed, Blue},AspectRatio->0.75];

mortalitynoevol=yD[tNOevolNC]+yDm[tNOevolNC]/.Flatten[NSolution[parset]];

parset=parvec/.pars/.mut->0.01;
finalfreq[parset]

tevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]
fevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red}},AspectRatio->0.75];

freqERADIC=Evaluate[(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))/.NSolution[parset]]/.t->tevolNC

FIGfreq=Plot[Evaluate[{(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{0,1}},PlotStyle->{{Thickness[0.008],Black}},AspectRatio->0.75];

fS=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005], Blue},AspectRatio->0.75];

fcontrol=Plot[2*c[t]/.pars,{t,0,tfig},Filling->Axis,FillingStyle->LightGray,PlotStyle->None,Frame->True,PlotRange->{{0,tfig},{0,1}},AspectRatio->0.75,Frame->True];


line1=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC-0.02}],Flatten[{tevolNC+5,freqERADIC+0.02}]}]}];
line2=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC+0.02}],Flatten[{tevolNC+5,freqERADIC-0.02}]}]}];


mortalityevol=yD[tevolNC]+yDm[tevolNC]/.Flatten[NSolution[parset]];
(*mortality=Graphics[{Text[Style[StringJoin[ToString[N[Round[100 mortalitynoevol,(10^(-1))] ]],"% (",ToString[N[Round[100 (mortalityevol),(10^(-1))]]],"%)"],12,FontFamily->"Helvetica"],{(t3min+t3max)/2,0.3}]}];*)

DeathRates=Graphics[{Pink,Rectangle[{t3max-tfig/20,0.7},{t3max,0.7+10*mortalitynoevol}],Red,Rectangle[{t3max,0.7},{t3max+tfig/20,0.7+10*mortalityevol}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.2/20,0.705}],Black,Rectangle[{t3max-tfig 1.4/20,0.8},{t3max-tfig 1.2/20,0.805}],Black,Rectangle[{t3max-tfig 1.4/20,0.9},{t3max-tfig 1.2/20,0.905}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.42/20,0.905}]}];

FigTOP3=Show[fS,FIGfreq,fnoevol,fevol,fSDashed,line1,line2,DeathRates,Frame->True,
TicksStyle->Directive[Black,FontFamily->"Helvetica",11],
LabelStyle->Directive[12,FontFamily->"Helvetica",Black],FrameTicks->{{N[{0.2,0.4,0.6,0.8}],N[{0.2,0.4,0.6,0.8}]},{{0,100,200,300,400,500},{0,100,200,300,400,500}}},FrameTicksStyle->{{Directive[FontOpacity->0,FontSize->0],Directive[FontOpacity->0,FontSize->0]},{Directive[FontSize->12],Directive[FontOpacity->0,FontSize->0]}}
]

Figure 3d: Evolution of virulence α, no control measures

pars={β->1.3, βm->1.3, (*TRANSMISSION OF THE MUTANT*)
α->0.005, αm->0.0, (*VIRULENCE OF THE MUTANT*)
κyP->1,κyPm->1, (* 1/TIME IN PRESYMPTOMATIC CLASS FOR THE MUTANT *)

f->0.4, fm->0.4, (* PROPORTION OF ASYMPTOMATIC CASES *)
κyE->0.25, (* 1/TIME IN EXPOSED CLASS *)
κyI->0.2, (* 1/TIME IN SYMPTOMATIC CLASS *)
κyA->0.11, (* 1/TIME IN ASYMPTOMATIC CLASS *)

(*Strength of control measures*)
c1->0,c2->0,c3->0
};

parset=parvec/.pars/.mut->0;
finalfreq[parset]

tNOevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]

fnoevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tNOevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red,Dashed}},AspectRatio->0.75];

fSDashed=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005],Dashed, Blue},AspectRatio->0.75];

mortalitynoevol=yD[tNOevolNC]+yDm[tNOevolNC]/.Flatten[NSolution[parset]];

parset=parvec/.pars/.mut->0.01;
finalfreq[parset]

tevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]
fevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red}},AspectRatio->0.75];

freqERADIC=Evaluate[(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))/.NSolution[parset]]/.t->tevolNC

FIGfreq=Plot[Evaluate[{(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{0,1}},PlotStyle->{{Thickness[0.008],Black}},AspectRatio->0.75];

fS=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005], Blue},AspectRatio->0.75];

fcontrol=Plot[2*c[t]/.pars,{t,0,tfig},Filling->Axis,FillingStyle->LightGray,PlotStyle->None,Frame->True,PlotRange->{{0,tfig},{0,1}},AspectRatio->0.75,Frame->True];


line1=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC-0.02}],Flatten[{tevolNC+5,freqERADIC+0.02}]}]}];
line2=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC+0.02}],Flatten[{tevolNC+5,freqERADIC-0.02}]}]}];


mortalityevol=yD[tevolNC]+yDm[tevolNC]/.Flatten[NSolution[parset]];
(*mortality=Graphics[{Text[Style[StringJoin[ToString[N[Round[100 mortalitynoevol,(10^(-1))] ]],"% (",ToString[N[Round[100 (mortalityevol),(10^(-1))]]],"%)"],12,FontFamily->"Helvetica"],{(t3min+t3max)/2,0.3}]}];*)

DeathRates=Graphics[{Pink,Rectangle[{t3max-tfig/20,0.7},{t3max,0.7+10*mortalitynoevol}],Red,Rectangle[{t3max,0.7},{t3max+tfig/20,0.7+10*mortalityevol}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.2/20,0.705}],Black,Rectangle[{t3max-tfig 1.4/20,0.8},{t3max-tfig 1.2/20,0.805}],Black,Rectangle[{t3max-tfig 1.4/20,0.9},{t3max-tfig 1.2/20,0.905}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.42/20,0.905}]}];

FigTOP4=Show[fS,FIGfreq,fnoevol,fevol,fSDashed,line1,line2,DeathRates,Frame->True,
TicksStyle->Directive[Black,FontFamily->"Helvetica",11],
LabelStyle->Directive[12,FontFamily->"Helvetica",Black],FrameTicks->{{N[{0.2,0.4,0.6,0.8}],N[{0.2,0.4,0.6,0.8}]},{{0,100,200,300,400,500},{0,100,200,300,400,500}}},FrameTicksStyle->{{Directive[FontOpacity->0,FontSize->0],Directive[FontOpacity->0,FontSize->0]},{Directive[FontSize->12],Directive[FontOpacity->0,FontSize->0]}}
]

Figure 3e: Evolution of transmission β, control measures

pars={β->1.3, βm->1.2 1.3, (*TRANSMISSION OF THE MUTANT*)
α->0.005, αm->0.005, (*VIRULENCE OF THE MUTANT*)
κyP->1,κyPm->1, (* 1/TIME IN PRESYMPTOMATIC CLASS FOR THE MUTANT *)

f->0.4, fm->0.4, (* PROPORTION OF ASYMPTOMATIC CASES *)
κyE->0.25, (* 1/TIME IN EXPOSED CLASS *)
κyI->0.2, (* 1/TIME IN SYMPTOMATIC CLASS *)
κyA->0.11, (* 1/TIME IN ASYMPTOMATIC CLASS *)

(*Strength of control measures*)
c1->0.6,c2->0.6,c3->0.6
};

parset=parvec/.pars/.mut->0;
finalfreq[parset]

tNOevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]

fnoevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tNOevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red,Dashed}},AspectRatio->0.75];

fSDashed=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005],Dashed, Blue},AspectRatio->0.75];

mortalitynoevol=yD[tNOevolNC]+yDm[tNOevolNC]/.Flatten[NSolution[parset]];

parset=parvec/.pars/.mut->0.01;
finalfreq[parset]

tevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]
fevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red}},AspectRatio->0.75];

freqERADIC=Evaluate[(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))/.NSolution[parset]]/.t->tevolNC

FIGfreq=Plot[Evaluate[{(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{0,1}},PlotStyle->{{Thickness[0.008],Black}},AspectRatio->0.75];

fS=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005], Blue},AspectRatio->0.75];

fcontrol=Plot[2*c[t]/.pars,{t,0,tfig},Filling->Axis,FillingStyle->LightGray,PlotStyle->None,Frame->True,PlotRange->{{0,tfig},{0,1}},AspectRatio->0.75,Frame->True];

line1=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC-0.02}],Flatten[{tevolNC+5,freqERADIC+0.02}]}]}];
line2=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC+0.02}],Flatten[{tevolNC+5,freqERADIC-0.02}]}]}];

mortalityevol=yD[tevolNC]+yDm[tevolNC]/.Flatten[NSolution[parset]];
(*mortality=Graphics[{Text[Style[StringJoin[ToString[N[Round[100 mortalitynoevol,(10^(-1))] ]],"% (",ToString[N[Round[100 (mortalityevol),(10^(-1))]]],"%)"],12,FontFamily->"Helvetica"],{(t3min+t3max)/2,0.3}]}];*)

DeathRates=Graphics[{Pink,Rectangle[{t3max-tfig/20,0.7},{t3max,0.7+10*mortalitynoevol}],Red,Rectangle[{t3max,0.7},{t3max+tfig/20,0.7+10*mortalityevol}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.2/20,0.705}],Black,Rectangle[{t3max-tfig 1.4/20,0.8},{t3max-tfig 1.2/20,0.805}],Black,Rectangle[{t3max-tfig 1.4/20,0.9},{t3max-tfig 1.2/20,0.905}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.42/20,0.905}]}];

FigBOT1=Show[fcontrol,fS,fSDashed, FIGfreq,fnoevol,fevol,line1,line2,DeathRates,Frame->True,
TicksStyle->Directive[Black,FontFamily->"Helvetica",11],
LabelStyle->Directive[12,FontFamily->"Helvetica",Black],FrameTicks->{{N[{0.2,0.4,0.6,0.8}],N[{0.2,0.4,0.6,0.8}]},{{0,100,200,300,400,500},{0,100,200,300,400,500}}},FrameTicksStyle->{{Directive[FontSize->12],Directive[FontOpacity->0,FontSize->0]},{Directive[FontSize->12],Directive[FontOpacity->0,FontSize->0]}}
]

Figure 3f: Evolution of asymptomatic fraction f, control measures

pars={β->1.3, βm->1.3, (*TRANSMISSION OF THE MUTANT*)
α->0.005, αm->0.005, (*VIRULENCE OF THE MUTANT*)
κyP->1,κyPm->1, (* 1/TIME IN PRESYMPTOMATIC CLASS FOR THE MUTANT *)

f->0.4, fm->0.3, (* PROPORTION OF ASYMPTOMATIC CASES *)
κyE->0.25, (* 1/TIME IN EXPOSED CLASS *)
κyI->0.2, (* 1/TIME IN SYMPTOMATIC CLASS *)
κyA->0.11, (* 1/TIME IN ASYMPTOMATIC CLASS *)

(*Strength of control measures*)
c1->0.6,c2->0.6,c3->0.6
};

parset=parvec/.pars/.mut->0;
finalfreq[parset]

tNOevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]

fnoevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tNOevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red,Dashed}},AspectRatio->0.75];

fSDashed=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005],Dashed, Blue},AspectRatio->0.75];

mortalitynoevol=yD[tNOevolNC]+yDm[tNOevolNC]/.Flatten[NSolution[parset]];

parset=parvec/.pars/.mut->0.01;
finalfreq[parset]

tevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]
fevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red}},AspectRatio->0.75];

freqERADIC=Evaluate[(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))/.NSolution[parset]]/.t->tevolNC

FIGfreq=Plot[Evaluate[{(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{0,1}},PlotStyle->{{Thickness[0.008],Black}},AspectRatio->0.75];

fS=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005], Blue},AspectRatio->0.75];

fcontrol=Plot[2*c[t]/.pars,{t,0,tfig},Filling->Axis,FillingStyle->LightGray,PlotStyle->None,Frame->True,PlotRange->{{0,tfig},{0,1}},AspectRatio->0.75,Frame->True];

line1=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC-0.02}],Flatten[{tevolNC+5,freqERADIC+0.02}]}]}];
line2=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC+0.02}],Flatten[{tevolNC+5,freqERADIC-0.02}]}]}];

mortalityevol=yD[tevolNC]+yDm[tevolNC]/.Flatten[NSolution[parset]];
(*mortality=Graphics[{Text[Style[StringJoin[ToString[N[Round[100 mortalitynoevol,(10^(-1))] ]],"% (",ToString[N[Round[100 (mortalityevol),(10^(-1))]]],"%)"],12,FontFamily->"Helvetica"],{(t3min+t3max)/2,0.3}]}];*)

DeathRates=Graphics[{Pink,Rectangle[{t3max-tfig/20,0.7},{t3max,0.7+10*mortalitynoevol}],Red,Rectangle[{t3max,0.7},{t3max+tfig/20,0.7+10*mortalityevol}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.2/20,0.705}],Black,Rectangle[{t3max-tfig 1.4/20,0.8},{t3max-tfig 1.2/20,0.805}],Black,Rectangle[{t3max-tfig 1.4/20,0.9},{t3max-tfig 1.2/20,0.905}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.42/20,0.905}]}];

FigBOT2=Show[fcontrol,fS,FIGfreq,fnoevol,fevol,fSDashed,line1,line2,DeathRates,Frame->True,
TicksStyle->Directive[Black,FontFamily->"Helvetica",11],
LabelStyle->Directive[12,FontFamily->"Helvetica",Black],FrameTicks->{{N[{0.2,0.4,0.6,0.8}],N[{0.2,0.4,0.6,0.8}]},{{0,100,200,300,400,500},{0,100,200,300,400,500}}},FrameTicksStyle->{{Directive[FontOpacity->0,FontSize->0],Directive[FontOpacity->0,FontSize->0]},{Directive[FontSize->12],Directive[FontOpacity->0,FontSize->0]}}
]

Figure 3g: Evolution of pre-symptomatic phase κ, control measures

pars={β->1.3, βm->1.3, (*TRANSMISSION OF THE MUTANT*)
α->0.005, αm->0.005, (*VIRULENCE OF THE MUTANT*)
κyP->1,κyPm->2/3, (* 1/TIME IN PRESYMPTOMATIC CLASS FOR THE MUTANT *)

f->0.4, fm->0.4, (* PROPORTION OF ASYMPTOMATIC CASES *)
κyE->0.25, (* 1/TIME IN EXPOSED CLASS *)
κyI->0.2, (* 1/TIME IN SYMPTOMATIC CLASS *)
κyA->0.11, (* 1/TIME IN ASYMPTOMATIC CLASS *)

(*Strength of control measures*)
c1->0.6,c2->0.6,c3->0.6
};

parset=parvec/.pars/.mut->0;
finalfreq[parset]

tNOevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]

fnoevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tNOevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red,Dashed}},AspectRatio->0.75];

fSDashed=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005],Dashed, Blue},AspectRatio->0.75];

mortalitynoevol=yD[tNOevolNC]+yDm[tNOevolNC]/.Flatten[NSolution[parset]];

parset=parvec/.pars/.mut->0.01;
finalfreq[parset]

tevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]
fevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red}},AspectRatio->0.75];

freqERADIC=Evaluate[(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))/.NSolution[parset]]/.t->tevolNC

FIGfreq=Plot[Evaluate[{(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{0,1}},PlotStyle->{{Thickness[0.008],Black}},AspectRatio->0.75];

fS=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005], Blue},AspectRatio->0.75];

fcontrol=Plot[2*c[t]/.pars,{t,0,tfig},Filling->Axis,FillingStyle->LightGray,PlotStyle->None,Frame->True,PlotRange->{{0,tfig},{0,1}},AspectRatio->0.75,Frame->True];

line1=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC-0.02}],Flatten[{tevolNC+5,freqERADIC+0.02}]}]}];
line2=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC+0.02}],Flatten[{tevolNC+5,freqERADIC-0.02}]}]}];

mortalityevol=yD[tevolNC]+yDm[tevolNC]/.Flatten[NSolution[parset]];
(*mortality=Graphics[{Text[Style[StringJoin[ToString[N[Round[100 mortalitynoevol,(10^(-1))] ]],"% (",ToString[N[Round[100 (mortalityevol),(10^(-1))]]],"%)"],12,FontFamily->"Helvetica"],{(t3min+t3max)/2,0.3}]}];*)

DeathRates=Graphics[{Pink,Rectangle[{t3max-tfig/20,0.7},{t3max,0.7+10*mortalitynoevol}],Red,Rectangle[{t3max,0.7},{t3max+tfig/20,0.7+10*mortalityevol}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.2/20,0.705}],Black,Rectangle[{t3max-tfig 1.4/20,0.8},{t3max-tfig 1.2/20,0.805}],Black,Rectangle[{t3max-tfig 1.4/20,0.9},{t3max-tfig 1.2/20,0.905}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.42/20,0.905}]}];

FigBOT3=Show[fcontrol,fS,FIGfreq,fnoevol,fevol,fSDashed,line1,line2,DeathRates,Frame->True,
TicksStyle->Directive[Black,FontFamily->"Helvetica",11],
LabelStyle->Directive[12,FontFamily->"Helvetica",Black],FrameTicks->{{N[{0.2,0.4,0.6,0.8}],N[{0.2,0.4,0.6,0.8}]},{{0,100,200,300,400,500},{0,100,200,300,400,500}}},FrameTicksStyle->{{Directive[FontOpacity->0,FontSize->0],Directive[FontOpacity->0,FontSize->0]},{Directive[FontSize->12],Directive[FontOpacity->0,FontSize->0]}}
]

Figure 3h: Evolution of virulence α, control measures

pars={β->1.3, βm->1.3, (*TRANSMISSION OF THE MUTANT*)
α->0.005, αm->0.00, (*VIRULENCE OF THE MUTANT*)
κyP->1,κyPm->1, (* 1/TIME IN PRESYMPTOMATIC CLASS FOR THE MUTANT *)

f->0.4, fm->0.4, (* PROPORTION OF ASYMPTOMATIC CASES *)
κyE->0.25, (* 1/TIME IN EXPOSED CLASS *)
κyI->0.2, (* 1/TIME IN SYMPTOMATIC CLASS *)
κyA->0.11, (* 1/TIME IN ASYMPTOMATIC CLASS *)

(*Strength of control measures*)
c1->0.6,c2->0.6,c3->0.6
};

parset=parvec/.pars/.mut->0;
finalfreq[parset]

tNOevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]

fnoevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tNOevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red,Dashed}},AspectRatio->0.75];

fSDashed=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005],Dashed, Blue},AspectRatio->0.75];

mortalitynoevol=yD[tNOevolNC]+yDm[tNOevolNC]/.Flatten[NSolution[parset]];

parset=parvec/.pars/.mut->0.01;
finalfreq[parset]

tevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]
fevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red}},AspectRatio->0.75];

freqERADIC=Evaluate[(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))/.NSolution[parset]]/.t->tevolNC

FIGfreq=Plot[Evaluate[{(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{0,1}},PlotStyle->{{Thickness[0.008],Black}},AspectRatio->0.75];

fS=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005], Blue},AspectRatio->0.75];

fcontrol=Plot[2*c[t]/.pars,{t,0,tfig},Filling->Axis,FillingStyle->LightGray,PlotStyle->None,Frame->True,PlotRange->{{0,tfig},{0,1}},AspectRatio->0.75,Frame->True];

line1=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC-0.02}],Flatten[{tevolNC+5,freqERADIC+0.02}]}]}];
line2=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC+0.02}],Flatten[{tevolNC+5,freqERADIC-0.02}]}]}];

mortalityevol=yD[tevolNC]+yDm[tevolNC]/.Flatten[NSolution[parset]];
(*mortality=Graphics[{Text[Style[StringJoin[ToString[N[Round[100 mortalitynoevol,(10^(-1))] ]],"% (",ToString[N[Round[100 (mortalityevol),(10^(-1))]]],"%)"],12,FontFamily->"Helvetica"],{(t3min+t3max)/2,0.3}]}];*)

DeathRates=Graphics[{Pink,Rectangle[{t3max-tfig/20,0.7},{t3max,0.7+10*mortalitynoevol}],Red,Rectangle[{t3max,0.7},{t3max+tfig/20,0.7+10*mortalityevol}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.2/20,0.705}],Black,Rectangle[{t3max-tfig 1.4/20,0.8},{t3max-tfig 1.2/20,0.805}],Black,Rectangle[{t3max-tfig 1.4/20,0.9},{t3max-tfig 1.2/20,0.905}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.42/20,0.905}]}];

FigBOT4=Show[fcontrol,fS,FIGfreq,fnoevol,fevol,fSDashed,line1,line2,DeathRates,Frame->True,
TicksStyle->Directive[Black,FontFamily->"Helvetica",11],
LabelStyle->Directive[12,FontFamily->"Helvetica",Black],FrameTicks->{{N[{0.2,0.4,0.6,0.8}],N[{0.2,0.4,0.6,0.8}]},{{0,100,200,300,400,500},{0,100,200,300,400,500}}},FrameTicksStyle->{{Directive[FontOpacity->0,FontSize->0],Directive[FontOpacity->0,FontSize->0]},{Directive[FontSize->12],Directive[FontOpacity->0,FontSize->0]}}
]

Altogether

Show[
GraphicsGrid[{{FigTOP1,FigTOP2,FigTOP3,FigTOP4},{FigBOT1,FigBOT2,FigBOT3,FigBOT4}},FrameLabel->label,Spacings->{-10,10}],
Graphics[Text[Style["(a)",16,FontFamily->"Helvetica"],{50,-33}]],
Graphics[Text[Style["(b)",16,FontFamily->"Helvetica"],{384,-33}]],
Graphics[Text[Style["(c)",16,FontFamily->"Helvetica"],{734,-33}]],
Graphics[Text[Style["(d)",16,FontFamily->"Helvetica"],{1084,-33}]],
Graphics[Text[Style["(e)",16,FontFamily->"Helvetica"],{50,-315}]],
Graphics[Text[Style["(f)",16,FontFamily->"Helvetica"],{384,-315}]],
Graphics[Text[Style["(g)",16,FontFamily->"Helvetica"],{734,-315}]],
Graphics[Text[Style["(h)",16,FontFamily->"Helvetica"],{1084,-315}]],
Graphics[Text[Style["Transmission",16,FontFamily->"Helvetica"],{180,10}]],
Graphics[Text[Style["Asymptomatic fraction",16,FontFamily->"Helvetica"],{ 520,10}]],
Graphics[Text[Style["Pre-symptomatic phase",16,FontFamily->"Helvetica"],{870,10}]],
Graphics[Text[Style["Virulence",16,FontFamily->"Helvetica"],{1220,10}]],
Graphics[Text[Style["Time (days)",16,FontFamily->"Helvetica"],{700,-580}]],
Graphics[Rotate[Text[Style["Frequency",16,FontFamily->"Helvetica"],{-20,-270}],Pi/2]]
]

Export["Fig3.pdf",%];

Figure 3alt - Main change is more transmission in infectious stage (βI = βP)

β adjusted to keep R0 and doubling time of the disease within range observed, and mutant effect on κP increased to visualize.

Code

Common settings

Clear["Global`*"]

Off[General::spell1]
Off[General::spell]
Off[NDSolve::nlnum]

tfinal=500;

tfig=500;

label={"Time (days)","Frequency"};
label=None;

parvec={β, βm,α, αm,κyP,κyPm,f,fm,κyE, κyI,κyA, c1,c2,c3,mut};

Clear[finalfreq,NSolution]
finalfreq[vec_]:=finalfreq[vec]=Block[{β=vec[[1]], βm=vec[[2]],α=vec[[3]], αm=vec[[4]],κyP=vec[[5]],κyPm=vec[[6]],f=vec[[7]],fm=vec[[8]],κyE=vec[[9]], κyI=vec[[10]],κyA=vec[[11]],c1=vec[[12]],c2=vec[[13]],c3=vec[[14]],mut=vec[[15]]},

(*REMAINING PARAMETERS*)
βyA=β/10;βyP=β;βyI=β;
κyAm=κyA;κyIm=κyI;
βyAm=βm/10;βyPm=βm;βyIm=βm;

(* E: Exposed (non-infectious) cases *)
(* A: asymptomatic cases *)
(* P: pre-symptomatic cases *)
(* S: symptomatic cases *)
(* f: proportion of fsymptomatic hosts *)
(* c: amount of distancing applied *)

(*INOCULATION*)
start=(10^(-5));

t=.;
c=.;

c[t_]:=Evaluate[If[t<t1min,0,If[t<t1max,c1,If[t<t2min,0,If[t<t2max,c2,If[t<t3min,0,If[t<t3max,c3,0]]]]]]];
h[t]=(1-c[t])(βyA yA[t]+βyP yP[t]+βyI yI[t]);
hm[t]=(1-c[t])(βyAm yAm[t]+βyPm yPm[t]+βyIm yIm[t]);

sys={S'[t]==-(h[t]+hm[t]) S[t],
yE'[t]==h[t] S[t]-κyE yE[t],
yA'[t]==f κyE yE[t]-κyA yA[t],
yP'[t]==(1-f) κyE yE[t]-κyP yP[t],
yI'[t]==κyP yP[t]-(κyI+α) yI[t],
yR'[t]==κyA yA[t]+κyI yI[t],
yD'[t]==α yI[t]};  (*yD serves as a counter of the cumulative number of deaths*)

sysm={
yEm'[t]==hm[t] S[t]-κyE yEm[t],
yAm'[t]==fm κyE yEm[t]-κyAm yAm[t],
yPm'[t]==(1-fm) κyE yEm[t]-κyPm yPm[t],
yIm'[t]==κyPm yPm[t]-(κyIm+αm) yIm[t],
yRm'[t]==κyAm yAm[t]+κyIm yIm[t],
yDm'[t]==αm yIm[t]};  (*yDm serves as a counter of the cumulative number of deaths*)

(*FIRST BOUT OF CONTROL*)
t1min=75;
t1max=150;
(*SECOND BOUT OF CONTROL*)
t2min=225;
t2max=300;
(*THIRD BOUT OF CONTROL*)
t3min=375;
t3max=450;

init={S[0]==1,yE[0]==start(1-mut),yA[0]==0,yP[0]==0,yI[0]==0,yR[0]==0,yD[0]==0};
initm={yEm[0]==start(mut),yAm[0]==0,yPm[0]==0,yIm[0]==0,yRm[0]==0,yDm[0]==0};

var={S,yE,yA,yP,yI,yR,yD};
varm={yEm,yAm,yPm,yIm,yRm,yDm};

NSolution[vec]=
NDSolve[
Flatten[{sys,sysm,init,initm}],Flatten[{var,varm}],{t,0,tfinal},AccuracyGoal->100];

(*R0*)
Print["R0 of the wild type = ",((βyA f)/(κyA))+(((1-f) (βyP (α+κyI)+βyI κyP))/((α+κyI) κyP))];

(*Rm*)
Print["R0 of the mutant type = ",((βyAm f)/(κyAm))+(((1-f) (βyPm (αm+κyIm)+βyIm κyPm))/((αm+κyIm) κyPm))];

(*CASE MORTALITY*)
Print["Case mortality = ",(1-f)(α/(κyI+α))];
Print["Case mortality (mutant) = ",(1-fm)((αm)/(κyI+αm))];

(*CUMULATIVE MORTALITY*)
Print["Cumulative mortality (total) = ",yD[tfinal]+yDm[tfinal]/.Flatten[NSolution[vec]]];
Print["Cumulative mortality (mutant) = ",yDm[tfinal]/.Flatten[NSolution[vec]]];

]

Figure 3a: Evolution of transmission β, no control measures

pars={β->0.53, βm->1.2 0.53, (*TRANSMISSION OF THE MUTANT*)
α->0.005, αm->0.005, (*VIRULENCE OF THE MUTANT*)
κyP->1,κyPm->1, (* 1/TIME IN PRESYMPTOMATIC CLASS FOR THE MUTANT *)

f->0.2, fm->0.2, (* PROPORTION OF ASYMPTOMATIC CASES *)
κyE->0.25, (* 1/TIME IN EXPOSED CLASS *)
κyI->0.2, (* 1/TIME IN SYMPTOMATIC CLASS *)
κyA->0.11, (* 1/TIME IN ASYMPTOMATIC CLASS *)

(*Strength of control measures*)
c1->0,c2->0,c3->0
};

The doublng time of the disease when rare, as calculated from the eigenvalue (“poly”) of the stability matrix:

Sort[Solve[(κE (-S α βP κA+f S α βP κA-S βP κA κI+f S βP κA κI-f S α βA κP+α κA κP-S βI κA κP+f S βI κA κP-f S βA κI κP+κA κI κP)+(-f S α βA κE-S α βP κE+f S α βP κE+α κA κE-S βP κA κE+f S βP κA κE-f S βA κE κI-S βP κE κI+f S βP κE κI+κA κE κI+α κA κP+α κE κP-f S βA κE κP-S βI κE κP+f S βI κE κP+κA κE κP+κA κI κP+κE κI κP) λ+(α κA+α κE-f S βA κE-S βP κE+f S βP κE+κA κE+κA κI+κE κI+α κP+κA κP+κE κP+κI κP) (λ^2)+(α+κA+κE+κI+κP) (λ^3)+(λ^4)/.S->1/.βA->β/10/.βP->β/.βI->β/3/.κA->κyA/.κE->κyE/.κP->κyP/.κI->κyI/.pars)==0,λ]]
Solve[(Exp[λ t]==2)/.Last[%],t]

parset=parvec/.pars/.mut->0;
finalfreq[parset]

tNOevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]

fnoevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tNOevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red,Dashed}},AspectRatio->0.75];

fSDashed=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005],Dashed, Blue},AspectRatio->0.75];

mortalitynoevol=yD[tNOevolNC]+yDm[tNOevolNC]/.Flatten[NSolution[parset]];

parset=parvec/.pars/.mut->0.01;
finalfreq[parset]

tevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]
fevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red}},AspectRatio->0.75];

freqERADIC=Evaluate[(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))/.NSolution[parset]]/.t->tevolNC

FIGfreq=Plot[Evaluate[{(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{0,1}},PlotStyle->{{Thickness[0.008],Black}},AspectRatio->0.75];

fS=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005], Blue},AspectRatio->0.75];

fcontrol=Plot[2*c[t]/.pars,{t,0,tfig},Filling->Axis,FillingStyle->LightGray,PlotStyle->None,Frame->True,PlotRange->{{0,tfig},{0,1}},AspectRatio->0.75,Frame->True];


line1=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC-0.02}],Flatten[{tevolNC+5,freqERADIC+0.02}]}]}];
line2=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC+0.02}],Flatten[{tevolNC+5,freqERADIC-0.02}]}]}];

mortalityevol=yD[tevolNC]+yDm[tevolNC]/.Flatten[NSolution[parset]];
(*mortality=Graphics[{Text[Style[StringJoin[ToString[N[Round[100 mortalitynoevol,(10^(-1))] ]],"% (",ToString[N[Round[100 (mortalityevol),(10^(-1))]]],"%)"],12,FontFamily->"Helvetica"],{(t3min+t3max)/2,0.3}]}];*)

DeathRates=Graphics[{Pink,Rectangle[{t3max-tfig/20,0.7},{t3max,0.7+10*mortalitynoevol}],Red,Rectangle[{t3max,0.7},{t3max+tfig/20,0.7+10*mortalityevol}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.2/20,0.705}],Black,Rectangle[{t3max-tfig 1.4/20,0.8},{t3max-tfig 1.2/20,0.805}],Black,Rectangle[{t3max-tfig 1.4/20,0.9},{t3max-tfig 1.2/20,0.905}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.42/20,0.905}]}];

FigTOP1=Show[fS,FIGfreq,fnoevol,fevol,fSDashed,line1,line2,DeathRates,Frame->True,
TicksStyle->Directive[Black,FontFamily->"Helvetica",11],
LabelStyle->Directive[12,FontFamily->"Helvetica",Black],
FrameTicks->{{N[{0.2,0.4,0.6,0.8}],N[{0.2,0.4,0.6,0.8}]},{{0,100,200,300,400,500},{0,100,200,300,400,500}}},FrameTicksStyle->{{Directive[FontSize->12],Directive[FontOpacity->0,FontSize->0]},{Directive[FontSize->12],Directive[FontOpacity->0,FontSize->0]}}
]

Figure 3b: Evolution of asymptomatic fraction f,  no control measures

pars={β->0.53, βm->0.53, (*TRANSMISSION OF THE MUTANT*)
α->0.005, αm->0.005, (*VIRULENCE OF THE MUTANT*)
κyP->1,κyPm->1, (* 1/TIME IN PRESYMPTOMATIC CLASS FOR THE MUTANT *)

f->0.2, fm->0.1, (* PROPORTION OF ASYMPTOMATIC CASES *)
κyE->0.25, (* 1/TIME IN EXPOSED CLASS *)
κyI->0.2, (* 1/TIME IN SYMPTOMATIC CLASS *)
κyA->0.11, (* 1/TIME IN ASYMPTOMATIC CLASS *)

(*Strength of control measures*)
c1->0,c2->0,c3->0
};

parset=parvec/.pars/.mut->0;
finalfreq[parset]

tNOevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]

fnoevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tNOevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red,Dashed}},AspectRatio->0.75];

fSDashed=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005],Dashed, Blue},AspectRatio->0.75];

mortalitynoevol=yD[tNOevolNC]+yDm[tNOevolNC]/.Flatten[NSolution[parset]];

parset=parvec/.pars/.mut->0.01;
finalfreq[parset]

tevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]
fevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red}},AspectRatio->0.75];

freqERADIC=Evaluate[(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))/.NSolution[parset]]/.t->tevolNC

FIGfreq=Plot[Evaluate[{(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{0,1}},PlotStyle->{{Thickness[0.008],Black}},AspectRatio->0.75];

fS=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005], Blue},AspectRatio->0.75];

fcontrol=Plot[2*c[t]/.pars,{t,0,tfig},Filling->Axis,FillingStyle->LightGray,PlotStyle->None,Frame->True,PlotRange->{{0,tfig},{0,1}},AspectRatio->0.75,Frame->True];


line1=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC-0.02}],Flatten[{tevolNC+5,freqERADIC+0.02}]}]}];
line2=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC+0.02}],Flatten[{tevolNC+5,freqERADIC-0.02}]}]}];


mortalityevol=yD[tevolNC]+yDm[tevolNC]/.Flatten[NSolution[parset]];
(*mortality=Graphics[{Text[Style[StringJoin[ToString[N[Round[100 mortalitynoevol,(10^(-1))] ]],"% (",ToString[N[Round[100 (mortalityevol),(10^(-1))]]],"%)"],12,FontFamily->"Helvetica"],{(t3min+t3max)/2,0.3}]}];*)

DeathRates=Graphics[{Pink,Rectangle[{t3max-tfig/20,0.7},{t3max,0.7+10*mortalitynoevol}],Red,Rectangle[{t3max,0.7},{t3max+tfig/20,0.7+10*mortalityevol}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.2/20,0.705}],Black,Rectangle[{t3max-tfig 1.4/20,0.8},{t3max-tfig 1.2/20,0.805}],Black,Rectangle[{t3max-tfig 1.4/20,0.9},{t3max-tfig 1.2/20,0.905}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.42/20,0.905}]}];

FigTOP2=Show[fS,FIGfreq,fnoevol,fevol,fSDashed,line1,line2,DeathRates,Frame->True,
TicksStyle->Directive[Black,FontFamily->"Helvetica",11],
LabelStyle->Directive[12,FontFamily->"Helvetica",Black],FrameTicks->{{N[{0.2,0.4,0.6,0.8}],N[{0.2,0.4,0.6,0.8}]},{{0,100,200,300,400,500},{0,100,200,300,400,500}}},FrameTicksStyle->{{Directive[FontOpacity->0,FontSize->0],Directive[FontOpacity->0,FontSize->0]},{Directive[FontSize->12],Directive[FontOpacity->0,FontSize->0]}}
]

Figure 3c: Evolution of pre-symptomatic phase κ, no control measures

pars={β->0.53, βm->0.53, (*TRANSMISSION OF THE MUTANT*)
α->0.005, αm->0.005, (*VIRULENCE OF THE MUTANT*)
κyP->1,κyPm->1/3, (* 1/TIME IN PRESYMPTOMATIC CLASS FOR THE MUTANT *)

f->0.2, fm->0.2, (* PROPORTION OF ASYMPTOMATIC CASES *)
κyE->0.25, (* 1/TIME IN EXPOSED CLASS *)
κyI->0.2, (* 1/TIME IN SYMPTOMATIC CLASS *)
κyA->0.11, (* 1/TIME IN ASYMPTOMATIC CLASS *)

(*Strength of control measures*)
c1->0,c2->0,c3->0
};

parset=parvec/.pars/.mut->0;
finalfreq[parset]

tNOevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]

fnoevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tNOevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red,Dashed}},AspectRatio->0.75];

fSDashed=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005],Dashed, Blue},AspectRatio->0.75];

mortalitynoevol=yD[tNOevolNC]+yDm[tNOevolNC]/.Flatten[NSolution[parset]];

parset=parvec/.pars/.mut->0.01;
finalfreq[parset]

tevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]
fevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red}},AspectRatio->0.75];

freqERADIC=Evaluate[(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))/.NSolution[parset]]/.t->tevolNC

FIGfreq=Plot[Evaluate[{(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{0,1}},PlotStyle->{{Thickness[0.008],Black}},AspectRatio->0.75];

fS=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005], Blue},AspectRatio->0.75];

fcontrol=Plot[2*c[t]/.pars,{t,0,tfig},Filling->Axis,FillingStyle->LightGray,PlotStyle->None,Frame->True,PlotRange->{{0,tfig},{0,1}},AspectRatio->0.75,Frame->True];


line1=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC-0.02}],Flatten[{tevolNC+5,freqERADIC+0.02}]}]}];
line2=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC+0.02}],Flatten[{tevolNC+5,freqERADIC-0.02}]}]}];


mortalityevol=yD[tevolNC]+yDm[tevolNC]/.Flatten[NSolution[parset]];
(*mortality=Graphics[{Text[Style[StringJoin[ToString[N[Round[100 mortalitynoevol,(10^(-1))] ]],"% (",ToString[N[Round[100 (mortalityevol),(10^(-1))]]],"%)"],12,FontFamily->"Helvetica"],{(t3min+t3max)/2,0.3}]}];*)

DeathRates=Graphics[{Pink,Rectangle[{t3max-tfig/20,0.7},{t3max,0.7+10*mortalitynoevol}],Red,Rectangle[{t3max,0.7},{t3max+tfig/20,0.7+10*mortalityevol}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.2/20,0.705}],Black,Rectangle[{t3max-tfig 1.4/20,0.8},{t3max-tfig 1.2/20,0.805}],Black,Rectangle[{t3max-tfig 1.4/20,0.9},{t3max-tfig 1.2/20,0.905}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.42/20,0.905}]}];

FigTOP3=Show[fS,FIGfreq,fnoevol,fevol,fSDashed,line1,line2,DeathRates,Frame->True,
TicksStyle->Directive[Black,FontFamily->"Helvetica",11],
LabelStyle->Directive[12,FontFamily->"Helvetica",Black],FrameTicks->{{N[{0.2,0.4,0.6,0.8}],N[{0.2,0.4,0.6,0.8}]},{{0,100,200,300,400,500},{0,100,200,300,400,500}}},FrameTicksStyle->{{Directive[FontOpacity->0,FontSize->0],Directive[FontOpacity->0,FontSize->0]},{Directive[FontSize->12],Directive[FontOpacity->0,FontSize->0]}}
]

Figure 3d: Evolution of virulence α, no control measures

pars={β->0.53, βm->0.53, (*TRANSMISSION OF THE MUTANT*)
α->0.005, αm->0.0, (*VIRULENCE OF THE MUTANT*)
κyP->1,κyPm->1, (* 1/TIME IN PRESYMPTOMATIC CLASS FOR THE MUTANT *)

f->0.2, fm->0.2, (* PROPORTION OF ASYMPTOMATIC CASES *)
κyE->0.25, (* 1/TIME IN EXPOSED CLASS *)
κyI->0.2, (* 1/TIME IN SYMPTOMATIC CLASS *)
κyA->0.11, (* 1/TIME IN ASYMPTOMATIC CLASS *)

(*Strength of control measures*)
c1->0,c2->0,c3->0
};

parset=parvec/.pars/.mut->0;
finalfreq[parset]

tNOevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]

fnoevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tNOevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red,Dashed}},AspectRatio->0.75];

fSDashed=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005],Dashed, Blue},AspectRatio->0.75];

mortalitynoevol=yD[tNOevolNC]+yDm[tNOevolNC]/.Flatten[NSolution[parset]];

parset=parvec/.pars/.mut->0.01;
finalfreq[parset]

tevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]
fevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red}},AspectRatio->0.75];

freqERADIC=Evaluate[(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))/.NSolution[parset]]/.t->tevolNC

FIGfreq=Plot[Evaluate[{(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{0,1}},PlotStyle->{{Thickness[0.008],Black}},AspectRatio->0.75];

fS=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005], Blue},AspectRatio->0.75];

fcontrol=Plot[2*c[t]/.pars,{t,0,tfig},Filling->Axis,FillingStyle->LightGray,PlotStyle->None,Frame->True,PlotRange->{{0,tfig},{0,1}},AspectRatio->0.75,Frame->True];


line1=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC-0.02}],Flatten[{tevolNC+5,freqERADIC+0.02}]}]}];
line2=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC+0.02}],Flatten[{tevolNC+5,freqERADIC-0.02}]}]}];


mortalityevol=yD[tevolNC]+yDm[tevolNC]/.Flatten[NSolution[parset]];
(*mortality=Graphics[{Text[Style[StringJoin[ToString[N[Round[100 mortalitynoevol,(10^(-1))] ]],"% (",ToString[N[Round[100 (mortalityevol),(10^(-1))]]],"%)"],12,FontFamily->"Helvetica"],{(t3min+t3max)/2,0.3}]}];*)

DeathRates=Graphics[{Pink,Rectangle[{t3max-tfig/20,0.7},{t3max,0.7+10*mortalitynoevol}],Red,Rectangle[{t3max,0.7},{t3max+tfig/20,0.7+10*mortalityevol}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.2/20,0.705}],Black,Rectangle[{t3max-tfig 1.4/20,0.8},{t3max-tfig 1.2/20,0.805}],Black,Rectangle[{t3max-tfig 1.4/20,0.9},{t3max-tfig 1.2/20,0.905}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.42/20,0.905}]}];

FigTOP4=Show[fS,FIGfreq,fnoevol,fevol,fSDashed,line1,line2,DeathRates,Frame->True,
TicksStyle->Directive[Black,FontFamily->"Helvetica",11],
LabelStyle->Directive[12,FontFamily->"Helvetica",Black],FrameTicks->{{N[{0.2,0.4,0.6,0.8}],N[{0.2,0.4,0.6,0.8}]},{{0,100,200,300,400,500},{0,100,200,300,400,500}}},FrameTicksStyle->{{Directive[FontOpacity->0,FontSize->0],Directive[FontOpacity->0,FontSize->0]},{Directive[FontSize->12],Directive[FontOpacity->0,FontSize->0]}}
]

Figure 3e: Evolution of transmission β, control measures

pars={β->0.53, βm->0.53 1.2, (*TRANSMISSION OF THE MUTANT*)
α->0.005, αm->0.005, (*VIRULENCE OF THE MUTANT*)
κyP->1,κyPm->1, (* 1/TIME IN PRESYMPTOMATIC CLASS FOR THE MUTANT *)

f->0.2, fm->0.2, (* PROPORTION OF ASYMPTOMATIC CASES *)
κyE->0.25, (* 1/TIME IN EXPOSED CLASS *)
κyI->0.2, (* 1/TIME IN SYMPTOMATIC CLASS *)
κyA->0.11, (* 1/TIME IN ASYMPTOMATIC CLASS *)

(*Strength of control measures*)
c1->0.6,c2->0.6,c3->0.6
};

parset=parvec/.pars/.mut->0;
finalfreq[parset]

tNOevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]

fnoevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tNOevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red,Dashed}},AspectRatio->0.75];

fSDashed=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005],Dashed, Blue},AspectRatio->0.75];

mortalitynoevol=yD[tNOevolNC]+yDm[tNOevolNC]/.Flatten[NSolution[parset]];

parset=parvec/.pars/.mut->0.01;
finalfreq[parset]

tevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]
fevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red}},AspectRatio->0.75];

freqERADIC=Evaluate[(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))/.NSolution[parset]]/.t->tevolNC

FIGfreq=Plot[Evaluate[{(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{0,1}},PlotStyle->{{Thickness[0.008],Black}},AspectRatio->0.75];

fS=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005], Blue},AspectRatio->0.75];

fcontrol=Plot[2*c[t]/.pars,{t,0,tfig},Filling->Axis,FillingStyle->LightGray,PlotStyle->None,Frame->True,PlotRange->{{0,tfig},{0,1}},AspectRatio->0.75,Frame->True];

line1=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC-0.02}],Flatten[{tevolNC+5,freqERADIC+0.02}]}]}];
line2=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC+0.02}],Flatten[{tevolNC+5,freqERADIC-0.02}]}]}];

mortalityevol=yD[tevolNC]+yDm[tevolNC]/.Flatten[NSolution[parset]];
(*mortality=Graphics[{Text[Style[StringJoin[ToString[N[Round[100 mortalitynoevol,(10^(-1))] ]],"% (",ToString[N[Round[100 (mortalityevol),(10^(-1))]]],"%)"],12,FontFamily->"Helvetica"],{(t3min+t3max)/2,0.3}]}];*)

DeathRates=Graphics[{Pink,Rectangle[{t3max-tfig/20,0.7},{t3max,0.7+10*mortalitynoevol}],Red,Rectangle[{t3max,0.7},{t3max+tfig/20,0.7+10*mortalityevol}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.2/20,0.705}],Black,Rectangle[{t3max-tfig 1.4/20,0.8},{t3max-tfig 1.2/20,0.805}],Black,Rectangle[{t3max-tfig 1.4/20,0.9},{t3max-tfig 1.2/20,0.905}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.42/20,0.905}]}];

FigBOT1=Show[fcontrol,fS,fSDashed, FIGfreq,fnoevol,fevol,line1,line2,DeathRates,Frame->True,
TicksStyle->Directive[Black,FontFamily->"Helvetica",11],
LabelStyle->Directive[12,FontFamily->"Helvetica",Black],FrameTicks->{{N[{0.2,0.4,0.6,0.8}],N[{0.2,0.4,0.6,0.8}]},{{0,100,200,300,400,500},{0,100,200,300,400,500}}},FrameTicksStyle->{{Directive[FontSize->12],Directive[FontOpacity->0,FontSize->0]},{Directive[FontSize->12],Directive[FontOpacity->0,FontSize->0]}}
]

Figure 3f: Evolution of asymptomatic fraction f, control measures

pars={β->0.53, βm->0.53, (*TRANSMISSION OF THE MUTANT*)
α->0.005, αm->0.005, (*VIRULENCE OF THE MUTANT*)
κyP->1,κyPm->1, (* 1/TIME IN PRESYMPTOMATIC CLASS FOR THE MUTANT *)

f->0.2, fm->0.1, (* PROPORTION OF ASYMPTOMATIC CASES *)
κyE->0.25, (* 1/TIME IN EXPOSED CLASS *)
κyI->0.2, (* 1/TIME IN SYMPTOMATIC CLASS *)
κyA->0.11, (* 1/TIME IN ASYMPTOMATIC CLASS *)

(*Strength of control measures*)
c1->0.6,c2->0.6,c3->0.6
};

parset=parvec/.pars/.mut->0;
finalfreq[parset]

tNOevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]

fnoevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tNOevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red,Dashed}},AspectRatio->0.75];

fSDashed=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005],Dashed, Blue},AspectRatio->0.75];

mortalitynoevol=yD[tNOevolNC]+yDm[tNOevolNC]/.Flatten[NSolution[parset]];

parset=parvec/.pars/.mut->0.01;
finalfreq[parset]

tevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]
fevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red}},AspectRatio->0.75];

freqERADIC=Evaluate[(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))/.NSolution[parset]]/.t->tevolNC

FIGfreq=Plot[Evaluate[{(yEm[t]+yAm[t]+yPm[t]+yIm[t])/(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{0,1}},PlotStyle->{{Thickness[0.008],Black}},AspectRatio->0.75];

fS=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005], Blue},AspectRatio->0.75];

fcontrol=Plot[2*c[t]/.pars,{t,0,tfig},Filling->Axis,FillingStyle->LightGray,PlotStyle->None,Frame->True,PlotRange->{{0,tfig},{0,1}},AspectRatio->0.75,Frame->True];

line1=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC-0.02}],Flatten[{tevolNC+5,freqERADIC+0.02}]}]}];
line2=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC+0.02}],Flatten[{tevolNC+5,freqERADIC-0.02}]}]}];

mortalityevol=yD[tevolNC]+yDm[tevolNC]/.Flatten[NSolution[parset]];
(*mortality=Graphics[{Text[Style[StringJoin[ToString[N[Round[100 mortalitynoevol,(10^(-1))] ]],"% (",ToString[N[Round[100 (mortalityevol),(10^(-1))]]],"%)"],12,FontFamily->"Helvetica"],{(t3min+t3max)/2,0.3}]}];*)

DeathRates=Graphics[{Pink,Rectangle[{t3max-tfig/20,0.7},{t3max,0.7+10*mortalitynoevol}],Red,Rectangle[{t3max,0.7},{t3max+tfig/20,0.7+10*mortalityevol}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.2/20,0.705}],Black,Rectangle[{t3max-tfig 1.4/20,0.8},{t3max-tfig 1.2/20,0.805}],Black,Rectangle[{t3max-tfig 1.4/20,0.9},{t3max-tfig 1.2/20,0.905}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.42/20,0.905}]}];

FigBOT2=Show[fcontrol,fS,FIGfreq,fnoevol,fevol,fSDashed,line1,line2,DeathRates,Frame->True,
TicksStyle->Directive[Black,FontFamily->"Helvetica",11],
LabelStyle->Directive[12,FontFamily->"Helvetica",Black],FrameTicks->{{N[{0.2,0.4,0.6,0.8}],N[{0.2,0.4,0.6,0.8}]},{{0,100,200,300,400,500},{0,100,200,300,400,500}}},FrameTicksStyle->{{Directive[FontOpacity->0,FontSize->0],Directive[FontOpacity->0,FontSize->0]},{Directive[FontSize->12],Directive[FontOpacity->0,FontSize->0]}}
]

Figure 3g: Evolution of pre-symptomatic phase κ, control measures

pars={β->0.53, βm->0.53, (*TRANSMISSION OF THE MUTANT*)
α->0.005, αm->0.005, (*VIRULENCE OF THE MUTANT*)
κyP->1,κyPm->1/3, (* 1/TIME IN PRESYMPTOMATIC CLASS FOR THE MUTANT *)

f->0.2, fm->0.2, (* PROPORTION OF ASYMPTOMATIC CASES *)
κyE->0.25, (* 1/TIME IN EXPOSED CLASS *)
κyI->0.2, (* 1/TIME IN SYMPTOMATIC CLASS *)
κyA->0.11, (* 1/TIME IN ASYMPTOMATIC CLASS *)

(*Strength of control measures*)
c1->0.6,c2->0.6,c3->0.6
};

parset=parvec/.pars/.mut->0;
finalfreq[parset]

tNOevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]

fnoevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tNOevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red,Dashed}},AspectRatio->0.75];

fSDashed=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005],Dashed, Blue},AspectRatio->0.75];

mortalitynoevol=yD[tNOevolNC]+yDm[tNOevolNC]/.Flatten[NSolution[parset]];

parset=parvec/.pars/.mut->0.01;
finalfreq[parset]

tevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]
fevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red}},AspectRatio->0.75];

freqERADIC=Evaluate[(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))/.NSolution[parset]]/.t->tevolNC

FIGfreq=Plot[Evaluate[{(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{0,1}},PlotStyle->{{Thickness[0.008],Black}},AspectRatio->0.75];

fS=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005], Blue},AspectRatio->0.75];

fcontrol=Plot[2*c[t]/.pars,{t,0,tfig},Filling->Axis,FillingStyle->LightGray,PlotStyle->None,Frame->True,PlotRange->{{0,tfig},{0,1}},AspectRatio->0.75,Frame->True];

line1=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC-0.02}],Flatten[{tevolNC+5,freqERADIC+0.02}]}]}];
line2=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC+0.02}],Flatten[{tevolNC+5,freqERADIC-0.02}]}]}];

mortalityevol=yD[tevolNC]+yDm[tevolNC]/.Flatten[NSolution[parset]];
(*mortality=Graphics[{Text[Style[StringJoin[ToString[N[Round[100 mortalitynoevol,(10^(-1))] ]],"% (",ToString[N[Round[100 (mortalityevol),(10^(-1))]]],"%)"],12,FontFamily->"Helvetica"],{(t3min+t3max)/2,0.3}]}];*)

DeathRates=Graphics[{Pink,Rectangle[{t3max-tfig/20,0.7},{t3max,0.7+10*mortalitynoevol}],Red,Rectangle[{t3max,0.7},{t3max+tfig/20,0.7+10*mortalityevol}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.2/20,0.705}],Black,Rectangle[{t3max-tfig 1.4/20,0.8},{t3max-tfig 1.2/20,0.805}],Black,Rectangle[{t3max-tfig 1.4/20,0.9},{t3max-tfig 1.2/20,0.905}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.42/20,0.905}]}];

FigBOT3=Show[fcontrol,fS,FIGfreq,fnoevol,fevol,fSDashed,line1,line2,DeathRates,Frame->True,
TicksStyle->Directive[Black,FontFamily->"Helvetica",11],
LabelStyle->Directive[12,FontFamily->"Helvetica",Black],FrameTicks->{{N[{0.2,0.4,0.6,0.8}],N[{0.2,0.4,0.6,0.8}]},{{0,100,200,300,400,500},{0,100,200,300,400,500}}},FrameTicksStyle->{{Directive[FontOpacity->0,FontSize->0],Directive[FontOpacity->0,FontSize->0]},{Directive[FontSize->12],Directive[FontOpacity->0,FontSize->0]}}
]

Figure 3h: Evolution of virulence α, control measures

pars={β->0.53, βm->0.53, (*TRANSMISSION OF THE MUTANT*)
α->0.005, αm->0.00, (*VIRULENCE OF THE MUTANT*)
κyP->1,κyPm->1, (* 1/TIME IN PRESYMPTOMATIC CLASS FOR THE MUTANT *)

f->0.2, fm->0.2, (* PROPORTION OF ASYMPTOMATIC CASES *)
κyE->0.25, (* 1/TIME IN EXPOSED CLASS *)
κyI->0.2, (* 1/TIME IN SYMPTOMATIC CLASS *)
κyA->0.11, (* 1/TIME IN ASYMPTOMATIC CLASS *)

(*Strength of control measures*)
c1->0.6,c2->0.6,c3->0.6
};

parset=parvec/.pars/.mut->0;
finalfreq[parset]

tNOevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]

fnoevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tNOevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red,Dashed}},AspectRatio->0.75];

fSDashed=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005],Dashed, Blue},AspectRatio->0.75];

mortalitynoevol=yD[tNOevolNC]+yDm[tNOevolNC]/.Flatten[NSolution[parset]];

parset=parvec/.pars/.mut->0.01;
finalfreq[parset]

tevolNC=LengthWhile[Flatten[Table[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,1,tfig,1}]],#>(10^(-7))&]
fevol=Plot[Evaluate[{(yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{10^-7,1}},PlotStyle->{{Thickness[0.005],Red}},AspectRatio->0.75];

freqERADIC=Evaluate[(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))/.NSolution[parset]]/.t->tevolNC

FIGfreq=Plot[Evaluate[{(((yEm[t]+yAm[t]+yPm[t]+yIm[t]))/((yE[t]+yA[t]+yP[t]+yI[t]+yEm[t]+yAm[t]+yPm[t]+yIm[t])))}/.NSolution[parset]],{t,0,tevolNC},PlotRange->{{0,tfig},{0,1}},PlotStyle->{{Thickness[0.008],Black}},AspectRatio->0.75];

fS=Plot[Evaluate[S[t]/.NSolution[parset]],{t,0,tfig},PlotRange->{{0,tfig},{0,1}},PlotStyle->{Thickness[0.005], Blue},AspectRatio->0.75];

fcontrol=Plot[2*c[t]/.pars,{t,0,tfig},Filling->Axis,FillingStyle->LightGray,PlotStyle->None,Frame->True,PlotRange->{{0,tfig},{0,1}},AspectRatio->0.75,Frame->True];

line1=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC-0.02}],Flatten[{tevolNC+5,freqERADIC+0.02}]}]}];
line2=Graphics[{Thickness[0.007],Black,Line[{Flatten[{tevolNC-5,freqERADIC+0.02}],Flatten[{tevolNC+5,freqERADIC-0.02}]}]}];

mortalityevol=yD[tevolNC]+yDm[tevolNC]/.Flatten[NSolution[parset]];
(*mortality=Graphics[{Text[Style[StringJoin[ToString[N[Round[100 mortalitynoevol,(10^(-1))] ]],"% (",ToString[N[Round[100 (mortalityevol),(10^(-1))]]],"%)"],12,FontFamily->"Helvetica"],{(t3min+t3max)/2,0.3}]}];*)

DeathRates=Graphics[{Pink,Rectangle[{t3max-tfig/20,0.7},{t3max,0.7+10*mortalitynoevol}],Red,Rectangle[{t3max,0.7},{t3max+tfig/20,0.7+10*mortalityevol}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.2/20,0.705}],Black,Rectangle[{t3max-tfig 1.4/20,0.8},{t3max-tfig 1.2/20,0.805}],Black,Rectangle[{t3max-tfig 1.4/20,0.9},{t3max-tfig 1.2/20,0.905}],Black,Rectangle[{t3max-tfig 1.4/20,0.7},{t3max-tfig 1.42/20,0.905}]}];

FigBOT4=Show[fcontrol,fS,FIGfreq,fnoevol,fevol,fSDashed,line1,line2,DeathRates,Frame->True,
TicksStyle->Directive[Black,FontFamily->"Helvetica",11],
LabelStyle->Directive[12,FontFamily->"Helvetica",Black],FrameTicks->{{N[{0.2,0.4,0.6,0.8}],N[{0.2,0.4,0.6,0.8}]},{{0,100,200,300,400,500},{0,100,200,300,400,500}}},FrameTicksStyle->{{Directive[FontOpacity->0,FontSize->0],Directive[FontOpacity->0,FontSize->0]},{Directive[FontSize->12],Directive[FontOpacity->0,FontSize->0]}}
]

Altogether

Show[
GraphicsGrid[{{FigTOP1,FigTOP2,FigTOP3,FigTOP4},{FigBOT1,FigBOT2,FigBOT3,FigBOT4}},FrameLabel->label,Spacings->{-10,10}],
Graphics[Text[Style["(a)",16,FontFamily->"Helvetica"],{50,-33}]],
Graphics[Text[Style["(b)",16,FontFamily->"Helvetica"],{384,-33}]],
Graphics[Text[Style["(c)",16,FontFamily->"Helvetica"],{734,-33}]],
Graphics[Text[Style["(d)",16,FontFamily->"Helvetica"],{1084,-33}]],
Graphics[Text[Style["(e)",16,FontFamily->"Helvetica"],{50,-315}]],
Graphics[Text[Style["(f)",16,FontFamily->"Helvetica"],{384,-315}]],
Graphics[Text[Style["(g)",16,FontFamily->"Helvetica"],{734,-315}]],
Graphics[Text[Style["(h)",16,FontFamily->"Helvetica"],{1084,-315}]],
Graphics[Text[Style["Transmission",16,FontFamily->"Helvetica"],{180,10}]],
Graphics[Text[Style["Asymptomatic fraction",16,FontFamily->"Helvetica"],{ 520,10}]],
Graphics[Text[Style["Pre-symptomatic phase",16,FontFamily->"Helvetica"],{870,10}]],
Graphics[Text[Style["Virulence",16,FontFamily->"Helvetica"],{1220,10}]],
Graphics[Text[Style["Time (days)",16,FontFamily->"Helvetica"],{700,-580}]],
Graphics[Rotate[Text[Style["Frequency",16,FontFamily->"Helvetica"],{-20,-270}],Pi/2]]
]