Summary:

This article aims to mathematically describe a co-evolutionary dynamics of virulence and drug resistance of a pathogen in a population infected by drug-sensitive and drug-resistant strains. The proposed model assumes that there is a correlation between virulence and the rate at which infecteds receive treatment. Co-evolution of the pathogen virulence and drug resistance under the virulence-treatment relationship is investigated. My comments and questions are as follows.

1. I think the introduction section could benefit from adding specific examples of viral, bacterial or fungal infections to which the proposed hypothesis/model, that is, infected individuals with increased virulence receive treatment more quickly and their recovery is thus faster, could be applied.

2. The model 1a-1c assumes that all infecteds receive treatment though at different rates depending on the virulence of the strain type ‘i’. A fraction (1 - rhoi) recovers and a fraction rhoi transitions to the resistant state. Here, if I understand correctly, recovery occurs immediately upon treatment administration, is positively correlated with virulence (more virulent implies more rapid treatment and thus recovery) and if there is no treatment (gamma = 0), there is no recovery. I wonder if having the natural recovery of drug-resistant strain (traded-off with virulence) in the model might be reasonable since people can spontaneously recover from a large number of viral and bacterial infections. Also, the model assumes that all infecteds receive treatment; perhaps having treated and untreated classes would be interesting to incorporate and investigate. For the case of self-medication/prophylaxis, a paper here https://doi.org/10.1111/j.0014-3820.2005.tb01008.x might be worth looking at.

3. The model assumes that drug treatment directly affects recovery. I could imagine that there might be indirect effect of the drug through e.g. decreasing the infectious pathogen load/virulence and thus potentially accelerating recovery. Is that something that could be investigated via the model and affect the co-evolutionary dynamics?

4. In the relationships for virulence and transmission of the resistant strains, alphai^R = a*alpha and betai^R = a*beta, the assumption is that a >1 or 0 <= b < 1, so that drug-resistance results in higher virulence and decreased transmission. It might be useful to elaborate more on these assumptions and the choice of their values.

5. The trade-off between virulence and treatment rate (p. 19, Fig. S3 A) assumes a linear positive relationship between the two. I imagine that this relationship could be sigmoid so that those infected with low to mid-virulent pathogens would receive treatment at slower rates. I would be curious how this would affect the outcome.

6. On Fig. S3 A, do these points represent the values used in the simulations with multiple strains? Perhaps it would be good to have a more detailed description of the used methods also in the text so that one does not need to look into the R code.

7. In the equation S9c, rho is assumed to be a function of virulence. Is this assumption used anywhere else in the manuscript or simulations?

The present manuscript offers a mathematical view on why the resistant and more virulent strain can be selected in epidemics. I have found interesting that the author highlighted the important aspect that this can be in particular attributed to public health policies.

Although I like and appreciate mathematics, I think it would be beneficial to give some ground for the performed simulations: for example, to adjust some values and time scale for a particular pathogen in mind. Otherwise, it is difficult to understand, e.g., what the time units for Fig.1 (is the x-axis in ages like for HIV or TB or in weeks like for influenza?). It would be also great to see more connection of the framework with real examples/diseases outlined in the introduction. Even though the discussion is quite detailed, I wish to see the same for Introduction section. In this respect, there is an imbalance between Introduction and Discussion sections.

Main remarks: L32: I have encountered a small contradiction in the framework. As the author writes about incubation period vs illness onset, but his main model does not include the "E" compartment. How the results would be changed or not changed if one would consider SEIR model instead of SIR model and how important the inclusion of the "E" compartment in view of the evolutionary dynamics? L63: I think that the examples make sense, but, e.g., TB is not modeled by a simple SIR model. In this view, the paragraph is more suited for introduction part. Then the author writes about parasites in L70, but all examples above except of Plasmodium falciparum are not parasites. Hence there is some mixture of terms, and it can be confusing.

Minor: L28: Would it be possible to clarify the term "pleiotropic action" for a general reader? L30: "the latter trait is usually binary rather than continuous." – the author does not explain much on why it is important. At the moment that sentence seems out of place. L32: I would advice to change the word order: the incubation period and the illness onset, because the former comes first. Introduction: I have noticed the gap in the literature dates, and it looks suspicious, because there should be some scientific works in the last few years. For example, Sylvain Gandon probably tackled some relevant problems, but it is just my guess. For example, the work Todd L. Parsons et al 2018 RSIF (doi:10.1098/rsif.2018.0135) looks relevant at least for the notion. L85: I would write "a disease-free equilibrium" instead of "a trivial state".

This paper is an interesting addition to the literature on the interaction between epidemiological processes. The author states and then illustrates than differential treatment of individuals based on the virulence of their infecting pathogen strains can cause changes in the absolute and relative levels of virulence in the resistant and sensitive strains.

The analysis appears robust and the conclusions are interesting. The analysis of the impact of noise in transmission rates is appreciated, but a more general consideration of the robustness of the results to changes in parameter values would add to the manuscript, as if the effect occurs only in a small region of the parameter space (in this, admittedly, highly abstracted model), it is less likely to have large biological relevance in the real world.

On the topic of modelling drug resistance as a continuous trait, the author admits that that the inclusion would render the model unable to show the effect stated (due to modelling the modelling framework). While this is not an inditement of the work in this paper, it might be worth the author’s time to add to this section to describe ways that continuous modelling of virulence could be performed. Just off hand, it appears that this could be done using (for instance) something analogous to the integro-partial differential equation framework of Delgado-Eckert et al in The evolution of virulence in RNA viruses under a competition-colonization trade-off (https://arxiv.org/pdf/1001.3101.pdf), albeit at the cost of the simple mathematical formulation, that is a strength of this paper. Something to consider.

However, overall I enjoyed reading the paper and it seems like a valuable addition to the literature!

General notes:

Line 53 & 54: “by a drug resistant host” - it is the pathogen that is drug resistant, not the host
The section after equations 2a and 2b is confusing. It eventually becomes clear the author is using x-bar as a generic stand in for the weighted mean of any parameter in the equations, but this would be best stated explicitly, so the reader is not attempting to find and work out what trait x and y actually correspond to.
The equilibria are defined have no strain subscripts at line 84, but equilibria 2 and 3 are given with subscripts. Is this meant to imply that each strain persists with its own equilibrium level given by their own strain specific parameters? If so, this would be better explicitly stated.
Equation 4 - as far as I can see, the subscript m is not defined either in the text or the appendix. It first appears after the calculation of the eigenvalues in appendix A, but is not commented on. Its meaning should be explicitly stated.
Line below equation 4 - apparition is presumably mean to be appearance?
Line 89 - I assume the dichotomy being referred to is which of the eigenvalues dominates and therefore determines the fate of the mutant? But this should be clarified.
Line 148 - misspelling of average