Should the rate term in the basic epidemiology models be second-order?

Global leaders are taking unprecendented steps with long-term financial and societal impact to prevent the spread of COVID-19, relying on early data and epidemiological models of differing complexity. The standard "compartment" models of the spread of disease through a population follow the law of mass action; that is, the rate of conversion of individuals from one compartment to another is directly proportional to the number of individuals in the former. The basic S-I-R model, which tracks susceptible, infectious, and recovered people, is functionally identical to autocatalytic reactions studied in chemical kinetics. However, if we account for heterogeneity in individuals' susceptibility to catch the disease, then the most susceptible individuals are likely to be removed from the pool early; this is analogous to evaporative cooling in chemical kinetics. Here, we use standard tools from statistical physics to account for this effect, which suggest that the rate of appearance of new cases should be proportional to the \emph{square} of the number of susceptible people in the population, not the first power. We validate this finding with both Monte Carlo (stochastic) and many-compartment differential equation (deterministic) simulations. This leads to quite different predictions of the peak infection rate and final outbreak size, with the second-order model giving milder predictions of each. This has dramatic effects on the extrapolation of a disease's progression when the data collected comes from the early, exponential-like region, and may remove a systematic bias from such models.