Probability distribution of the number of infected people in a room Suppose we have a room of $N$ people. The number of people who have the common cold is $k$, which is equal to $1$ at the start. Now, $h$ handshakes occur completely randomly. If an infected person shakes hands with a uninfected person, the uninfected person gets infected and can then go on to infect other people if they shake hands with them. Otherwise, nothing happens. What is the probability distribution of $k$ after all the handshakes occur?
My attempts to solve this:
The probability of an additional person being infected with the cold in 1 handshake would be represented by
$\frac{2k(N-k)}{N(N-1)}$
, but aside from this, I don't know where to go next.
I also looked at hypergeometric distributions and compound probability distributions. Would either be related?
 A: Let $\mathbf{K} = \{ K_h | h \in \mathbb{N}_{0+} \}$ denote the stochastic time-series showing the number of infected people after each handshake, and let $K_0 = 1$ at the start of the series.  This is a Markov chain that falls within the category of discrete "pure birth" processes".  A single random handshake gives the transition probabilities:
$$p_{k,k+r} \equiv \mathbb{P}( K_{h+1} = k+r | K_h = k ) = \begin{cases}
1-\frac{2k(N-k)}{N(N-1)} & & \text{if } r=0, \\[6pt]
\frac{2k(N-k)}{N(N-1)} & & \text{if } r=1, \\[8pt]
0 & & \text{otherwise}. \\[6pt]
\end{cases}$$
Thus, the transition matrix for the chain is:
$$\mathbf{P} \equiv \begin{bmatrix}
1-\frac{2}{N} & \frac{2}{N} & \cdots & 0 & 0 & 0 \\ 
0 & \frac{4(N-2)}{N(N-1)} & \cdots & 0 & 0 & 0 \\  
\vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ 
0 & 0 & \cdots & 1-\frac{4(N-2)}{N(N-1)} & \frac{4(N-2)}{N(N-1)} & 0 \\ 
0 & 0 & \cdots & 0 & 1-\frac{2}{N} & \frac{2}{N} \\ 
0 & 0 & \cdots & 0 & 0 & 1 \\ 
\end{bmatrix}.$$
After $h$ random handshakes, the probability that $k$ people are infected is:
$$\mathbb{P}(K_{h} = k) = [\mathbf{P}^{h}]_{1,k}.$$
You can compute this probability by programming the transition probability matrix into an appropriate piece of mathematical software (e.g., R) and then obtaining the first row of the appropriate power of the matrix.  If you would like to try to get a closed-form expression for the probability, I would recommend deriving the eigen-decomposition or Jordan decomposition of the matrix to see if this simplifies the problem.

Computing the probability vector: It is quite simple to program this Markov chain in R.  In the code below I create a function to compute the vector of probabilities (or log-probabilities) for arbitrary input values for the number of people and the number of handshakes. 
#Load required library
library(expm);

#Create a function to compute the probability vector
COMPUTE_PROBS <- function(N, h, log.p = FALSE) {

    #Define the transition probability matrix
    P <- matrix(0, nrow = N, ncol = N);
    for (k in 1:N)     { P[k,k]   <- 1 - 2*k*(N-k)/(N*(N-1)); }
    for (k in 1:(N-1)) { P[k,k+1] <- 1 - P[k,k]; }

    #Compute probability vector
    PPP <- expm::'%^%'(P,h);
    if (log.p) { PPP <- log(PPP); }
    PPP[1, ]; }

We can use this function to compute the vector of probailities (or log-probabilities) for arbitrary values of N and h.  Here is an example using some chosen values for the parameters.
#Compute an example of this probability vector
N <- 40;
h <- 80;
PROBS <- COMPUTE_PROBS(N,h);

#Plot the probability mass function
library(ggplot2);
DATA  <- data.frame(Infected = 1:N, Probability = PROBS);
THEME <- theme(plot.title    = element_text(hjust = 0.5, size = 14, face = 'bold'),
               plot.subtitle = element_text(hjust = 0.5, face = 'bold'));
ggplot(aes(x = Infected, y = Probability), data = DATA) +
    geom_bar(stat = 'identity', fill = 'blue') + THEME +
    ggtitle('PMF of Number of Infected People') +
    labs(subtitle = paste0('(', N, ' people and ', h, ' handshakes)')) +
    xlab('Number of Infected People');



Monte-Carlo simulation: We can confirm that the above result is correct by comparing the theoretical probabilities to Monte-Carlo simulations of the process.  To do this, we can program a simulation function in  R.  (Hat tip to user2974951 for suggesting this approach, and writing the initial code.)  In the code below I create a function to simulate outcomes of the chain and take empirical estimates of the vector of probabilities for arbitrary input values for the number of people and the number of handshakes.
#Create a function to simulate the chain
SIMULATE_CHAIN <- function(N, h, times = 10^5) {

    #Set the simulation vector
    SIM <- rep(0, times);

    #Run simulations
    for (s in 1:times) {

        #Compute initial vector of infected people
        INFECTED <- c(1, rep(0, N-1));

        #Implement random handshakes
        for (i in 1:h) {
            H <- sample(1:N, size = 2, replace = FALSE);
            if (INFECTED[H[1]] == 0 & INFECTED[H[2]] == 1) { INFECTED[H[1]] <- 1 }
            if (INFECTED[H[1]] == 1 & INFECTED[H[2]] == 0) { INFECTED[H[2]] <- 1 } }
        SIM[s] <- sum(INFECTED); }

    SIM; }

Using this function we can simulate the Markov chain and take empirical estimates of the probabilities of each outcome.  The plot confirms the same shape we obtained in our theoretical analysis, which confirms that the calculations are correct.
#Simulate the chain
set.seed(1)
SIMS <- SIMULATE_CHAIN(N,h);

#Estimate the probability vector
PROBS_EST <- rep(0,N);
for (i in 1:N) { PROBS_EST[i] <- sum(SIMS == i)/length(SIMS); } 

#Plot the probability mass function
DATA  <- data.frame(Infected = 1:N, Probability = PROBS_EST);
THEME <- theme(plot.title    = element_text(hjust = 0.5, size = 14, face = 'bold'),
               plot.subtitle = element_text(hjust = 0.5, face = 'bold'));
ggplot(aes(x = Infected, y = Probability), data = DATA) +
    geom_bar(stat = 'identity', fill = 'red') + THEME +
    ggtitle('Monte-Carlo estimate of PMF of Number of Infected People') +
    labs(subtitle = paste0('(', N, ' people and ', h, ' handshakes)')) +
    xlab('Number of Infected People') + ylab('Estimated Probability');


A: Here is an alternative to Ben's answer using simulations in R, using his parameters.
Edit: fixed the bug.
N=40 #number of people
h=80 #handshakes
k=1 #number of infected people at the start
n=1e5 #number of simulations

result=rep(NA,n)
for (r in 1:n) {
  initial=rep(0,N) #N healthy people
  initial[1:k]=1 #k infected
  for (t in 1:h) {
    random2=sample(1:N,2) #two random people
    if (initial[random2[1]]==1 | initial[random2[2]]==1) {
      initial[random2[1]]=initial[random2[2]]=1 #now both infected
    }
  }
  result[r]=sum(initial)
}

which looks like this

