Infection in a village

Question

Consider the following problem:

Suppose a lonely wanderer infected with a virus came into an isolated village with $M$ villagers and stayed there. Every week each of the infected villagers coughs onto $n$ random other villagers (each of them chosen uniformly and independently among everyone) and then develops antibodies becoming immune to it. All villagers who are coughed upon become infected if they are not immune. Nobody left or entered the village after the arrival of the lonely wanderer. Consider time to be discrete and measured in weeks. We say, that the virus survives as long as someone is infected with it. For what $n$ is the expected time of its survival the longest?

The extremum clearly is not achieved in the border cases here.

Indeed, if $n = 0$ the lonely wanderer becomes immune before being able to infect anyone else, thus the virus will survive only for $1$ week.

If $n \to \infty$ the probability that the lonely wanderer infects everyone in the first week tends to $1$. Thus the expected time of the survival of the virus tends to $2$ in this case.

So, we must look for optimal $n$ somewhere in between. However, I have no idea how to find it (or even its asymptotic for large $M$)…

At the first glance the problem looked to be somewhat similar to two well studied problems: branching processes (villagers infected by a given infected villager - their descendants in terms of branching processes) and coupon collector problem (uninfected villagers as coupons to be collected). However, it is different from both of them (the number of ‘descendants’ changes each turn here, which makes it different from a Galton-Watson branching process, and the number of ‘coupons collected per turn’ depends on the number of ‘coupons collected on the previous turn’, which makes it different from a classical coupon collector) and methods, similar to the ones used to solve them, are unlikely to work here.

Answer 1

This answer uses programming to solve the above question, using principles from probability and Markov chains. No closed mathematical formula is provided, although results are shown for $m < 100$ and curve fitting is used to estimate the average survival time for larger numbers.

Initially, I simulated $1000$ runs of different configurations for $3 \lt m \lt 100$ and $0 \leqslant n \leqslant 4$, plotting the results using $95$% confidence intervals:

For the configuration in which $n = 1$, we find that, with $X$ being the week in which the disease last makes a new victim, the following probabilities apply:

$$\begin{cases} P(X = 1) = \frac{1}{m} \\ P(X = x) = \frac{m}{x - 1} \frac{m - x + 1}{m} \frac{x}{m} P(X = x - 1) \end{cases} \iff P(X = x) = \frac{x}{m^x} \frac{(m - 1)!}{(m - x)!}$$

For $n > 1$, Markov chains can be used. For any combination of $m, n$, let us define a state $\{k, l, w\}$ as a tuple containing the number of immune people $k$, the number of infectuous people $l$ and the current week $w$, with $k + l \leqslant m + 1$. From this state, we can transition to any of the states $\{k + l, 0, w + 1\}, \{k + l, 1, w + 1\}, \ldots, \{k + l, \min(n l, m + 1 - k - l), w + 1\}$.

The probability of each possible state transition can be calculated using recursion. The probability that $l$ infectuous people affect $0 \leqslant p \leqslant \min(n l, m + 1 - k - l)$ while $k$ people are immune, can be determined by first considering all possible actions of a single person. This person coughs onto $n$ random people, infecting $i \leqslant \min(n, p)$ people. To do so, this person must randomly choose $i$ healthy people and $n - i$ immune people, which happens with probability:

$$P(i) = \underbrace{\frac{m - k - l + 1}{m} \frac{m - k - l}{m - 1} \cdots \frac{m - k - l - i + 2}{m - i + 1}}_{\text{healthy people}} \\ \underbrace{\frac{k + l - 1}{m - i} \frac{k + l - 2}{m - i - 1} \cdots \frac{k + l - n + i}{m - n + 1}}_{\text{immune people}} {n \choose i}$$

Using this expression, the probability that a total of $p$ people are infected from a given state $\{k, l, w\}$ can be written as:

$$P(\{k + l, p, w + 1\} | \{k, l, w\}) = \sum_{i = 0}^{\min(p, n)} P(i) P(\{k + l, p - i, w + 1\} | \{k + i, l - 1, w\})$$

This formula can be solved recursively until $p = 0$ or $l = 0$. Considering the visited states, it is worth noting that we always start in the state $\{0, 1, 1\}$. Assuming that $m \geqslant n$, this state transitions to the state $\{1, n, 2\}$ with probability $1$. We can then calculate the probability of each successive state, resulting in different Markov chains.

Considering all chains that end in a state in which the number of infectuous people $l$ equals $0$, we can directly calculate the expected time of survival by considering the variable $w$. Using a Python-based implementation to do this for $3 \lt m \lt 100$ and $0 \leqslant n \leqslant 4$ results in:

The results show that the expected surival time is highest for $n = 2$ when $m \leqslant 43$, while it is highest for $n = 1$ when $m \geqslant 44$. Square root curve fitting results in $f(m) \approx 1.25 \sqrt{m} + 0.70$ for $n = 1$, while logarithmic curve fitting results in $f(m) \approx 2.68 \ln(m) - 1.12$ for $n = 2$. Values around the turning point are as follows:

$$\begin{array}{c|c|c} m & n = 1 & n = 2 \\ \hline \text{39} & 8.51 & 8.65 \\ \text{40} & 8.61 & 8.72 \\ \text{41} & 8.71 & 8.79 \\ \text{42} & 8.80 & 8.86 \\ \text{43} & 8.90 & 8.92 \\ \text{44} & 9.00 & 8.99 \\ \text{45} & 9.09 & 9.05 \\ \text{46} & 9.18 & 9.12 \\ \text{47} & 9.27 & 9.18 \\ \text{48} & 9.36 & 9.24 \end{array}$$

The full Python code has been made available here.

Answer 2

Putting that at week $w$ we have the following numbers of infected, susceptible (infectable), and immunized $$ I(w) = {\rm infected}\quad S(w) = {\rm infectable}\quad H(w) = {\rm immunized} $$ then, allowing for the numbers to be rationals, we shall have $$ \left\{ \matrix{ I(0) = 1\quad S(0) = m\quad H(0) = 0 \hfill \cr I(w) + S(w) + H(w) = m + 1 \hfill \cr H(w) = I(w - 1) \hfill \cr I(w+1) = {n \over {m + 1}}S(w) \hfill \cr} \right. $$

Solving for $I$ (excluding the banal case $n=0$) $$ I(w) + {{m + 1} \over n}I(w + 1) + I(w - 1) = m + 1\quad \left| \matrix{ \;I(0) = 1 \hfill \cr \;I(1) = n{m \over {m + 1}} \hfill \cr} \right. $$ and recasting it by including the initial conditions, so that the recursion is valid for any $0 \le w$, under the acception that $w <0 \; \Rightarrow \; I(w)=0 $ $$ I(w) + {n \over {m + 1}}I(w - 1) + {n \over {m + 1}}I(w - 2) + \left[ {0 = w} \right]\left( {n - 1} \right) = n\quad \left| {\;0 \le w} \right. $$ We can solve the above via ordinary generating function $$ F(z) = \sum\limits_{0\, \le \,w} {I(w)\,z^{\,w} } $$ readily obtaining $$ \eqalign{ & F(z) = {{{n \over {1 - z}} - \left( {n - 1} \right)} \over {1 + {n \over {m + 1}}\left( {z + z^{\,2} } \right)}} = {{m + 1} \over n}{{1 + \left( {n - 1} \right)z} \over {\left( {1 - z} \right)\left( {z^{\,2} + z + {{m + 1} \over n}} \right)}} = \cr & = {{m + 1} \over n}{{1 + \left( {n - 1} \right)z} \over {\left( {1 - z} \right)\left( {z + {1 \over 2} + {{\sqrt {1 - 4\left( {m + 1} \right)/n} } \over 2}} \right) \left( {z + {1 \over 2} - {{\sqrt {1 - 4\left( {m + 1} \right)/n} } \over 2}} \right)}} = \cr & = {A \over {\left( {1 - z} \right)}} + {B \over {\left( {z + {1 \over 2} + {{\sqrt {1 - 4\left( {m + 1} \right)/n} } \over 2}} \right)}} + {C \over {\left( {z + {1 \over 2} - {{\sqrt {1 - 4\left( {m + 1} \right)/n} } \over 2}} \right)}} \cr} $$ where, however, $A,B,C$ are rather complicated expressions of $m,n$.

Normally it would be $n \le m$, so that the square roots will turn into immaginary values.