Abstract
Individual-level variation in susceptibility to infection and transmissibility of infection can affect population-level dynamics in epidemic outbreaks. Prior work has incorporated independent variation in susceptibility or transmissibility of individuals into epidemic compartmental models. Here, we develop and assess a mathematical framework that includes covariation in susceptibility and transmissibility. We show that uncorrelated variation in susceptibility and transmissibility leads to an effective transmissibility distribution that has a constant coefficient of variation such that the epidemic dynamics match those with variation in susceptibility alone, providing a baseline for comparison across different correlation structures. Increasing the correlation between susceptibility and transmissibility increases both the speed and strength of the outbreak – and is indicative of outbreaks which might be strongly structured by contact rate variation. In contrast, negative correlations between susceptibility and transmissibility lead to overall weaker outbreaks – with the caveat that the strength of effective transmission increases over time. In either case, correlations can shift the transmissibility distribution, thereby modifying the speed of the epidemic as the susceptible population is depleted. Overall, this work demonstrates how (often unaccounted) covariation in susceptibility and transmission can shape the course of outbreaks and final outbreak sizes.
Highlights
Developed models incorporating susceptibility and transmissibility covariation.
Identified eigendistributions of the force of infection.
Uncorrelated transmissibility reduces to variation in susceptibility alone.
Positive correlations lead to increases in the basic reproduction number.
Positive correlations give faster, stronger, more likely outbreaks.
Effective transmission rates increase over time with negative correlations.
1 Introduction
Individuals differ in response to infection: some people may be more likely to get sick than others, and some people may be more likely to transmit an infection on to others. Variation in susceptibility to infection has been introduced into susceptible-infectious-recovered (SIR) (and related) compartmental epidemic models (Kermack & McKendrick, 1927) to account for intrinsic heterogeneity (Rose et al., 2021; Gomes et al., 2022), extrinsic differences based on, e.g., age-dependent contact rates (Davies et al., 2020; Lovell-Read et al., 2022), or differences in prior immunity (Gart, 1972). Prior research has found that variation in susceptibility reduces the epidemic burden (i.e., outbreak size) relative to the homogeneous model (Ball, 1985; Dushoff & Levin, 1995; Coutinho et al., 1999; Dushoff, 1999; Dwyer et al., 2000; Novozhilov, 2008; Novozhilov, 2012; Karev & Novozhilov, 2019; Britton et al., 2020; Rose et al., 2021; Gomes et al., 2022; Tuschhoff & Kennedy, 2024). Hence, the distribution of heterogeneity in susceptibility and the epidemic burden can jointly vary: when susceptible individuals become infected, the joint variation leads to the redistribution and ‘sculpting’ of the susceptibility distribution. The sculpting leads to epidemic slowdowns relative to that of the homogeneous case, reflecting a fundamental difference in the nonlinearity of incident infections. As shown in Rose et al., 2021, the susceptibility distribution is sculpted toward eigendistributions e.g., including gamma distributions with a constant coefficient of variation. Therefore, outbreaks may appear similar during the early stages, but heterogeneity in susceptibility can slow the speed of the epidemic, leading to lower final outbreak sizes (Rose et al., 2021; Gomes et al., 2022), akin to similar impacts of awareness-induced behaviour change (Eksin et al., 2019).
Variability in transmission has also been studied in epidemic models. In practice, the number of secondary infections caused by a focal infected individual can exhibit significant variability, linked to pathogen type, host features, environmental context, and mode of transmission (Lloyd-Smith et al., 2005; Wong & Collins, 2020; Meehan et al., 2023; Murayama et al., 2023). The extent of variation is often measured in terms of an effective ‘dispersion’ parameter of a negative binomial distribution fit to secondary cases (Lloyd-Smith et al., 2005; Lloyd-Smith, 2007; Endo et al., 2020). During the COVID-19 pandemic, household surveys suggest significant variation in both susceptibility and transmissibility (Anderson et al., 2023); and contact survey data has shown that age-dependent variation in contact rates can be a key driver of variation in transmission (Zhang et al., 2020; Quilty et al., 2024) and susceptibility (Britton et al., 2020). Notably, the presence of heterogeneity in transmission may make diseases harder to control at the outset (Frieden & Lee, 2020; Goyal et al., 2022) but can make them more vulnerable to control via mitigation aimed at reducing the relatively small fraction of superspreading events (Frieden & Lee, 2020; Sneppen et al., 2021). Variation both in susceptibility and transmissibility has previously been introduced to epidemic compartmental models. For instance, in a parasite-host system, susceptibility values were fit to dose-response data, revealing that the transmission rate is lower with more heterogeneity (Dwyer et al., 1997). Additionally, models with uncorrelated variation in susceptibility and transmissibility have previously been explored – showing that power-law distributions in transmission can arise from epidemic dynamics unfolding given initial gamma distributed susceptibility and transmissibility (Novozhilov, 2012).
Together, variation in both susceptibility and transmissibility can shape disease dynamics, but compartmental epidemic models have not yet accounted for the effects of potential covariation on dynamics. For example, individuals who interact more with others may be more likely to become infected and more likely to infect others (e.g., if they continue to interact at similar rates when infectious). In this way contact rates could be considered as being equivalent to correlated variation in susceptibility and transmission – but, we caution that many other factors may also be in play. Moreover, individuals may be more vulnerable to infection due to health, behavioural, social, and/or genetic factors that mean they have limited interactions when infectious (e.g., with trained health-care providers) and are therefore less likely to infect others. Here, we advance an approach that is agnostic with respect to contact variation in order to consider both the sign and strength of potential correlations between suceptibility and transmissibility more generally. In this manuscript, we provide a mathematical framework that includes individual-level variation in both susceptibility and transmissibility. The framework allows for comparisons between different model implementations of variation in susceptibility and transmissibility and makes explicit the consequences of covariation between susceptibility and transmissibility on outbreak dynamics.
2 Model Framework
2.1 Epidemiological dynamics of models with susceptibility and transmissibility
We extend the model framework developed in Rose et al., 2021 to incorporate heterogeneity in both susceptibility and transmissibility into SIR-like epidemic models. To do so, we consider the following population compartment states: susceptible (S), infected (I), and recovered (R). We assume that each individual in the population has a fixed intrinsic susceptibility value, ε, as well as fixed intrinsic transmissibility value, δ. Hence, the S-I-R compartments are functions of susceptibility (ε) and transmissibility (δ) that we denote as S(t, ε, δ), I(t, ε, δ), and R(t, ε, δ). We denote S(t), I(t), R(t) to represent the respective population densities of all susceptible, infected and recovered individuals. Then we can define sub-population densities: the susceptible sub-population density with intrinsic susceptibility ε and intrinsic transmissibility δ is given by
the infected sub-population density is
and the recovered sub-population density is
Note that these definitions are joint densities such that at any time, t, then 
To calculate the mean susceptibility and transmissibility over time, we consider the marginal distributions relevant to the disease dynamics:
The other relevant marginal distribution is the potential transmissibility distribution in the susceptible population. This distribution is indicative of the remaining infectivity of the population who might have the potential to be infectious in the future. During the exponential growth phase of the epidemic individuals are drawn from the susceptible population at varying rates that depend on an individual’s susceptibility to infection. As individuals become infected, their contribution to effective transmissibility is drawn from the potential transmissibility distribution, and thus, the effective transmissibility distribution (Equation 5) is “filled in” by the potential transmissibility distribution given by:
Given the susceptibility distribution (Equation 4) and transmissibility distributions (Equations 5-6), we can define the mean susceptibility and transmissibility. The mean susceptibility to infection within the susceptible population is:
the mean effective transmissibility of individuals within the infected population is:
and the mean potential transmissibility of individuals in the susceptible population is:
Then the force of infection to the susceptible population with susceptibility level ε is:
where β is the baseline transmission rate. For each subpopulation with (ε, δ), the SIR model equations with susceptibility and transmissibility heterogeneity can be written as:
where γ is the recovery rate of all infected individuals. See Appendix A for the derivation of Equation 11 from discrete model variables. See Figure 1 for a visual representation of the model framework using discrete susceptibility and transmissibility variables.
(A) Positive correlation, no correlation, and negative correlation between individual susceptibility and transmissibility. (B) S-I-R compartments are discretized into subpopulations distributed according to susceptibility (εi) and transmissibility (δj). Here, λi(t) = λ(t, εi) is the force of infection to the susceptible population with susceptibility, εi, and γ is the mean recovery rate for all infected individuals.
Integrating with respect to the continuous variables, ε and δ, we can obtain the total population incidence:
From here, differential equations for the following variables can be identified: the mean susceptibility 
the mean potential transmissibility 
and the mean effective transmissibility 
. This model can be described by the following 6-dimensional system of ordinary differential equations (where, for convenience, we set CS(t) = covS(ε, δ)(t) to represent the covariance between susceptibility and transmissibility in the susceptible subpopulation):
Here, the dependence on t is implicit for all time-dependent variables except for the covariance between susceptibility and transmissibility in the susceptible subpopulation, CS(t), and the variance, 
, of the susceptibility distribution, gS(t, ε). This system is not closed for arbitrary starting joint distributions in the susceptible population, as CS(t) and 
may change over time, impacting the mean susceptibility and transmissibility as the epidemic progresses.
2.2 Associations of correlations with epidemic strength and dispersion
Using the dynamical system presented in Equation 13, we can define the basic reproduction number (ℛ0) and dispersion (κ) of epidemics as a function of the correlation coefficient (ρ) between susceptibility and transmissibility. First, in a fully susceptible population, S = 1, the basic reproduction number is given by:
where 
is the mean effective transmissibility during initial exponential growth, and 
is the initial mean susceptibility. Note that without correlations (i.e., ρ = 0) between initial susceptibility and transmissibility values, then 
, but with correlations (i.e., ρ ≠ 0), they can differ. During exponential growth, the joint distribution of the susceptible population is fixed, which implies that the derivatives, 
. Hence, both ε and 
remain fixed during the exponential growth phase. Since the covariance of the joint distribution of the susceptible population, CS(t), also remains fixed during the exponential growth phase, then 
is constant. Hence, in order to obtain the reproduction number, we need to find δI (0) relative to δS(0). In order for 
during exponential growth, then 
must converge to 
, and the reproduction number is
Then the reproduction number can be written equivalently as:
where the correlation coefficient between susceptibility and transmissibility (during exponential growth) is given by 
.
Second, to define the dispersion, κ(t), of an epidemic we first let the mean susceptibility, 
, serve as a dimensionless progress variable. Rose et al., 2021 used 
ds as a dimensionless progress variable, noting that ϕ(t) is proportional to the cumulative infectious force. We compute the mean susceptibility 
, potential transmissibility 
, and effective transmissibility 
as the epidemic progresses over time. We also compute the variances 
and 
in order to obtain the corresponding squared coefficients of variation over time.
From Equation 13, we follow the analysis in Rose et al., 2021 and have that 
. Using differential notation, d(ln(x)) = dx/x, then the square of the coefficient of variation for the susceptibility distribution (CV 2) is given by:
For given initial distributions, we compute the squared coefficient of variation for the susceptibility distribution as well as the squared coefficients of variation for the potential and effective transmissibility distributions over time.
3 Results
3.1 Uncorrelated gamma-distributed susceptibility and transmissibility
We first examine the dynamics of Equation 13 when susceptibility and transmissibility are uncorrelated. In this case, the covariance, CS, is zero. Hence, from Equation 13, 
and so 
remains constant over time. If the initial potential and effective transmissibility values are equal, i.e., 
, then 
. Since 
is constant, then 
remains constant and equal to the initial potential value, 
. Therefore, without correlations, Equation 13 further simplifies to
This nonlinear dynamical system is equivalent to prior work on variation in susceptibility alone (Rose et al., 2021; Gomes et al., 2022). Here, the variance in susceptibility 
denoted with explicit time t to emphasize that it may change over the course of the outbreak. In this work, we also introduce a variation of Equation 18, which we term the ‘reduced model’, in which we approximate the variance term 
using analytical results from (Rose et al., 2021). For uncorrelated gamma distributions, the reduced model sets the variance 
where kε is the shape parameter of the gamma distributed susceptibility distribution gS(ε, t). For uncorrelated Gaussian distributions, the reduced model sets the variance 
.
To examine how the initial distributions changes through the epidemic dynamics, we first consider initially gamma-distributed susceptibility and transmissibility values and examine the squared coefficient of variation in susceptibility, 
(see Equation 17). For initially uncorrelated gamma distributions, we find that:
, where k is the shape parameter of the susceptibility distribution, gS(0, ε). Here, the mean transmissibility, 
, is a multiplicative factor that modifies the force of infection. However, since 
is constant, the system collapses to the system with variation in susceptibility alone. Hence, initially gamma-distributed susceptibility distributions remain gamma-distributed with mean: 
(Rose et al., 2021; Gomes et al., 2022). That is, the same power-law relationship between susceptibility (ε) and the susceptible population (S(t)) holds here with uncorrelated gamma-distributed variation in susceptibility and transmissibility.
We can analyze the change in the joint distribution of susceptibility and transmissibility in the susceptible population, fS(t, ε, δ), through the epidemic dynamics. We find that fS(t, ε, δ) satisfies the partial differential equation:
see Appendix B.1. Note that we can integrate Equation 19 over all transmissibility values so that the marginal susceptibility distribution gS(t, ε) satisfies:
It has been shown (see Section S3 in Rose et al., 2021) that distributions of the exponential family, including initially uncorrelated gamma distributions, with shape parameter k, of the form:
satisfy the PDE given in Equation 20. We show that gamma distributions that evolve with a fixed shape parameter and Gaussian distributions that evolve with a fixed variance satisfy Equation 20 (see Appendix B.2).
To verify this analysis, we can compare the simulations of the discrete model given in Equation 11 with the uncorrelated reduced model in Equation 18, which sets the variance 
Details of model parameterization and simulation are given in Appendix C-E.
For initially uncorrelated gamma-distributed susceptibility and transmissibility, the dynamics of incident infections (Figure 2A, blue) agree with Equation 18 (Figure 2A, green dashed). They also agree with the case of variation in susceptibility alone (Figure 2A, dashed black). Consistent with results of Rose et al., 2021, variation in susceptibility slows down incident infections compared to the conventional SIR model (Figure 2A, gray). As predicted, the κ for both susceptibility transmissibility remain constant over time (Figure 2B,C). We show the initial joint distribution, fS(0, ε, δ), in the susceptible population with uncorrelated susceptibility (ε; x-axis) and potential transmissibility (δ; y-axis) (Figure 2D). In Figure 2E, we compare the susceptibility distribution, gS(t, ε), at two time points during exponential growth: t0 = 0 days (black circle) and t1 = 10 days (violet circle). The distribution remains constant during exponential growth when susceptible depletion is negligible. In Figure 2F, we show the potential and effective transmissibility distributions at these time points. The epidemic is initialized with a few infected individuals, and the effective transmissibility distribution, hI (t, δ), in the infected population is determined by potential transmissibility, hS(t, δ), in the susceptible population. Thus, to remove transients, we make the initial transmissibility distributions equal: hI (0, δ) = hS(0, δ). Then, we show that hI (0, δ) = hI (t1, δ) remains fixed during exponential growth (Figure 2F, gray matches dashed violet). These results indicate that initially uncorrelated gamma distributions remain gamma-distributed such that the mean susceptibility satisfies 
, and the effective mean transmissibility 
remains constant such that the transmissibility distribution hI (t, δ) is constant over the course of the epidemic.
(A) Incident infections. (B) The coefficient of variation (squared) for susceptibility remains constant over time for gamma-distributed susceptibility values. (C) The coefficient of variation (squared) in transmissibility remains constant over time for gamma-distributed transmissibility values. (D) Initial joint distribution (uncorrelated gamma distributions) for susceptibility values (ε) and potential transmissibility values (δ). (E) Comparing susceptibility distributions at time points: t0 = 0 and t1 = 50 days. (F) Potential and effective transmissibility distributions at time points: t0 = 0 and t1 = 50 days. The transmission rate is equal to β = 0.2, and the recovery rate is equal to γ = 1/10 such that the basic reproduction number is ℛ0 = 2.0. Initial gamma distribution shape parameters: kε = 3, kδ = 10. The reduced model refers to Equation 18 with 
.
3.2 Uncorrelated Gaussian-distributed susceptibility and transmissibility
As a prelude to introducing correlations we consider initially uncorrelated truncated Gaussian distributions with low (Figure S1) and high (Figure 3) variance, such that ε, δ > 0. For the low variance case we show incident infections from discrete model simulations (using Equation 11) in Figure S1A. The coefficient of variation (squared) for susceptibility is constant during the initial exponential growth phase of the outbreak, but increases during susceptible depletion because the mean decreases faster than the variance (see Figure S1B). Here, we set σε(t) = σε(0) in simulations of the reduced model (using Equation 18) and observe increases in κ(t) (Figure S1B), meaning that the mean susceptibility decreases. As predicted by our analysis, in reducing Equation 13 in the absence of covariation to Equation 18, the coefficient of variation (squared) for transmissibility remains constant over time (Figure S1C). Despite the differences in κ of susceptibility between the full and reduced models, the reduced model can still approximate incident infections. Consistent with results from Rose et al., 2021, Gaussian and gamma distributions are eigendistributions with respect to the epidemic dynamics. Moreover, gamma distributions have constant κ, whereas Gaussian distributions have constant variance (approximately, considering that Gaussian distributions have proper support on the whole real line).
Similar to Figure S1 but with higher variances in the initial susceptibility and transmissibility distribution. (A) Incident infections. (B) κ for susceptibility. (C) κ for transmissibility. (D) Initial joint distribution of susceptibility and transmissibility in the susceptible population, fS(t0, ε, δ), is given by uncorrelated Gaussian distribution, truncated to have positive support in both ε and δ. (E) Susceptibility distributions at t0 and t1. (F) Potential transmissibility distribution at t0 matches the Effective transmissibility distribution during exponential growth at t1. Parameters: The transmission rate is β = 0.2, and the recovery rate is γ = 1/10 such that the basic reproduction number is ℛ0 = 2.0. The variance values in the initial joint:
, The reduced model refers to Equation 18 with 
.
To better see the effects of correlations when incorporated, we also examine a high variance case in which we increase the variance for both the initial potential susceptibility and transmissibility distributions 
from 0.15 to 0.5; and 
from 0.05 to 0.35). In doing so, the effects of truncation on the bivariate distribution become more pronounced. Even when the initial variances are increased, we can see reasonable agreement, albeit less than with smaller variance, between the full and reduced model simulations (Figure 3A). For truncated Gaussian initial distributions, the κ in the reduced model simulations increase more than in the discrete model simulations (Figure 3B). Due to truncation, susceptibility variance decreases over time as well as the mean susceptibility. The reduced model κ diverges as this model is unable to capture this decrease in variance over time, as by definition 
Despite this, the reduced model in Equation 18 provides a reasonable approximation of incident infections in the discrete model simulations with the highest discrepancies observed during the decay phase (Figure 3A).
3.3 Correlations between susceptibility and transmissibility modify the basic reproduction number
Next, we introduce covariation by considering correlations between susceptibility and transmissibility and compare the ensuing epidemic dynamics with respect to the case without correlations. We vary the correlation coefficient (ρ) from negative to positive, in simulations we explore scenarios over the range of values from −0.6 to 0.6, and find that the speed, i.e., the exponential growth rate, increases with increasing correlation (Figure 4A). Recall that the basic reproductive number is dependent on the correlations between susceptibility and transmissibility (Equation 16). In the absence of correlations (ρ = 0), the basic reproduction number is 
, and the product of the initial mean susceptibility and transmissibility multiply the basic reproduction of the conventional SIR model, β/γ. Note that ℛ0 is an increasing function of the correlation coefficient, ρ, which is in agreement with simulations (Figure 4B). For ρ > 0, the more susceptible individuals are infected earlier and are also more transmissible than on average, causing more transmission during exponential growth such that the basic reproduction number is greater than in the uncorrelated case. For ρ < 0, the basic reproduction number is less than in the uncorrelated case because the more susceptible individuals are less transmissible on average (Figure 4B). For a given transmission rate, β, and recovery rate, γ, the initial speed and strength of the epidemic increases with the initial correlation coefficient between susceptibility and transmissibility, leading to larger outbreaks by larger initial correlation coefficients (Figure 4A).
(A) The speed of incident infections is the exponential growth rate which varies with the initial correlation coefficient. Positive correlations between susceptibility and transmissibility result in faster epidemic speeds with increased peak incident infections. (B) The basic reproduction increases with increasing initial correlation coefficient. Comparison of Equation 14 with computed 
, for the five simulations (colored dots corresponding to the scenarios in A) and Equation 16 with (approximately) fixed initial standard deviations of the initial joint susceptibility distribution (dashed line). Across all simulations, the means of the initial joint in susceptibility (i.e., εS(0) = δS(0) = 1), the transmission rate is equal to β = 0.2, and the recovery rate is equal to γ = 1/10 so that when there is no correlation between susceptibility and transmissibility (ρ = 0), then ℛ0 = 2, as expected in the standard SIR framework. Varying ρ = 0.6, ™0.30.3, 0, 0.3, 0.6, the parameter values of initial joint distribution are given by: 
See Figure S2 for the corresponding initial joint distributions.
3.4 Sensitivity of heterogeneous model outcomes to an outbreak index case
The introduction of population variability also raises questions about how the initiation of an outbreak may effect epidemic outcomes – the transmissibility and susceptibility of individuals in the the first chains of infection may have a large effect on how an outbreak takes off (see Goyal et al., 2022). To probe this, we first assess how variation in the initial distribution of the infected population may impact epidemic trajectory and timing; and second, utilize stochastic simulations to additionally assess variation in outbreak occurrence and epidemic trajectories. In Figure 5A, we show two potential points in susceptibility-transmissibility parameter space that an index infection could take, and examine the impact of these initial conditions on epidemic trajectories in Figure 5B. We find that differences in the transmissibility (but, not susceptibility) of the initial infection can impact the timing of the epidemic – essentially translating the epidemic trajectory in the time axis. When more infectious individuals seed an outbreak, the epidemic trajectory emerges earlier than if the initial case is less infectious than average – in which case the epidemic occurs more slowly than expected under the outbreak eigendistribution and baseline SIR models. As the infection has already occurred, the susceptibility of this (small) index infection does not play a role in ongoing transmission or the long-term epidemic trajectory given mass-action kinetics.
(A) Two different choices for the characteristics of the initial index case index with susceptibility ε = 1 and transmissibility values: δ = 0.5 (gray dot) and δ = 2 (black dot), superimposed atop the the initial joint distribution of the uncorrelated Gaussian case. (B) Using the discrete model, we fix the mean susceptibility ε = 1 and show that varying the transmissibility of the index cases chosen in A can shift the timing of epidemic onsets and peaks. (C) Using stochastic simulations with a randomly chosen index case, we show the probability of an outbreak of more than 50 infections occurring for the SIR model, and for models incorporating negative correlations (ρ = −0.60), no correlations (ρ = 0) and positive correlations (ρ = 0.60) between susceptibility and transmissibility. (D) from the same stochastic ensembles as in C, variation in the final outbreak size (given that outbreaks generate more than 50 infections) is shown, outliers are shown as jittered points. Asterisk’s (*) show the results from the corresponding deterministic simulations from Figure 4. Parameters: transmission rate is β = 0.2 and recovery rate is γ = 0.1. Stochastic simulations were initialized in a population of 10,000. Each stochastic ensemble consists of 1,000 trajectories.
However, the susceptibility (and transmissibility) of individuals in the first few chains of infection may be important. To examine this, we adapted our model to include individual transmission events – utilising stochastic Gillespie simulations (see Appendix E) for the baseline SIR, and the heterogeneous cases with ρ = −0.6, ρ = 0, ρ = 0.6 examined in Figure 4. In doing so, we show that incorporating heterogeneity can additionally alter the probability of an outbreak (here defined as more than 50 infections, see Figure 5C). While the SIR model and the uncorrelated (ρ = 0) model both have an initial reproduction number of ℛ0 = 2, they differ in the likelihood of an outbreak occurring. For the SIR model, the outbreak probability (given m initial infections) is expected as: 
(Southall et al., 2023). With one initial index infection (m = 1) this supports an analytic outbreak probability of p = 0.5, in close agreement with the proportion of stochastic SIR model simulations in which an outbreak occurred (0.483); but further from the uncorrelated model (0.412). The correlated models also differ in the observed proportion of outbreaks relative to the expected outbreak probability; 0.339 vs p = 0.371 for the ρ = −0.6 case, and 0.474 vs. p = 0.5887 for the ρ = 0.6 condition. Regardless of correlation, the introduction of heterogeneity lowers the expected outbreak probability relative to the SIR baseline. Utilizing a stochastic framework also allows us to assess variability in epidemic trajectories (see Figure 5D) whose average final sizes for simulations that ran to epidemic burnout (i.e., simulations which ended due to susceptible depletion, rather than fizzling out) are in good agreement with the deterministic simulations. Histograms of final outbreak size and outbreak duration for all epidemic trajectories are shown in Figure S3. Together with Figure 4 we observe that within our framework, positive correlations between susceptibility and transmissibility lead to epidemics that are more likely to occur, and are faster (with shorter duration), stronger (higher incident infections and final outbreak size) and less variable (final outbreak size interquartile range = 0.0112, for ρ = 0.6), while negative correlation outbreaks are on average less likely to occur, have longer duration and lower, but more variable final size (final outbreak size interquartile range = 0.0239, for ρ = −0.6).
3.5 The effects of correlations on epidemic progress
Initial correlations between susceptibility and transmissibility of the population impact the strength (Figure 4) and potential (Figure 5) of outbreaks. However, correlations have consequences throughout an outbreak. Hence, we next explore how correlations between susceptibility and transmissibility impact epidemic trajectories. To do so, we vary the correlation coefficient, ρ, between ρ = −0.6 and ρ = 0.6 and match the exponential growth rate of incident infections by adjusting the transmission rate, β, to ensure an equivalent basic reproduction number (ℛ0 = 2) (see Figure 6A).We compare the epidemic dynamics during susceptible depletion using the progress variable, 
, i.e., the mean susceptibility (Figure 6B). For 
, and 
, scenarios with positive (light blue), negative (dark blue), and no correlations (blue) reach this susceptibility level at a similar rate; but trajectories diverge as these scenarios move toward 
. For negative correlations, we find the effective transmission rate increases over time, whereas for positive correlations, the effective transmission rate decreases over time (Figure 6C). The justification for this can be seen from the temporal evolution of the bivariate susceptibility-transmissibility distributions (shown by the changes in fS with respect to 
in Figure 6D). As the epidemics progress, the most susceptible individuals are more likely to be infected causing a reduction in mean susceptibility and effectively shifting the underlying distributions. The most salient effect here is the reduction of mean susceptibility over time as the most susceptible individuals become infected and then removed. With positive correlations, the most susceptible individuals are also the most transmissible; thus, as the most transmissible individuals are sculpted into the epidemic, the mean effective transmission rate of the remaining population declines. On the other hand, under negative correlations the effective transmission rates increase over time, as the least transmissible individuals are on average sculpted into the epidemic sooner such that the remaining population has higher average transmissibility. In Figure S4 and Figure S5 we examine these effects on population distributions by showing how the marginal susceptibility and transmissibility distributions in these particular models evolve with the mean susceptibility over time: At mean susceptibility, 
, the susceptible population is depleted to about 80% of the population (Figure S4D). The marginal distributions for susceptibility are similar across ρ = − 0.6 to ρ = 0.6, with small deviations which we attribute to differences in truncation associated with the different correlation coefficients (Figure S4E). We find the marginal distributions for transmissibility differ at 
with positive correlation the transmissibility distribution is shifted to the right, whereas with negative correlations the transmissibility distribution is shifted to the left relative to the transmissibility distribution without correlations (Figure S4F). Since more susceptible individuals are infected earlier during the epidemic with positive correlations, the more susceptible individuals are also more transmissible, leading to increases in the initial speed and strength of the outbreak. With negative correlations, individuals who are more susceptible are less transmissible, leading to decreases in the initial speed and strength. We also further examined the susceptibility and transmissibility marginal distributions over time, comparing across simulations using the epidemic progress variable at values: 
(Figure S5A, middle). For different correlation coefficients, the mean susceptibility, 
, decreases at different rates due to differences in the dynamics of the effective transmission rates (Figure S5A, bottom). We show the susceptibility (Figure S5B) and transmissibility distributions (Figure S5C) over the epidemic progress variable (panels going down): Without correlations (ρ = 0), the effective transmission rate remains constant (Figure S5A). For positive correlations (ρ = 0.6), the effective transmission rate decreases over time, whereas for negative correlations (ρ = − 0.6), the effective transmission rate increases over time. Hence, in either case, the transmissibility distributions tend toward the mean transmissibility of 
(Figure S5C), despite opposite tendencies in the effective transmission rate (Figure S5A, bottom panel). We note these differences in the temporal evolution of the marginal transmissibility distributions are offset by the differences in average transmission rate, β (which is 50% larger when ρ = −0.6 than when ρ = 0.6) in order to have the same ℛ0. However, the underlying differences in correlations cause the differences in effective transmission rate to evolve in opposite directions over time, which contribute to larger final outbreak sizes under negative correlations (Figure S4A).
(A) Incident infections. (B) Mean susceptibility decreases over time as the susceptible population is depleted. (C) Without correlations, the effective transmission rate remains constant. Positive correlations cause the transmissibility to decrease over time, whereas negative correlations cause transmissibility to increase over time. (D) bivariate distributions of fS over time for ρ = 0.6 (top row), ρ = 0 (middle row), ρ =− 0.6 (bottom row) at positions marked by 
in B (columns: left 
, middle 
, right 
). Dashed lines denote average susceptibility and transmissibility.
Overall, simulations of the full PDE model (Equation 11) agree with the qualitative analysis (see Supplementary Information). In particular, correlations modify the speed of susceptible depletion such that the epidemic slows down with positive correlations and speeds up with negative correlations between susceptibility and transmissibility. Consistent with results on heterogeneity in susceptibility, the final outbreak sizes are all less than in the conventional SIR model (see Figures 1-6). In this example, the uncorrelated case leads to about 60% of the initial susceptible population becoming infected, whereas the SIR model leads to about 80% of the susceptible population becoming infected (Figure S6A).
4 Discussion
We developed an epidemic model framework incorporating population-level covariation in both individual susceptibility and transmissibility. Our work investigates how susceptibility and transmissibility distributions are “sculpted” over the course of an epidemic, and how correlated variation may affect population-level dynamical outcomes. Consistent with prior findings (Rose et al., 2021; Gomes et al., 2022) initial gamma and Gaussian distributions are eigendistributions of the force of infection such that dynamics given uncorrelated susceptibility and transmissibility are equivalent to those in which average transmissibility is fixed and only susceptibility varies. Moving to exploring covariation, we identified a relation between ℛ0 and the correlation of initial potential transmissibility and susceptibility. Holding transmission and recovery rates constant, we found that when susceptibility and transmissibility are correlated (anticorrelated), then ℛ0 increases (decreases), epidemics initially grow faster (slower) and are more (less) likely to become outbreaks, and infect more (fewer) individuals. However, if instead ℛ0 is kept constant, we find models with covariation share the same initial epidemic speed, but differ in outcome. In order to keep ℛ0 fixed, and introducing negative correlations between susceptibility and transmissibility necessitates increasing the average transmission rate β (or decreasing the average recovery rate γ), leading to larger epidemics and additionally leading to an increase in the effective transmission rate over time, as individuals who are both highly susceptible and less transmissible are infected earlier on, leaving behind a less susceptible, but more transmissible population. Congruently and counterintuitively, when ℛ0 is fixed with positive correlations between susceptibility and transmissibility, epidemics are smaller than for the uncorrelated case and the effective transmission rate decreases over time.
This approach comes with caveats, insofar as we focus on inherent differences in individual susceptibility and transmissibility in a well-mixed population without vital dynamics and without the potential for reinfections. Going beyond inherent differences, recent work has highlighted that individual-level susceptibility and transmissibility can be associated with human behavior via risk-perceptive decision making (Salomon et al., 2021; Stolerman et al., 2023). Importantly, coupling informed human behaviour with disease dynamics can lead to conditions where dynamic changes in susceptibility (via changing behaviors) can explain epidemic peaks, oscillations, and shoulder behaviors (Weitz et al., 2020; Berestycki et al., 2023). The current model does not allow individual susceptibility or transmissibility to change in time, unlike Weitz et al., 2020 (while neglecting heterogeneity) and Berestycki et al., 2023 (while neglecting variability in transmissibility). Incorporating reinfection and vital dynamics might also enrich the observed dynamics and could allow one to probe differences in heritability of epidemiologically relevant life-history traits. While epidemic burnout is expected in well-mixed SIR models, even with vital dynamics (Parsons et al., 2024), population contact structure is also a highly relevant driver of disease dynamics (Keeling & Eames, 2005; Bansal et al., 2007; Funk et al., 2010; Prem et al., 2021). Our framework does not explicitly represent contact heterogeneity, however, variation in contact rates could be incorporated by assuming a positive correlation between susceptibility and transmissibility – with those who interact the most having the greatest potential to both catch and to transmit an infection. The degree to which contact rate variation factors as the strongest determinant in structuring the underlying joint distribution in susceptibility and transmissibility remains elusive and may differ between diseases and across contexts – however, our expectation is that for a naïve and fully susceptible population there are likely positive correlations between susceptibility and transmissibility. Future extensions might consider additional dynamical effects caused by incorporating additional parameter covariation with recovery rates, heterogeneity in vaccination (Saad-Roy et al., 2024), social determinants of health (Manna et al., 2024; Surasinghe et al., 2024), or with explicit population contact structure. Additionally, further investigation of how susceptibility and transmissibility distributions connect to other distributions of interest, such as the secondary attack rate (Anderson et al., 2023), is warranted.
There are also important questions to be addressed related to parameter inference and outbreak control. As we and others have shown, incorporating individual-level variation provides departures from baseline SIR dynamics (Dushoff, 1999; Novozhilov, 2008; Novozhilov, 2012; Karev & Novozhilov, 2019; Britton et al., 2020; Rose et al., 2021; Gomes et al., 2022; Anderson et al., 2023). In early outbreaks ℛ0 is one of the first parameters epidemiologists attempt to infer, yet our framing suggests ℛ0 might be entangled with covariation in susceptibility and transmissibility. For an identified value of ℛ0, we might expect different epidemic trajectories depending on the degree of covariation in the population. On the other hand, if ℛ0 is identified via average estimations of transmission and recovery rates, the degree of co-variation in the population may lead to mischaracterization of ℛ0. Beyond covariation in population-level susceptibility and transmissibility distributions, heterogeneities in population contact structure also factor into structuring transmission chains, which is not captured in our current models that assume well-mixed populations. Indeed, in network contexts ℛ0 is dependent on both the mean and variance of the degree distribution, as well as the correlation between vertex in- and out- degrees (Allard et al., 2023). Utilizing new inference approaches and data types will be required to identify the degree of covariation between relevant disease parameters e.g., (Kuylen et al., 2022; Anderson et al., 2023; Quilty et al., 2024; Tran-Kiem & Bedford, 2024; Tuschhoff & Kennedy, 2024). Beyond inference of ℛ0 as an early indicator of implementing control measures, there may be additional ramifications if susceptibility or transmissibility covary with infection severity. With a public health goal of minimizing severe outcomes across populations, then if severity is correlated with susceptibility and/or anti-correlated with transmissibility then stronger control measures may be required.
In closing, our work shows how covariation in heterogeneous levels of susceptibiliy and transmissibility scales up to population-level epidemics. Identifying dynamical hallmarks of covariation, and quantifying how multi-dimensional (dynamical) covariation drives population dynamics offer important future avenues to explore. In particular, given the relevance of heterogeneity to the liftoff and potential control of epidemics, we anticipate that decomposing the mechanistic basis of (co)variation of susceptibility and transmissibility will have both fundamental and applied relevance.
Data Availability
MATLAB (version 2023b and 2024a) code for the analysis performed in this manuscript is available at https://github.com/Jeremy-D-Harris/SIR_heterogeneity_project and is archived on Zenodo, doi:10.5281/zenodo.13891898.
https://github.com/Jeremy-D-Harris/SIR_heterogeneity_project
Authorship contribution statement
Jeremy D. Harris: Conceptualization, Formal analysis, Investigation, Software, Supervision, Visualization, Writing – original draft, Writing – review & editing. Esther Gallmeier: Conceptualization, Formal analysis, Investigation, Software, Visualization, Writing – review & editing. Jonathan Dushoff: Writing – review & editing. Stephen J. Beckett: Conceptualization, Investigation, Software, Validation, Visualization, Writing – original draft, Writing – review & editing. Joshua S. Weitz: Conceptualization, Funding acquisition, Project administration, Writing – review & editing.
Code availability
MATLAB (version 2023b and 2024a) code for the analysis performed in this manuscript is available at https://github.com/Jeremy-D-Harris/SIR heterogeneity project and is archived on Zenodo (Harris et al., 2024).
A Derivation of susceptibility and transmissibility from discrete model variables
First, we write the susceptible (S) infected (I), and recovered (R) populations in terms of discrete model variables: Si,j, Ii,j, and Ri,j, where i and j are indices for the discrete susceptibility value, εi, and transmissibility value, δj. Then the force of infection for the susceptible population with susceptibility level εi is
where I(t) is the total infected population and 
is the mean of the effective transmissibility distribution given by
We can write the discrete model equations as
where γ is the recovery rate for infected individuals. This discrete model forms the basis of the continuous model equations given in Equation 11 with the connection between discrete and continuous model variables given as follows: Si,j = S(t, εi, δj), Ii,j = I(t, εi, δj), and Ri,j = R(t, εi, δj). We can calculate the the total incidence:
Next, we derive differential equations for 
and 
found in Equation 13 in terms of discrete model variables. For 
, we examine the time derivative of total susceptibility,
Note that 
. On the other hand, 
. Hence,
which leads to
For 
, we examine the time derivative of the total potential transmissibility in the susceptible population,
Here, covS(ε, δ) is the covariance between susceptibility and transmissibility in the susceptible population. On the other hand, 
. Hence,
so that
For 
, we examine the time derivative the total effective transmissibility in the infected population,
On the other hand, 
. Equating sides and simplifying, we obtain
B Derivations associated with Equation 19 and Equation 20
B.1 Derivation of Equation 19
We derive the partial differential equation that describes the evolution of the susceptible sub-population density fS(t, ε, δ) with intrinsic susceptibility ε and intrinsic transmissibility δ. Rearranging the definition of fS(t, ε, δ) from Equation 1 and taking the partial derivative with respect to time, we obtain
Expansion of the left-hand side through the product rule and the use of Equation 11 give us
Recalling the definition of 
described in Equation 13, and expanding λ from Equation 10, we obtain
which after rearrangement yields the partial differential equation in Equation 19.
B.2 Gamma and Gaussian distributions satisfy Equation 20
We check that gamma distributions with fixed shape parameter and Gaussian distributions with constant variance satisfy Equation 20. Consider setting the marginal susceptibility distribution as a gamma distribution:
Fixing the susceptibility ε and the shape parameter k, we first take the partial derivative with respect to time t and apply the product rule. Subsequently, using Equation 11 (i.e., 
and the fact that 
for gamma distributions, this expression simplifies and we recover Equation 20:
Similarly, consider setting the marginal susceptibility distribution as a Gaussian distribution:
Fixing the susceptibility ε and the standard deviation σ (due to constant variance) and taking the partial derivative with respect to time t, and again recalling 
(see Equation 11), we show Gaussian distributions also satisfy Equation 20:
C Initial joint distributions of susceptibility and transmissibility
To incorporate variation in both susceptibility and transmissibility we use bivariate distributions to initialize our epidemiological models. We do so by first creating the initial joint bivariate distribution in S, using a built-in function from the Statistics and Machine Learning Toolbox (MATLAB version 2023b and 2024a) for the specific probability distributions used. Without correlations between susceptibility and transmissibility, the initial joint distributions in fI and fR are set equal to that in fS. When correlations are introduced, the initial conditions of fI and fR are approximated via the attracting eigendistribution during the exponential growth phase of the epidemic from a preliminary simulation. In this work, we consider three types of bivariate distributions: gamma, truncated Gaussian, and negative binomial.
We first consider independent gamma distributions for the initial joint distribution for susceptibility (ε) and transmissibility (δ). However, these independent gamma distributions do not allow for covariation in ε and δ. Hence, we consider the bivariate Gaussian distribution with 2×2 covariance matrix, allowing us to compare the effects of increasing covariation on epidemic dynamics. We increase the initial variance in susceptibility (σε(0)) and transmissibility (σδ(0)) to better see the effects of covariation. In doing so, the bivariate Gaussian distributions are truncated, because their support lies on the whole real plane ℝ 2. We ensure that mean values are set with 
and adjust variances to match the specified correlation coefficient. In practice, truncation of distributions with high variance in susceptibility and/or transmissibility can be hard to match with large (anti)correlations. In our work, we restrict analysis to correlations between −0.6 and 0.6.
D Model parameters
For all figures, except for Figure 4 which shows the reproduction number as a function of correlation coefficient, we set ℛ 0 =2.0, a value representative of several respiratory viruses such as flu and SARS. We assume that the average time to recovery is 10 days and is exponentially distributed so that the recovery rate is given by γ = 1/10. Hence, we set β = 0.2, with the exception of Figure 6 (and related Supplemental figures) where β is adjusted to compensate for the effects of correlations (between susceptibility (ε) and transmissibility (δ)) and therefore, match the effective exponential growth rate of the epidemics across simulations. The model parameters (descriptions, values, and ranges, thereof) are shown in Table D1.
Epidemic parameters and distribution parameters explored in models with individual traits of susceptibility (ε) and transmissibility (δ), with 
. Baseline values refer to those used in the reference SIR model. Ranges indicate that these parameters vary based on initial conditions which depend on the correlation coefficient. Specific parameter values are noted in Figures.
E Simulation methods
All simulations and analysis were performed using MATLAB (version 2023b and 2024a). All simulation code is available at https://github.com/Jeremy-D-Harris/SIR heterogeneity project and archived on Zenodo (Harris et al., 2024).
E.1 Deterministic simulations
To approximate the continuous susceptibility and transmissibility model variables, we use discrete variables composed of 100 uniformly spaced values between 0 and 6, such that the initial joint distributions we consider are seeded onto a uniform mesh of size 100 × 100. In visualization of the initial joint distributions we show only the range ∈ ([0, 3], [0, 3]), which represents ≈ 90% of the population. In all cases, distributions are chosen such that the initial population average susceptibility 
and potential transmissibility 
are equal to 1. Epidemic model simulations were numerically integrated using ode45 in MATLAB (Dormand & Prince, 1980; Shampine & Reichelt, 1997).
To implement initial conditions, we first create the initial joint distributions of ε and δ in the S, I, and R classes (see Appendix C). In all simulations except for in Figure 5, we let the total population be N = 1, as in Rose et al., 2021 and let the total population initial conditions be: S = N, I = 0, R = 0, with a small perturbation in the direction of the eigenvector of the SIR model. (Adjusting the perturbation magnitude translates the dynamics in time.) The initial conditions that are passed into the ode45 function are calculated from Equation 1 – Equation 3.
E.2 Stochastic simulations
To analyze the outbreak potential of epidemics with different underlying susceptibility and transmissibility characteristics we utilize a stochastic simulation approach using the Gillespie algorithm (Gillespie, 1976, 1977). We initialize simulations with a population of 10,000 whose susceptibility and transmissibility values are seeded with probabilities taken from the joint probability distributions used in discrete model simulations to characterize a representative population with explicit individual-level variation. In each stochastic simulation run, one individual, chosen at random, is designated as the index infection. For each of the initial distributions we analyze (SIR, ρ = − 0.6, ρ = 0, ρ = 0.6) we run the stochastic simulation 1,000 times to obtain ensembles of epidemic trajectories; and denote a threshold of 50 infections to represent the occurrence an outbreak.
Supplemental Material
From equations Equation 11 and Equation 13, the joint distribution in the infected population satisfies the partial differential equation:
Integrating over ε, the effective transmissibility distribution satisfies the partial differential equation:
In the case of uncorrelated susceptibility and transmissibility values, 
, which means that the effective transmissibility distribution remains constant and equal to the potential transmissibility distribution in the susceptible population, if they are initially equal. The mean effective transmissibility remains constant, here equal to 
(Figure S5C; medium blue dot). For positive correlations, the mean effective transmissibility is greater than in the case without correlations (Figure S5C; light blue dot), whereas for negative correlations, the mean effective transmissibility is less than in the case without correlations (Figure S5C; dark blue dot). For positive (negative) correlations between susceptibility and transmissibility, initial incident infections are comprised of more (less) transmissible individuals. Thus, for positive correlations, 
such that 
monotonically decreases toward 
. whereas for negative correlations, 
such that 
monotonically increases toward 
.
Supplemental Figures
Population dynamics with variation in susceptibility and transmissibility. (A) Incident infections. (B) Coefficient of Variation (squared) in susceptibility remains constant during exponential growth and increases over time as the susceptible population decreases. (C) When transmissibility and susceptibility are uncorrelated, Coefficient of Variation (squared) in transmissibility remains constant over time for Gaussian distributions. (D) Initial joint distributions (uncorrelated Gaussian distributions) of susceptibility values (ε) and potential transmissibility values (δ). (E) Susceptibility distributions remain constant during exponential growth, shown at two time points: t0 = 0 and t1 = 50 days. (F) Potential and effective transmissibility distributions at the time points: t0 = 0 and t1 = 50 days. The transmission rate is equal to β = 0.2, and the recovery rate is equal to γ = 1/10 such that the basic reproduction number is ℛ0 = 2.0. The variance values in the initial joint: 
. The reduced model refers to Equation 18 with 
.
Corresponding to simulations in Figure 4, where the joint distribution in I arises from the eigendistribution: (A) ρ = −0.6, (B) ρ = −0.3, (C) ρ = 0, (D) ρ = 0.3, (E) ρ = 0.6.
Count histograms from the four ensembles of 1,000 simulations of epidemic trajectories shown in Figure 5C,D. Final outbreak size (A) and duration (days) (B) for the SIR model, where 51.7% of simulations generated ≤ 50 infections. Final outbreak size (C) and duration (days) (D) for the model with negative correlation between susceptibility and potential transmissibility (ρ = − 0.6), where 66.1% of simulations generated ≤ 50 infections. Final outbreak size (E) and duration (days) (F) for the model with no correlation between susceptibility and potential transmissibility (ρ = 0), where 58.8% of simulations generated ≤ 50 infections. Final outbreak size (G) and duration (days) (H) for the model with positive correlation between susceptibility and potential transmissibility (ρ = 0.6), where 52.6% of simulations generated ≤ 50 infections. Parameters: transmission rate is β = 0.2 and recovery rate is γ = 0.1. Stochastic simulations were initialized in a population of 10,000.
(A) Incident infections. (B) Mean susceptibility decreases over time as the susceptible population is depleted. (C) Without correlations, the effective transmission rate remains constant (medium blue). Positive correlations cause the transmissibility to decrease over time (light blue), whereas negative correlations cause transmissibility to increase over time (dark blue). (D) About 80% of the susceptible population is depleted. (E) Susceptibility distributions plotted at 
(F) Transmissibility distribution at time points corresponding to the progress variable, 
For positive (negative) correlations, the mean transmissibility is greater (less) than the mean transmissibility without correlations. Across all simulations, the recovery rate is γ = 1/10 and the basic reproduction number is ℛ 0 = 2.0. The transmission rates vary: β = 0.254, 0.2, 0.165 for negative correlation (ρ = 0.6), no correlation (ρ = 0), and positive correlation (ρ = 0.6), respectively. The initial variance values are: 
and 
.
(A) Population dynamics including incident infections (η), the mean susceptibility 
(both redrawn from Figure S4A,B), and the effective transmission rate 
. (B) Susceptibility distributions associated with four mean susceptibility values going down the rows:
,(C) Effective transmissibility distributions at the same four time points corresponding to the mean susceptibility values. Parameter values are the same as in Figure S4. Corresponding coefficients of variation (squared) are shown in Figure S6.
(A) Cumulative infections showing different final outbreak sizes for different correlations between susceptibility and transmissibility. Parameter values are the same as in Figure S4 and Figure S5. The coefficients of variation (squared) for susceptibility (B) and transmissibility (C).
Acknowledgements
We thank Tapan Goel for constructive feedback and Raunak Dey for code review. This work was supported by grants from the National Science Foundation (no. 2200269) to JSW, the Simons Foundation (no. 930382) to JSW, and support from the Charities in Aid Foundation, the Marier Cunningham Foundation, and the Chaires Blaise Pascal program of the Île-de-France region (Blaise Pascal Institute Chair of Excellence awarded to JSW).
Footnotes
References
