Virulence evolution of pathogens that can grow in reservoir environments

Virulence evolution of pathogens that can grow in reservoir environments
by Aakash Pandey

Article abstract

Many pathogens reside in environmental reservoirs within which they can reproduce and from which they can infect hosts. These facultative pathogens experience different selective pressures in host-associated environments and reservoir environments. Heterogeneous selective pressures have the potential to influence the virulence evolution of these pathogens. Previous research has examined how environmental transmission influences the selective pressures shaping the virulence of pathogens that cannot reproduce in environmental reservoirs, yet many pathogens of humans, crop plants, and livestock can reproduce in these environments. We build on this work to examine how reproduction in reservoirs influences disease dynamics and virulence evolution in a simple facultative pathogen model. We use adaptive dynamics to examine the evolutionary dynamics of facultative pathogens under potential trade-offs between transmission and virulence, shedding and virulence, and reservoir persistence and virulence. We then perform critical function analysis to generalize the results independent of specific trade-off assumptions. We determine that diverse virulence strategies, sometimes resulting from evolutionary bistability or evolutionary branching conditions, are expected for facultative pathogens. Our findings motivate research establishing which trade-offs most strongly influence the virulence evolution of facultative pathogens.

Introduction

This file presents the codes used to generate the figures presented in the accompanying article (Pandey, A., Mideo, N., & Platt, T. G. (2021). Virulence evolution of pathogens that can grow in reservoir environments. The American Naturalist). Here we explore the epidemiological and evolutionary outcomes for a simple facultative pathogen model:

dS

dt

=b-βSI-

γSϕP

(τ+P)

-dS

dI

dt

=βSI+

γSϕP

(τ+P)

-(d+ρ+α)S

dP

dt

=θI-σP+rP1-

P

K

(1)The susceptible hosts (S) are born at a constant rate (b) and die at a per capita rate (d), with direct infection following mass action incidence with rate constant β. The recovery rate is ρ, and we assume lifelong immunity once the host is recovered. The contact rate between susceptible hosts and pathogens in the environmental reservoir is an increasing but saturating function of P, described by ϕP/(P+τ). The term γ represents the probability of environmental transmission on contact with pathogens in the reservoir. Free-living pathogen population dynamics depend on shedding from infected hosts at rate θ, death of free-living pathogen propagules at rate σ, and density-dependent (specifically, logistic) reproduction of pathogen propagules at per capita rate r. K represents the carrying capacity in the reservoir.

Figure 2a. Differences in threshold parameters:
R
0d
vs λ

λ is the demographic reproduction number and is derived from the Jacobian of Infected- Free-living pathogen system (equation 1).
R
0d
represents the number of secondary infections that arise from a single infected individual before the pathogen reproduces in the reservoir. The code below shows the effect of virulence (α) and pathogen growth rate in environment (r) on these two metrics through a contour plot. Disease spreads in population if either of
R
0d
and λ is greater than 1. Please refer to the paper and its appendix section for more details on these metrics.

Clear["Global`*"](* Parameter values *)Pars={b1,d 0.2,γ0.1,σ0.5,θ0.1,β0.1};(* Define λ *)λ=1/2(β*b/(d*(d+α))+r/σ+Sqrt[(β*b/(d*(d+α))-r/σ)^2+4*γ*b/(d*σ)*θ/(d+α)] )/.Pars;(* Define R0 *)R0d=β*b/(d*(d+α))+γ*b/(d*σ)*θ/(d+α) /.Pars;Show[ ContourPlot[ λ,{α,0,2},{r,0,1}, Contours{1}, RegionFunctionFunction[ {α,r,R0d},R0d>1 ], ContourShading{Darker[Orange,.008]} ], ContourPlot[ λ,{α,0,2},{r,0,1}, Contours{1}, RegionFunctionFunction[ {α,r,R0d},R0d<1 ], ContourShading{Lighter[Orange,1/2]} ], ContourPlot[ R0d,{α,0,2},{r,0,1}, Contours{1}, RegionFunctionFunction[ {α,r,R0d},R0d>1 ], ContourShading{Darker[Orange,1/5]} ], FrameLabel{"Virulence (α)","Free-living pathogen birth rate (r)"}, Frame{True,True,False,False}, FrameStyleDirective[Black,15], FrameTicks{{0.5,1.0,1.5,2.0},{0.2,0.4,0.6,0.8,1.0}}, LabelStyleDirective[Black,15], Epilog Text[Style["λ > 1",Medium],{0.25,0.6}], Text[Style[Subscript[R,"0d"]>1,Medium],{0.25,0.55}], Text[Style["λ > 1",Medium],{1,0.6}], Text[Style[Subscript[R,"0d"]<1,Medium],{1,0.55}], Text[Style["λ < 1",Medium],{1,0.25}], Text[Style[Subscript[R,"0d"]<1,Medium],{1,0.2}]  ]

Out[]=

Figure 2b. Effect of free-living pathogens and reservoir carrying capacity in endemic disease prevalence

Clear["Global`*"](*Parameter values *)Pars={b1,γ 0.3,d 0.2,ρ 0,β0.01,θ0.1,σ0.1,α0.3};(* System of equations *)dS= b-(β*S[t]*Inf[t])-(γ*S[t]*(P[t]/(P[t]+1)))-(d*S[t]) /.Pars;dI=(β*S[t]*Inf[t])+(γ*S[t]*(P[t]/(P[t]+1)))-(d+ρ+α)*Inf[t] /.Pars;dP=(θ*Inf[t])-σ P[t]+r P[t] (1-P[t]/Ka) /.Pars;(* Equilibrium solutions *)sol2=Solve[{dS0,dI0,dP0},{S[t],Inf[t],P[t]}];(* Plotting disease prevalence at endemic equilibrium solution *)ContourPlot[ {Inf[t]/(S[t]+Inf[t])}/.sol2[[5]], {r,0.1,0.5}, {Ka,2,20}, FrameLabel{"Growth rate (r)","Carrying capacity (K)"}, LabelStyleDirective[Black,15], FrameStyleDirective[Black,15], FrameTicks {{5,10,15,20},None},{{0.1,0.2,0.3,0.4,0.5},None} , PlotLegendsBarLegend[Automatic,LegendMarkerSize140,LegendFunction"Frame",LegendLabel->"Prevalence"], ColorFunctionGrayLevel ]

Out[]=

Prevalence

Fig 3. PIPs under transmission-virulence trade-off

Here we generate the pairwise-invasibility-plot (PIP) to show the effect of transmission-virulence trade-off. You can change the value of h to change the strength of trade-off.

Clear["Global`*"]Pars={b1,γ 0.35,d 0.2,ρ 0, v 0.02, u3,g1,f1,c1,h.3,r0.1, Ka15,σ0.05,θ0.2};(*Here, we define a non-linear saturating relationship between transmission and virulence that represents thetransmission-virulence trade-off of our model.*)α[α1_]:=α1;β[α1_]:=c*α1/(c*α1+h);(* System of equations *)dSfull= b-(β[x1]*S*Inf)-(γ*S*(P/(P+1)))-(d*S);dIfull =(β[x1]*S*Inf)+(γ*S*(P/(P+1)))-(d+ρ+α[x1])*Inf;dPfull =(θ*Inf)-σ P+r P (1-P/Ka);dS= dSfull//.Pars;dI =dIfull//.Pars;dP =dPfull//.Pars;(* Equilibrium solution *)Equ=Solve[{dS0, dI0, dP0},{S, Inf, P}];(* Endemic equilibrium. There are multiple solutions and the 5th solution here gives real and positive values *)Equ1=Equ[[5]] ;EqS=S//.Equ1 ;EqI=Inf//.Equ1 ;EqP=P//.Equ1 ;(*Checking stability condition for endemic equilibrium*)Jac=D[{dS, dI, dP}, {{S, Inf, P}}];JacEq=Jac/.Equ1;JacEq/.Pars/.x11//Eigenvalues(* Here we define the next generation matrix. F represents the 'transmission' matrix and minV represents the negative oftransition matrix. The next generation is the product of F and the inverse of minV. Please see supplementary file for derivation of the NGMs *)F= {{β[xm]*S,γ*S*(1/(P+1))},{θ,r(1-P/Ka)}};minV={{d+ρ+α[xm],0},{0,σ}};NGM=F.(Inverse[minV]); (* Invasion fitness = maximum eigenvalue of the NGM *)EVNGM=Eigenvalues[NGM] /.Pars //.S->EqS//.Inf->EqI//.PEqP;(* To generate the Pairwise Invasibility Plots (PIPs), we take the maximum of the eigenvalues from EVNGM. This gives the invasion fitness.The contourplot is then generated that separates the trait space where the invasion fitness is greater than 1 (white region)vs where the invasion fitnessis less than 1 (dark region). Regions where the invasion fitness is greater than 1, the mutant can invade and replace the resident. *)ContourPlot Evaluate[ Re[EVNGM[[2]]], {x1,0,2},{xm,0,2} ], Contours{1}, ContourShading{Black,White}, FrameLabel{"Resident Virulence","Mutant Virulence"}, LabelStyleDirective[Black,15], FrameStyleDirective[Black,15], FrameTicks{{{0.5,1,1.5,2},None},{{0.5,1,1.5,2},None}}, Epilog {Brown,PointSize[Large],Point[{0.169,0.168}]}, Inset Plot x(x+.3),{x,0,2}, PlotStyle{Black}, Frame{True,True,False,False}, FrameLabel{"α","β"}, LabelStyleDirective[Black,15], FrameTicks{{0.5,1,1.5,2.0},{0.2,0.4,0.6,0.8,1.0}} , {1.85,0.7},{1.85,0.7}   

Out[]=

{-0.745087-0.732386,-0.745087+0.732386,-0.0779302-2.17973×

-18

10

}

Out[]=

Fig 4. PIPs for transmission-virulence + shedding-virulence trade-offs

Here we generate the pairwise - invasibility - plot (PIP) to show the combined effect of transmission - virulence and shedding - virulence trade - offs . The procedures are similar to figure 3.

Clear["Global`*"](* Parameters *)Pars={b1,γ 0.35,d 0.2,ρ 0, v 0.02, u3,g1,f1,c1,h.4,r0.1, Ka15,σ0.05};(* Trade-offs. Here we define the shedding-virulence trade-off in addition to the transmission-virulence tradeoffs shown in 3.*)α[α1_]:=α1;β[α1_]:=c*α1/(c*α1+h);θ[α1_]:=g*α1/(g*α1+f); (* System of equations *)dSfull= b-(β[x1]*S*Inf)-(γ*S*(P/(P+1)))-(d*S);dIfull =(β[x1]*S*Inf)+(γ*S*(P/(P+1)))-(d+ρ+α[x1])*Inf;dPfull =(θ[x1]*Inf)-σ P+r P (1-P/Ka);dS= dSfull//.Pars;dI =dIfull//.Pars;dP =dPfull//.Pars;(* Solving for equilibrium solution *)Equ=Solve[{dS0, dI0, dP0},{S, Inf, P}];(* Picking endemic equilibrium *)Equ1=Equ[[5]] ;EqS=S//.Equ1 ;EqI=Inf//.Equ1 ;EqP=P//.Equ1 ;(*Checking stability condition for endemic equilibrium*)Jac=D[{dS, dI, dP}, {{S, Inf, P}}];JacEq=Jac/.Equ1;JacEq/.Pars/.x11//Eigenvalues(*NGM and invasion fitness *)F= {{β[xm]*S,γ*S*(1/(P+1))},{θ[xm],r(1-P/Ka)}};minV={{d+ρ+α[xm],0},{0,σ}};Inverse[minV];NGM=F.(Inverse[minV]); (* Invasion fitness = Maximum eigenvalue of the NGM *)EVNGM=Eigenvalues[NGM] /.Pars //.S->EqS//.Inf->EqI//.PEqP;ContourPlot Evaluate[ Re[EVNGM[[2]]], {x1,0,2},{xm,0,2} ], Contours{1}, ContourShading{Black,White}, FrameLabel{"Resident Virulence","Mutant Virulence"}, LabelStyleDirective[Black,15], FrameStyleDirective[Black,15], FrameTicks{{{0.5,1,1.5,2},None},{{0.5,1.0,1.5,2},None}}, Epilog Brown,PointSize[0.025],Point[{0.32,0.32}], Inset Show Plot x(x+.4),{x,0,2}, PlotStyleBlack, BackgroundWhite , Plot y(y+1),{y,0,2}, PlotStyleBlue, PlotRangeFull, Frame{True,True,False,False}, FrameTicks{{0.5,1,1.5,2.0},{0.4,0.8}}, FrameLabel{"α"}, LabelStyleDirective[Black,12] , {1.93,1},{1.93,1} , Inset[ Style[β,Directive[Black,12]], {1.5,0.87} ], Inset[ Style[θ,Directive[Blue,12]], {1.5,0.72} ]  

Out[]=

{-0.741691+0.71007,-0.741691-0.71007,-0.106593+0.}

Out[]=

Fig 5a. Persistence-virulence trade-off

Here, we look at the effect of persistence-virulence trade-off in virulence outcomes. The code is slightly modified from Figure 4 for faster results. For figure b, use

σ(α)=

a

-uα

e

+v

with a=0.3, u=3, v=1. All other parameters are the same as in 5a.

Clear["Global`*"]Pars={b2,γ 0.15,d 0.2,ρ 0, v 0.07, u3.5,g1,f1,c1,h.4,r.3, Ka15,β.12,θ0.3};α[α1_]:=α1;(*System of equation*)dSfull= b-(β*S*Inf)-(γ*S*(P/(P+1)))-(d*S);dIfull =(β*S*Inf)+(γ*S*(P/(P+1)))-(d+ρ+α[x1])*Inf;dPfull =(θ*Inf)-σ [x1]*P+r P (1-P/Ka);dS= dSfull//.Pars;dI =dIfull//.Pars;dP =dPfull//.Pars;(*Solving for equilibrium solutions *)Equ=Solve[{dS0, dI0, dP0},{S, Inf, P}];(* Passing trade-off function. Although it doesn't matter where we pass this trade-off function in this specific problem, but it matters in multiple trade-off functions. Passing it after solving for equilibrium soulion makes it faster in computation. *)σ[α1_]:=Exp[(-u*α1)]+v;(* Picking endemic equilibrium *)Equ1=Equ[[5]] ;EqS=S//.Equ1 ;EqI=Inf//.Equ1 ;EqP=P//.Equ1 ;(*Stability of endemic equilibrium. If the real part of the maximum eigenvalue is negative, then the equilibrium is stable. *)Jac=D[{dS, dI, dP}, {{S, Inf, P}}];JacEq=Jac/.Equ1;JacEq/.Pars/.x11//Eigenvalues;(*NGM and Invasion fitness *)F= {{β*S,γ*S*(1/(P+1))},{θ,r(1-P/Ka)}};minV={{d+ρ+α[xm],0},{0,σ[xm]}};Inverse[minV];NGM=F.(Inverse[minV]); EVNGM=Eigenvalues[NGM] /.Pars //.S->EqS//.Inf->EqI//.PEqP;ContourPlot[ Evaluate[ Re[ EVNGM[[2]] ], {x1,0,2},{xm,0,2} ], Contours{1}, ContourShading{Black,White}, FrameLabel{"Resident Virulence","Mutant Virulence"}, LabelStyleDirective[Black,15], FrameStyleDirective[Black,15], FrameTicks{{{0.5,1.0,1.5,2.0},None},{{0.5,1.0,1.5,2.0},None}}, Epilog {Red,PointSize[0.025],Point[{0.34,0.34}]}, {Brown,PointSize[0.025],Point[{1.64,1.64}]}, Inset[ Plot[ (Exp[-3*x]+0.02), {x,0,2}, Frame{True,True,False,False}, FrameLabel{"α","σ"}, FrameTicks{{0.5,1,1.5,2},{0.2,0.4,0.6,0.8,1.0}}, PlotStyleOrange, FrameStyleDirective[Black,12], BackgroundWhite, ImageSize160 ], {1.5,0.45} ]  ]

Fig 6a. Bifurcation diagram under multiple trade-offs with changing u (persistence-virulence strength)

In this section, we generate the so called bifurcation diagrams that represent different types of singular strategies. Singular strategies are those strategy values where the fitness gradient vanishes. The techniques of finding different singular strategies are given as comments with the code below. Here, we look at the effect of changing persistence-virulence trade-off by changing the parameter u.

Note: This section is computationally intensive and generates large outputs (> 1 GB).

Fig 6b-d. PIPs for multiple trade-offs

Here, we use persistence-virulence trade-off in addition to the shedding-virulence and transmission-virulence trade-offs. Change the parameters (e.g. u for changing persistence-virulence trade-off) as mentioned in the main text to get different PIPs.

Clear["Global`*"]Pars={b1,γ 0.35,d 0.2,ρ 0, v 0.02, u4,g1,f1,c1,h0.1,r0.07, Ka15};(* Here, 2 trade-offs are defined before solving for equilibria. Third is passed after the solution to make the code run faster. You can change the position and see its effect. *)α[α1_]:=α1;β[α1_]:=c*α1/(c*α1+h);θ[α1_]:=g*α1/(g*α1+f); (* Define system of equations *)dSfull= b-(β[x1]*S*Inf)-(γ*S*(P/(P+1)))-(d*S);dIfull =(β[x1]*S*Inf)+(γ*S*(P/(P+1)))-(d+ρ+α[x1])*Inf;dPfull =(θ[x1]*Inf)-σ[x1] P+r P (1-P/Ka);dS= dSfull//.Pars;dI =dIfull//.Pars;dP =dPfull//.Pars;(* Solve for equilibrium solution *)Equ=Solve[{dS0, dI0, dP0},{S, Inf, P}];(* Define the third trade-off *)σ[α1_]:=Exp[(-1*u*α1)]+0.02 /.Pars;(* Select endemic equilibrium solution *)Equ1=Equ[[5]] ;EqS=S//.Equ1 ;EqI=Inf//.Equ1 ;EqP=P//.Equ1 ;(* Check for stability of endemic equilibrium *)Jac=D[{dS, dI, dP}, {{S, Inf, P}}];JacEq=Jac/.Equ1;JacEq/.Pars/.x11//Eigenvalues(* Define the next-generation matrix *)F= {{β[xm]*S,γ*S*(1/(P+1))},{θ[xm],r(1-P/Ka)}};minV={{d+ρ+α[xm],0},{0,σ[xm]}};Inverse[minV];NGM=F.(Inverse[minV]); (* Find the eigenvalues. Maximum eigenvalue is the desired invasion fitness *)EVNGM=Eigenvalues[NGM] /.Pars //.S->EqS//.Inf->EqI//.PEqP;(* Use the maximum eigenvalue to generate the PIP *)ContourPlot[ Evaluate[ EVNGM[[2]] ], {x1,0,2},{xm,0,2}, Contours{1}, Method{MaxBend0.01}, FrameLabel{"Resident Virulence","Mutant virulence"}, LabelStyleDirective[Black,15], ContourShading{Black,White}, FrameStyleDirective[Black,15], FrameTicks{{{0.5,1.0,1.5,2.0},None},{{0.5,1.0,1.5,2.0},None}}, Epilog {Brown,PointSize[0.025],Point[{0.18,0.18}]}, {Red,PointSize[0.025],Point[{0.42,0.42}]}, {Brown,PointSize[0.025],Point[{1.21,1.21}]}, Inset[ Style["u = 4",Directive[Black,15],FormatTypeTraditionalForm], {0.3,1.8} ]  ]

Out[]=

{-0.779923+0.783034,-0.779923-0.783034,-0.0861077+0.}

Out[]=

Fig 7a. Bifurcation diagram with varying reservoir carrying capacity

Similar to the section 6a, we now generate the bifurcation diagram where we look at the effect of reservoir carrying capacity on different singular strategies.

Note: This section is computationally intensive and generates large outputs (> 1 GB).

Clear["Global`*"]Pars={b1,γ0.35,d0.2,ρ0,v0.02,g1,f1,c1,h.1,r0.1,u3};β[α_]:=c*α/(c*α+h);θ[α_]:=g*α/(g*α+f);dSf=b-β[α]*S*Inf-γ*S*P/(P+1)-d*S /.Pars;dIf=β[α]*S*Inf+γ*S*P/(1+P)-(d+α)*Inf /.Pars;dPf=θ[α]*Inf-σ[α]*P+r*P*(1-P/Ka) /.Pars;Equ=Solve[{dSf0,dIf0,dPf0},{S,Inf,P}];σ[α_]:=Exp[(-u*α)]+v;Eq1=Equ[[5]];E

Article abstract

Introduction

Figure 2a. Differences in threshold parameters: R0d vs λ

Figure 2b. Effect of free-living pathogens and reservoir carrying capacity in endemic disease prevalence

Fig 3. PIPs under transmission-virulence trade-off

Fig 4. PIPs for transmission-virulence + shedding-virulence trade-offs

Fig 5a. Persistence-virulence trade-off

Fig 6a. Bifurcation diagram under multiple trade-offs with changing u (persistence-virulence strength)

Fig 6b-d. PIPs for multiple trade-offs

Fig 7a. Bifurcation diagram with varying reservoir carrying capacity

Figure 2a. Differences in threshold parameters:
R
0d
vs λ