Conditions for Deriving $R_0$ for SIR Model Using Survival Function Method I'm taking a look at the SIR model given by the system of differential equations 
\begin{align} \frac{dS}{dt} & = -  \beta S I \\ \frac{dI}{dt} & =  \beta S I - \gamma I \\ \frac{dR}{dt}& = \gamma I \end{align}
where $S$, $I$, and $R$ represent the number of susceptible, infectious, and removed individuals in a population, and we have that $S + I + R = 1$.
I am trying to justify the derivation of $R_0$, the basic reproductive number, as being equal to $\frac{\beta}{\gamma}$, using its definition as the expected number of secondary infections caused by a single index case in a population where all other individuals are susceptible. 
To do this, I am trying to apply the survival function method outlined on page two of this article http://mysite.science.uottawa.ca/rsmith43/R0Review.pdf. This method, as far as I understand it, computes the value of $R_0$ as $$R_0 = \displaystyle\int_0^{\infty} F(\tau) b(\tau) d\tau$$ where $F(\tau)$ is the probability that the initially infectious individual is still infectious at time $\tau$, and $b(\tau)$ is the rate at which this infectious individual infects susceptibles at time $\tau$. 
My first guess for an approach is based on the fact that the initially infectious individual will deterministically remain infectious for $\frac{1}{\gamma}$ time in the SIR model. I would take $$F(\tau)\ = \left\{
     \begin{array}{ll}
       1 & : \tau <= \frac{1}{\gamma}\\
       0 & : \tau > \frac{1}{\gamma}
     \end{array}
   \right.
$$
And then to interpret $b(\tau) = \beta$ for all values of $\tau$, which will give us $$R_0 = \displaystyle\int_0^{\infty} \beta F(\tau) d\tau = \displaystyle\int_0^{\frac{1}{\gamma}} \beta d\tau = \frac{\beta}{\gamma}$$
This produces the desired result, but I feel that one step might be unjustified. When I assume that $b(\tau) = \beta$, I feel like I am implicitly assuming that, for the entire infectious period of the index case, the proportion of individuals infected by that index case have almost no impact on the availability of susceptibles to be infected.
I guess I feel that I would like conditions on $\beta$, $\gamma$, and the population size such that the primary index case does not infect a substantial enough fraction of the population to slow down its own infection rate. For example, if $\beta$ was very large, it could be possible for the index case to have already depleted a large portion of the susceptibles before exiting its own infectious period.
Do you guys know of any such conditions? Additionally, is it overkill to assume that such an issue could be possible in a somewhat realistic epidemic process on a large population?
 A: Here are a few perspectives; they reflect my opinion, so they may be coloured slightly by my exposure to mathematical biology.
1) $R_0$ is thought of as an expected number of secondary infections caused by a typical infectious individual during its infectious lifetime.  This definition makes no reference to the size of the total population, and being a non-mathematical definition, we can take the mathematical "version" of its definition to be one in which we assume the population is "sufficiently large", so that your depletion issue can never feasibly happen.
2) The survival function method is more of a heuristic construction than a mathematical one.  As an example, the SIR model does not come equipped with a distribution of the duration of infection.  Really, the function $F(\tau)$ should be strictly decreasing, satisfy $F(0)=1$, and have a very long tail, since although we may have never observed a human that is sick with a cold be sick forever, this is not strictly impossible.  However, what you have chosen is a good approximation if the tail falls off very quickly after $\tau=\frac{1}{\gamma}$, which it very well might, if "most" people clear their infection after this time.
3) The survival function method shares a lot in common with the mathematically more rigorous next-generation method.  This defines $R_0$ as the expected number of secondary infections caused by a small quantity of infected individual introduced into a population at equilibrium.  For the finite-dimensional SIR model, this can be computed as follows: 
let $F=\beta$ be the infection rate, and suppose $V_0$ individuals are initially infected.  If we track only this subpopulation of individuals, then their population size will asymptotically fall off, since they will eventually clear their infection.  Let $V(t,V_0)$ denote their population at time $t\geq 0$.  For the SIR model, it can be shown rigorously that this is equal to $V(t,v_0)=e^{-\gamma t}V_0$.
Now, let's count how many new infections these people are responsible for.  The infection rate is $\beta S(t)I(t)$, and there are $V(t,V_0)$ people infected at time $t$.  If we assume that $S(t)\approx S^*$ (this is justfied; see the papers I cite later.), where $S^*$ is the equilibrium susceptible-only population size, then the rate of new infections at time $t$ is approximately equal to $\beta S^*V(t,V_0)=\beta e^{-\gamma t}V_0$.  Integrating this function from $0$ to $\infty$ counts the total number of new infections that the $V_0$ initial infectious people are responsible for.  This is $$\int_0^\infty \beta e^{-\gamma t}V_0 dt = \frac{\beta}{\gamma}V_0.$$  If we average over an acceptable range of initial infectious populations, $V_0$, the result is that $V_0$ is cancelled, and we obtain $R_0=\frac{\beta}{\gamma}$.
Notice that this definition has a lot in common with the integral appearing in the survival function method.  It is the integral of an infection rate multiplied by a quantity that tracks the infection status of the population; the only difference is that this time, it isn't an ill-defined probability, but rather something that comes from the differential equations only.
Summary: The survival function method is not a mathematically rigorous construction; it's something that makes a lot of sense to epidemiologists, but for something simple like the SIR model that doesn't have a lot of extra probabilistic bells and whistles thrown in, it is pretty ill-defined (even though the intuition might give the "correct" answer).  The next-generation method is mathematically rigorous, and for simple models, might very well give the same "answer" as the survival function method does.
For a more thorough exposition on the next-generation method as it applies to ODE models, read the paper by P. van den Driessche and James Watmough (2005).  
