at a fixed per-capita rate, γ, independent of their group. As a result of recovery, individuals are imbued with life-long immunity to future infection. Individuals can also switch groups. 69 For now we use τ ij , φ ij and η ij to represent the per-capita rates at which susceptible, infective, 70 and recovered individuals, repsectively, switch from group j to i. We expand on the details 71 surrounding group switching later. We can summarize the description above using a system 72 of differential equations. Scaling time so that one time unit is equivalent to the average 73 lifetime of an individual in the population, and matching birth and death rates, we get
Note that the differential equations in (1) sum to zero, and so total population size N is 75 constant. This, along with the fact that group membership among recovered individuals is 76 of no consequence, allows us to omit (1e) and (1f). We now use u i = S i /N and v i = I i /N to 77 denote the fraction of susceptible and infective individuals, respectively, in group i. Similarly, 78 we use u = u 0 + u 1 and v = v 0 + v 1 to denote the total fraction of susceptible and infective or may not hold. The linear stability analysis presented in Appendix A shows that the 127 well-known endemic equilibrium remains stable as long as
where R 0 = β 00 /(1 + γ) is the basic reproductive number [9] We consider the evolution of the level of pathogen exploitation of its host, denoted ξ > 0.
Exploitation affects disease transmission, with a greater ξ value corresponding to a greater 137 β ij . To reflect this, we now write β ij (ξ) where
In words, we are treating β ij as an increasing function of ξ that saturates at a value of [27] by assuming the trade-off faced by the pathogen involves recovery, e.g., through increased 152 viral load making it more likely that the pathogen is detected by the host's immune system. 153 We use an adaptive dynamics approach to model the evolution of pathogen exploitation 154 under the primary influence of natural selection [8, 28, 26] . We introduce a rare mutant 155 pathogen with exploitation trait ξ m into a resident pathogen population with exploitation 156 trait ξ. It is assumed that the resident system has reached equilibrium prior to introducing
As we do with the resident population, we can also track the proportion of individuals 162 infected with the mutant strain engaging in prophylaxis. Denoting this proportion by y m , 163 we can describe its dynamics by
If the mutant strain becomes common, the approximate dynamics in (8) break down and we 165 say that the mutant has successfully invaded the resident population. Provided the system 166 is sufficiently close to an evolutionarily steady state, a mutant who successfully invades will 167 become the new resident [29] .
Since v m does not appear in (8b), we can first solve for the equilibrium valueȳ m of 169 y m , substitute that value into (8a), and study (8a) alone. If the right-hand side of (8a) is 170 positive (resp. negative), the mutant invades (resp. is eliminated) because it is favoured 171 (resp. disfavoured) by natural selection. The sign of the right-hand side of (8a) is the same 172 as the sign of the difference between
and unity, where β avgū new infections are created during this time. If this quantity is larger than one (resp. 179 smaller than one), then the mutant population will grow (resp. shrink). Writing W as we 180 have done in equation (9) also highlights the two trade-offs faced by the pathogen. One is 181 the transmission-recovery trade-off described above, captured through the β avg and γ terms.
By unpacking the β avg term, though, we can also see a second trade-off between the type of 183 infection created. Through theȳ m terms, the pathogen is able to influence whether a new 184 infection occurs in a protector or a non-protector. We can, therefore, use W as a payoff 185 function in the adaptive dynamics analysis.
Following the discussion above, the direction of evolution of ξ that is favoured by natural 
An equilibrium valueξ, i.e., one that satisfies condition (10) There are two special cases that can be analyzed with relative ease. The first special case 
The mutant strain is then able to invade (resp. is eliminated) if W (ξ m , ξ) exceeds (resp. is 206 less than) unity; equivalently, if R 0 (ξ m ) exceeds or is less than R 0 (ξ). More importantly, 207 the CSS level of exploitation,ξ, will maximize R 0 (ξ) and soξ = √ κ. This will serve as a 208 benchmark against which more general results will be compared.
The second special case assumes that prophylactic measures are cost-free, i.e., χ = 0. to the pathogen losing access to the trade-off in the type of infection created.
In general, the model cannot be explored analytically. However, if we are near the critical 219 cost χ c we can derive quasi-analytic results. When the cost χ is slightly below its critical 220 threshold χ c , we can approximate the CSS exploitation level asξ
If σ is positive (resp. negative), thenξ is above (resp. below) the benchmark value of √ κ. As 223 equation (14) shows, whether we are above or below this benchmark depends on how small 224 changes in cost lead to small changes in the proportion of susceptible protectors which, in 225 turn, lead to changes in the selection gradient acting on exploitation. 226 We can show with a quasi-analytic approach that equation (14) is always negative. Our 227 evidence relies on first choosing feasible values of our parameters. In particular, we need
To ensure that ε < 1, we then need to choose k > 1 R 0 −1 .
Using feasible parameters and working to zeroth order in χ − χ c , we use the computer 231 algebra software (CAS) Maple (version 2019.1) to investigate the sign of σ as described in 232 equation (14). We find that the requirement that R 0 > 1 necessarily restricts our choices of , and the region of parameter space where χ c > 0 and R 0 > 1. Dotted grey curves represent the roots of the right-hand side of equation (14) and panel B shows that these all occur on or below the black and dashed grey curves (note that vertical dotted grey lines are an artifact of jump discontinuities). Choosing parameter values in the region above the curves in panel A results in the right-hand side of equation (14) being negative, as noted in panel B. The implication is that the pathogen exploitation level (virulence) will decrease relative to the benchmark level of √ κ close to the cost threshold χ c .
Our results can be extended numerically for costs that are possibly much smaller than the 242 critical value, χ c , using a Matlab (version R2019a) procedure detailed in Appendix D. We 243 build the procedure around the observation that locally asymptotically stable equilibrium 244 solutions to dξ/dt = (∂W/∂ξ m )| ξm=ξ are also convergence-stable evolutionary equilibria as 245 defined by conditions (10) and (11), respectively. As a result, numerical iteration of this 246 differential equation can be used to find candidate CSS strategies. The evolutionary stability 247 of candidate CSS strategies can be confirmed with a centred finite-difference approximation 248 of (12). Since the error is on the order of the square of the distance between ξ values used 249 in the approximation, we consider any value within this error to satisfy the ESS condition.
The results of our numerical procedure confirm that the benchmark CSS level ofξ = √ κ 251 is obtained when χ = 0 and χ = χ c . Second, numerical results confirm the reduction in 252 the CSS level of pathogen exploitation for χ slightly smaller than χ c . Third, and most 253 important, numerical results indicate that the CSS exploitation rateξ changes in a simple 254 way as cost is reduced from its critical value to its natural lower limit at zero (figure 2). In 
choose five values of k (k = 1/(R 0 − 1) + 5, k = 1/(R 0 − 1) + 10, k = 1/(R 0 − 1) + 50, 264 k = 1/(R 0 − 1) + 100, and k = 1/(R 0 − 1) + 500) and five values of ε spread evenly between 265 the threshold R 0 /((k +1)(R 0 −1)) and 1, for a total of 500 combinations of parameter values. pathogen with a virulence ξ > √ κ. When the cost is above its critical value so that no one 298 is engaging in prophylaxis, then this mutant cannot invade a resident population at the CSS 299ξ = √ κ. If we then decrease the cost below its critical value so that individuals begin to 300 take prophylactic measures, the CSS exploitation level decreases away fromξ = √ κ and so 301 the mutant is still unable to invade the resident population.
The discrepancy between our prediction and those of Pharaon and Bauch is due to the Moreover, the pathogen's virulence is always below the level predicted in the absence of in-370 dividuals engaging in prophylactic behaviour. There is a balance, then, between the level of 371 cost that is optimal for minimizing prevalence and the level optimal for minimizing pathogen 372 virulence. It is important to recognize that efforts to minimize the cost of prophylaxis will 373 result in a pathogen strain that is more virulent than it otherwise might be. More broadly, 
which, when evaluated at the disease-free equilibrium (DFE) (ū,v) = (1, 0), has eigenvalues 487 λ 1 = −1 and λ 2 = β 00 − 1 − γ. The DFE is stable when both eigenvalues are negative, which 488 leads us to the condition that R 0 = β 00 /(1 + γ) < 1 for stability.
When R 0 exceeds this threshold, the DFE becomes unstable and the standard two-490 dimensional system moves towards an endemic equilibrium (ū,v) = (1/R 0 , (1 − 1/R 0 )/(1 + 491 γ)). To determine the region in which this equilibrium remains stable after we incorporate the 492 host behavioural dynamics, we need to consider the Jacobian matrix of our four-dimensional
where χ c = R 0 − k+1 k β 10 β 00 (R 0 − 1) + 1 and asterisks denote entries that are possibly non- (ū,v,x,ȳ), which we know has negative eigenvalues whenever R 0 > 1. The 2 × 2 matrix in 499 the bottom right is lower triangular, so its eigenvalues are the entries on the main diagonal.
The second of these entries is always negative, while the first is negative as long as χ > χ c .
This defines a critical cost threshold where for χ > χ c the endemic equilibrium (ū,v,x,ȳ) 502 is stable, while for χ < χ c our system tends towards an endemic equilibrium that contains 503 protectors in some non-zero quantities. tion G(r, s; χ). We know that below the critical cost threshold χ c the endemic equilibrium 509 (r,ŝ) = (û,v,x,ŷ,ξ) is stable, while above this threshold our system tends towards the 510 protector-free equilibrium (r,s) = (ū,v, 0, 0,ξ). The critical cost level represents a bifur-511 cation point where these two equilibria coincide and undergo an exchange of stability. To 512 study how our system reacts as we decrease the cost away from this threshold, we intro-513 duce a perturbation parameter δ = χ − χ c and take a first-order approximation to our new 514 equilibrium point (r,ŝ) = (ū + ρ 1 δ,v + ρ 2 δ, ρ 3 δ, ρ 4 δ,ξ + σδ). Knowing that this equilibrium 515 point must satisfy F(r,ŝ; χ) = 0 and G(r,ŝ; χ) = 0, our goal is to find expressions for the 516 perturbation coefficients ρ 1 , ρ 2 , ρ 3 , ρ 4 , and σ.
If we treat (r,ŝ) as a function of χ, we can make a first-order Taylor series approximation 
where subscripts denote partial derivatives and dr = (ρ 1 , ρ 2 , ρ 3 , ρ 4 ) and dŝ = σ are the 520 derivatives with respect to χ ofr andŝ, respectively. We know that (r,s) is an equilibrium 521 point, so the first term in equations (17a) and (17b) evaluates to zero. We also observe that 522 every term in F and G involving χ is multiplied by at least one of x or y, and so the partial 523 derivatives with respect to χ vanish when we evaluate at (r,s). This simplifies (17) 
with asterisks denoting entries that are possibly non-zero. Since J is a block triangular 527 matrix, the eigenvalues are given by the eigenvalues of the matrices on the main diagonal.
The 2 × 2 matrix in the upper left is the Jacobian matrix arising from the linearization of the 529 standard SIR model around the protector-free endemic equilibrium. Since the lower-right 530 3 × 3 block is lower triangular and has a zero entry on its main diagonal, we can see that zero 531 is an eigenvalue of J. This allows us to interpret [ dr dŝ ] as the eigenvector of J associated with 532 the zero eigenvalue, and so shows that there is a non-trivial solution for our perturbation 533 coefficients.
While an analytic expression for this eigenvector can be found, it is unwieldy. Of more 535 interest is the sign of the perturbation coefficient σ, as this tells us in which directionξ 536 moves as we decrease the cost below its critical value. The third row of (19) tells us thatx 537 is a free variable, and the last row tells us that there is a simple relationship between this 538 free variable andξ. In particular, if we consider finding the eigenvector [ dr dŝ ] by solving the 539 expression J [ dr dŝ ] = 0, then the last row of (19) tells us that
We know that the denominator of (20) is always negative sinceξ = √ κ is convergence stable 541 (see (11)), and we know that ρ 3 > 0 since the proportion of susceptible protectors increases as 542 the cost is decreased below its critical value. It follows, then, that the sign of σ is controlled 543 only by the numerator of (20) and so we arrive at equation (14) in the main text. Since the first of these quantities is negative and the second is positive, this shows that the 566 first root is the stable equilibrium: We now define the differential equations for u, v, x, and y, as well as the partial derivative Using this, we put together the matrix J described in Appendix B: If we also choose ε values so that χ c > 0 and k values so that ε < 1, we can generate a series If we use the ε and k values chosen above, we can also look at the value of the numerator of 663 σ as described in equation (14). This gives a second set of β max -κ expressions that we can This produces the plot shown in panel B of figure 1 and shows that all of these curves that 683 represent (14) We also need functions to define the resident system and the partial derivative with respect 718 to ξ m of the fitness function: 
( kappa + xi )*( epsil^2* xbar -1)* cost + ubar^2* xi^2* bmax^2*...
( epsil^2* xbar -2* xbar * epsil + 1)^2)*(1 + xi )^2*( kappa ...
+ xi )^2))*( -((( ubar *( epsil^2* xbar -2* xbar * epsil + 1)* bmax ...
-k * cost )* xi^2 -2* cost * k * kappa * xi -( ubar *( epsil^2* xbar ...
-2* xbar * epsil + 1)* bmax + k * kappa * cost )* kappa )*...
sqrt ( k^2*( kappa + xi )^2* cost^2 -2* k * ubar * xi * bmax *( kappa ...
+ xi )*( epsil^2* xbar -1)* cost + ubar^2* xi^2* bmax^2*...
( epsil^2* xbar -2* xbar * epsil + 1)^2) + ( ubar^2* bmax^2*...
( epsil^2* xbar -2* xbar * epsil + 1)^2 -2* k * ubar * cost *...
( epsil^2* xbar -1)* bmax + cost^2* k^2)* xi^3 + 3* k * cost *...
(( -epsil^2* xbar + 1)* ubar * bmax + k * cost )* kappa * xi^2 ...
+ 3*( -(1/3)* ubar^2* bmax^2*( epsil^2* xbar -2* xbar * epsil ...
+ 1)^2 -(1/3)* k * ubar * cost *( kappa -1)*( epsil^2* xbar ...
-1)* bmax + cost^2* k^2* kappa )* kappa * xi + k *( ubar * bmax *...
( epsil^2* xbar -1) + k * kappa * cost )* cost * kappa^2));
Finally, we use all of these to define a function that numerically approximates the CSS level 