Modeling sub-exponential epidemic growth dynamics through unobserved individual heterogeneity: a frailty model approach

1.   Introduction

The scaling of growth during the initial phase of an infectious disease outbreak is influenced by a complex interplay of factors, many of which are not directly observable. These factors include the nature of the contact network over which the disease spreads ^[1,2,3,4,5], variability in individuals' risk of infection, and changes in behavior or the implementation of control measures ^[6]. To place our study in context within the evolution of epidemiological modeling, it is important to understand how these factors contribute to the non-linear dynamics of disease spread. Under the assumption of homogeneous mixing, a scenario where each individual has an equal chance of interacting with others, the incidence curve is expected to grow exponentially during the initial epidemic phase, provided there are no changes in susceptibility, behavior, or interventions ^[7]. This exponential growth presumes a constant rate intrinsic growth rate, a condition rarely met in real-world scenarios. In contrast, empirical data frequently reveal a sub-exponential growth in the early transmission phase, suggesting a more complex and heterogeneous interaction pattern among individuals. This pattern can emerge from spatial limitations within the contact network partly associated with the transmission mode of the disease, early behavioral modifications, or the introduction of public health interventions ^[6]. Moreover, the sub-exponential dynamic is intricately linked to the concept of "generation interval contraction", a phenomenon observed when the average time between successive cases in a chain of transmission shortens as the epidemic progresses ^[8]. This occurs when a susceptible individual is exposed to multiple infectious contacts, creating a competitive scenario for transmission. Such dynamics underscore the critical role of underlying network structures and behavioral responses in shaping the trajectory of an epidemic, leading to the investigation of complex network effects in epidemiological modeling ^{[9,10,11,12,13,14]}.

The crucial task of accurately characterizing epidemic growth is underscored by the varied growth patterns observed in empirical data from both historical and recent outbreaks, as documented in the prior studies ^[15]. These diverse patterns challenge the conventional exponential growth assumption, underscoring the need for developing alternative modeling approaches that reflects the heterogeneity of epidemic processes. The two-parameter generalized-growth model (GGM) offers a nuanced method for estimating the epidemic growth scaling, employing the scaling of growth parameter $0 < p < 1$ ^[15]. The exploration of this parameter space enables a deeper understanding of the underlying mechanisms driving different epidemic trajectories. Recent research further validates the utility of incorporating these nuanced growth dynamics into epidemic models. Studies have shown that models accounting for the possibility of early sub-exponential growth, rather than just assuming a simple exponential increase, significantly enhance the models' fit to empirical data and improve the accuracy of short-term epidemic forecasts ^{[16,17,18,19,20]}. This advancement in modeling aligns more closely with observed epidemic patterns and provides a more reliable foundation for predicting the course of outbreaks, thereby informing more effective public health interventions.

While the generalized-growth model captures the nuances of early sub-exponential growth patterns identified in empirical outbreak data, this model also suggests simple underlying mechanisms for these patterns. However, a comprehensive mathematical framework that integrates the classic mechanistic theories of infectious diseases with the observed variability in early growth dynamics remains elusive. To bridge this gap, our work introduces the concept of frailty from survival analysis into epidemic modeling, aiming to reflect the observed heterogeneity in individual disease susceptibility and transmission. By incorporating individual-level unobserved heterogeneity into established mechanistic transmission models, we aim to provide a more accurate and nuanced understanding of epidemic spread, better aligning theoretical predictions with empirical observations

2.   Model initial epidemic growth in the presence of unobservable heterogeneous infectious contacts through random effect

Understanding the dynamics of epidemic growth is crucial for developing effective public health interventions. Traditional compartmental models often assume homogeneous mixing and uniform susceptibility within a population, leading to predictions of exponential growth during the initial phase of an epidemic. However, empirical observations frequently show sub-exponential growth, indicating that these assumptions may not capture the full complexity of real-world epidemic dynamics. One approach to addressing this discrepancy is the use of phenomenological models, which describe the observed growth patterns without necessarily detailing the underlying transmission mechanisms. To bridge the gap between phenomenological and mechanistic approaches, we introduce a frailty model, which incorporates individual-level heterogeneity into the epidemic growth dynamics. Frailty, a concept from survival analysis, represents unobserved variability in individual susceptibility and transmission potential. This variability can significantly impact the overall epidemic growth rate.

Population growth is described by the instantaneous increase in the total count $C(t)$ at any given time $t$. This is mathematically represented as:

Within the realm of infectious diseases, $C(t)$ signifies the cumulative count of infections at time $t.$ In the context of our model, $r(t)$ represents the instantaneous growth rate of the epidemic at time $t$. Unlike the constant growth rate $r$ in the exponential model, $r(t)$ can vary over time to reflect changes in the population's susceptibility and interaction patterns. This variation is influenced by the distribution of frailty among individuals, leading to a more nuanced understanding of epidemic growth.

2.1.   The growth rate $r(t)$ in phenomenological and in transmission models

The expression (2.1) has been connected with empirical models to describe population-based phenomena without "mechanically" modelling the transmission dynamics. For instance,

$1)$ The exponential growth equation $C^{\prime}(t) = rC(t)$ is commonly used to approximate the early stages of an outbreak, assuming a very large population and during the first few generations of disease transmission.

$2)$ The General Growth Model (GGM) as described in ^[6,15,21,22] is given by the growth equation

where the rate $r(t) = r/C(t)^{1-p}$ is a function that decreases over time $t$. Assuming the outbreak began with just one infected individual, meaning $C(0) = 1,$ there is a solution given by $C(t) = (1+r\theta t)^{1/\theta},$where $\theta = 1-p,$ leading to $r(t) = r/\left(1+r\theta t\right).$ The GGM returns to the exponential growth model when $\theta\rightarrow0.$

$3)$ The logistic growth model is described by the equation:

Here, $K = \lim_{t\rightarrow \infty}C(t)$ represents the maximum possible value, known as the carrying capacity. It is evident that the growth rate, $r(t) = r\left(1-\frac{C(t)}{K}\right)$ also diminishes over time in accordance with the logistic rule. This can be explicitly expressed as:

The exponential growth model is the special case when $K$ approaches infinity.

$4)$ The Richards model is represented by the equation:

Beginning with an initial condition of $C(0) = 1,$ the model can also be solved analytically for $C(t)$, leading to the expression:

In this model, when $K\rightarrow \infty$, the growth rate decreases over time as a negative exponential function $r(t)\rightarrow re^{-r\theta t}.$ It approaches constant value $r$ as $\theta\rightarrow0.$

In these models, the growth rate $r(t)$ phenomenologically incorporates a range of influences from both individual behaviors and environmental conditions that impact how the disease spreads. However, it is noteworthy that distinct "mechanical epidemic models" can yield identical phenomenon in terms of the cumulative growth. For instance, the logistic growth rate above may be derived from various deterministic transmission models. This includes, among others, the SI (Susceptible-Infected) and SIS (Susceptible-Infected-Susceptible) models ^[23].

2.2.   The exponential growth approximation

We begin our discussion by looking into the exponential growth dynamics, which is defined by a steady growth rate, denoted as $r(t) = r$. This approximation of exponential growth is typical for the initial phase of an epidemic in various deterministic and stochastic models, under the assumption that the population size is infinitely large, making the depletion in susceptible individuals within the population negligible. These models further assume that mixing and transmission among the population are uniform. The growth rate, $r$, is then a function of additional parameters that are specific to particular mechanistic models. We will consider the following examples.

In the Simple Epidemic model, as outlined by Bailey in ^[24], also known as the SI model, the count of infectious persons at any specific time $t$ is expressed by the equation $I^{\prime}(t) = \beta\left(1-I(t)/N\right) I(t)$. In this context, $I(t)$ is equivalent to $C(t)$, and $N$ denotes the entire population size. The parameter $\beta$ is the transmission rate per contact between susceptible and infectious individuals. As $N$ tends toward infinity, $C(t)$ converges to $i_{0}e^{-\beta t}.$ Thus, $r = \beta.$

In the SIS and SIR (Susceptible-Infected-Recovered) models, which are framed using ordinary differential equations, the rate of change in the number of infected individuals at any given moment (time $t$) is described by the formula $I^{\prime}(t) = \beta S(t)I(t)/N-\gamma I(t)$. In this context, $N$ represents the total population size, and $\beta$ has the same interpretation as it does in the SI model, indicating the transmission rate. $S(t)$ and $I(t)$ stand for the counts of susceptible and infected individuals at time $t$, respectively. The $\gamma$ parameter represents the rate at which infected individuals recover, moving either back to being susceptible or gaining immunity. If the reduction in susceptible individuals is minor, the fraction $S(t)/N$ can be approximated as $1$. This simplification leads to the equation $I^{\prime}(t) = (\beta-\gamma)I(t).$ If $\beta > \gamma,$ it suggests that the number of infected individuals, $I(t)$, could increase exponentially without limit. Many standard texts denote $r = \beta-\gamma$, although this assumes conditions of time-homogenity, often described as being "under equilibrium conditions". In this paper, we delve deeper into these equilibrium conditions and examine whether it is accurate to model the initial growth of a disease outbreak, represented by $C(t)$ with the equation $C^{\prime}(t) = (\beta-\gamma)C(t),$ from the outbreak's onset.

By introducing a new compartment, $E(t)$, which tracks the number of people who are infected but not yet capable of spreading the infection as of time $t$, the SIS and SIR models are expanded into the SEIS and SEIR models, respectively. These models use ordinary differential equations to calculate how the numbers of latent and infectious individuals evolve over time. Specifically, the equations $E^{\prime}(t) = \beta\frac{S(t)I(t)}{N}-\alpha E(t)$ and $I^{\prime}(t) = \alpha E(t)-\gamma I(t)$ describe the rates of change for latent and infectious individuals at any moment. Here, the new parameter $\alpha$ represents the rate at which latent individuals progress to being infectious. Under the same set of equilibrium conditions as before, there exists an implicit function between $r$ and the parameters $(\beta, \alpha, \gamma)$ through the equation $\left(r+\alpha\right) \left(r+\gamma\right) = \beta\alpha.$ This relationship is also a recognized finding in standard mathematical epidemiology textbooks such as ^[25] and others.

These examples assume homogeneity so that the infection spread rate, denoted as $\beta$, remains the same for all individuals. Additionally, the natural increase rate, $r$, is less than or equal to $\beta$. The equation $\left(r+\alpha\right) \left(r+\gamma\right) = \beta\alpha$ can be re-written as $r = \beta\frac{\alpha}{r+\alpha}\frac{r}{r+\gamma} < \beta$. The exception is the SI model, where there is no latent period and once individuals are infected, they stay infectious indefinitely, which means $\alpha = \gamma = 0$, leading to $r = \beta$.

2.3.   From stochastic and sampling perspectives

From a stochastic viewpoint, considering $C(t)$ as a counting process, the immediate increase of $C(t)$ given that $C(t) = n$ is explained as follows:

In this context, we posit that every one of the $n$ infected persons independently possesses a positive rate of $r_{i}(t) > 0,$ $i = 1, ..., n,$ which contributes to the emergence of a new infection. At the level of the entire population, the rate of immediate growth is represented by $h(t)$, which is $h(t) = \sum_{i = 1}^{n}r_{i}(t).$ Assuming $r_{i}(t) = r(t)$ for all $i = 1, ..., n,$ we arrive at $h(t) = nr(t).$

We further investigate the scenario of uninterrupted growth when the rate of change is constant $r(t) = r.$ This situation is linked to an exponentially distributed random variable $Y$ associated with a typical infected individual, characterized by the survival function $\overline{G}(t) = \Pr(Y > t) = e^{-rt}.$ The random variable $Y$ is well-defined and can be quantified based on the timing for each person who gets infected, thanks to its exponential distribution and the property of being memoryless.

The stochastic counterpart of the exponential growth $C^{\prime}(t) = rC(t)$ is the conditional probability

The marginal distribution of $C(t)$ at time $t$ is negative binomial, with mean value satisfying the exponential growth $E[C(t)] = i_{0}e^{rt},$ where $i_{0} = C(0)$ is the initial condition. Consequently, we arrive at the formulas:

where the first formula can be interpreted as $r$ representing the growth rate of the total number of infected individuals; and the second formula views $r$ as the hazard rate of the exponential distribution, considering it from the viewpoint of an average individual who is infected.

Yan and Chowell ^[23] highlight that, if the growth is exponential with constant rate $r,$ it can be expressed by adjusting the transmission rate $\beta$ based on the probability of individuals within the population being infectious at any particular moment, by considering a specific sample taken at time $t$. A binary indicator, $\Delta_{i}(t)$, is defined to indicate whether an individual $i$, who was infected prior to time $t$, is infectious at time $t$. This is represented by the probability $\Pr(\Delta_{i}(t) = 1)$, where a value of $1$ signifies the individual is infectious, and a value of $0$ indicating otherwise. Those individuals for whom $\Delta_{i}(t) = 1$ are grouped into what is called a "prevalence cohort". Under suitable equilibrium conditions, to be elaborated on in the following subsection, the relationship

is established, indicating that the probability $\Pr(\Delta_{i}(t) = 1)$ stays constant, regardless of the time of sampling, $t$.

The expressions $r = \beta r\left(r+\gamma\right) ^{-1}$ and $r = \alpha\beta r\left(r+\alpha\right) ^{-1}\left(r+\gamma\right) ^{-1}$ in the aforementioned SIS, SIR, SEIS and SEIR models are the special cases of (2.7) where the durations of both the latent and infectious phases are assumed to follow exponential distributions with rates $\alpha$ and $\gamma$, respectively. This underlying assumption is crucial for models that utilize ordinary differential equations.

We begin by defining the random variable $T_{I},$ which symbolizes the duration of infection among infected people. This duration follows a probability distribution with a density function denoted as $f_{I}(x)$, where $x$ measures time starting from when a person becomes infectious. If we ignore any latent period, assuming individuals are infectious right from the infection's start and $Y$ is independent of $T_{I}$, the probability $\Pr(Y > T_{I})$ can be calculated using the integral $\int_{0}^{\infty} e^{-rt}f_{I}(t)dt = L[f_{I}](r),$ where $L[f_{I}](r)$ is the Laplace transform applied to the probability density function $f_{I}(x)$. This leads us to understand that, if all infected individuals have the same constant hazard rate $r$ to produce new infections, then the chance of an individual causing a new infection at time $t$ while still infectious is given by

The probability $\Pr(\Delta_{i}(t) = 1)$ is established through a stochastic comparison between the random variable $Y$ and the infectious period $T_{I}$ of infected individuals. Specifically, in instances where the infectious period follows an exponential distribution with an average $\gamma^{-1}$, the expression $r = \beta\left(1-L[f_{I}](r)\right) = \beta r/\left(r+\gamma\right)$ simplyfies to $r = \beta-\gamma.$

Now we introduce a random variable $T_{E},$ representing the period of latency in individuals who are infected but not yet capable of transmitting the infection. We consider $T_{E}$ and $T_{I}$ to be sequential and independent of each other. Assuming $Y$ is also independent from both $T_{E}$ and $T_{I}$, we start with $L[f_{E}](r) = \Pr(Y > T_{E})$ and $\Pr(T_{E} < Y\leq T_{E}+T_{I}) = \Pr(Y > T_{E})\Pr(T_{E} < Y\leq T_{E}+T_{I}|Y > T_{E}).$ Given that the distribution of $Y$ is exponential and memoryless, $\Pr(T_{E} < Y\leq T_{E}+T_{I}|Y > T_{E}) = \Pr(Y\leq T_{I}) = 1-L[f_{I}](r)$. Thus, $\Pr(\Delta _{i}(t) = 1) = \Pr(T_{E} < Y\leq T_{E}+T_{I})$ and

Specifically, in instances where the both $T_{E}$ and $T_{I}$ ard exponentially distributed with averages $\alpha^{-1}$ and $\gamma^{-1},$ respectively, the formula is expressed as $r = \beta\left(1-L[f_{I} ](r)\right) = \beta r/\left(r+\gamma\right).$

2.4.   Multifaceted aspects of homogeneity

We have assumed that all the $C(t) = n$ infected individuals are homogeneous in the sense that they carry the same rate to produce a new infection independently from the other $n-1$ individuals. This rate $r$ is further factorized in (2.7) as the product of transmission rate $\beta$ and the sampling probability $\Pr(\Delta_{i}(t) = 1).$

2.4.1.   Time-homogeneity and equilibrium:

The constant growth rate implies that $\Pr(\Delta_{i}(t) = 1)$ is independent of time $t.$ This is time-homogeneity. It requires that the dynamic system is at a state of equilibrium.

In the case without the latent period, (2.8) becomes $r = \beta\left(1-L[f_{I}](r)\right)$ and can be re-written as $\beta\left(1-L[f]_{I} (r)\right) /r = 1.$ Integration by parts, one can show that $\left(1-L[f_{I}](r)\right) /r = \int_{0}^{\infty}e^{-rt}\overline{F}_{I}(t)dt,$ where $\overline{F}_{I}(x) = \Pr(T_{I} > x)$ is the survival function of the infectious period. Thus $r = \beta\left(1-L[f_{I}](r)\right)$ can be also written as $\beta\mu_{I}\int_{0}^{\infty}e^{-rt}\overline{F}_{I}(x)/\mu_{i}dt = 1.$ Denote $f_{w}(x) = \overline{F}_{I}(x)/\mu_{i}$ and $R_{0} = \beta\mu_{I},$ we have the equation

The condition $R_{0} > 1$ is implied, because $\int_{0}^{\infty}f_{w}(t)dt = 1$ and $L[f_{W}](r) = \int_{0}^{\infty}e^{-rt}f_{W}(t)dt < 1.$

The Eq (2.9) is recognized as the Euler-Lokta equation in mathematical biology involving two intrinsic parameters, $R_{0}$ and $r.$ $R_{0}$ is widely known as the basic reproduction number. This equation links these parameters by stabilizing the dynamic system in a state of equilibrium. Furthermore, the function $f_{w}(x)$ adheres to the characteristics of a probability density function and is widely acknowledged as the equilibrium distribution according to renewal theory in the context of stochastic processes.

Figure 1 shows how the equilibrium distribution works, using an example where there are $C(t) = 8$ individuals at time $t$, with $6$ of them being infectious. As a result, the probability $\Pr(\Delta_{i}(t) = 1) = 0.75.$ Each individual is represented by a beginning and an end point. For a given individual $i$ where $\Delta_{i}(t) = 1,$ two time-related random variables$,$ $U_{i}$ and $V_{i}$, are associated with them at time $t$. In renewal process, $U_{i}$ refers to the backward recurrence time, which measures the time backward from t to when the individual became infectious. On the other hand, $V_{i}$ represents the forward recurrence time, measuring the time from $t$ until the individual is no longer infectious. For equilibrium to be achieved, it is necessary that $U_{i}$ and $V_{i}$ have the same distribution as the distribution of $W_{i}$, with each $W_{i}$ being independent and identically distributed according to the equilibrium distribution, which has a probability density function denoted by $f_{w}.$ In addition, the distribution of the sum $T_{i}^{(B)} = U_{i}+V_{i}$ is stochastically longer than that of the infectious period $T_{I}.$ This occurs because the distribution of infectious periods for individuals in the prevalence cohort at time $t$ is biased towards longer durations. Individuals with longer infectious periods are more likely to be included in the cohort. The probability density of $T_{i}^{(B)}$ is represented by $xf_{I}(x)/\mu_{I},$ which is referred to as the length-biased distribution.

Thus, under equilibrium, the following pair of relationships are equivalent:

where the second line treats $W$ as the forward recurrent time $V_{i}$ and each of the infected individual at time $t$ has a probability of $1/R_{0}$ failing to produce a secondary infection during the infectious period. In the special case of exponentially distributed infectious periods with $\mu _{I} = \gamma^{-1},$ the equilibrium distribution is identical to the original exponential distribution. Therefore $L[f_{W}](r) = L[f_{I}](r) = \gamma/\left(r+\gamma\right),$ leading to $r = \beta-\gamma = \gamma(R_{0}-1)$.

Without further elaboration, adding a latent period $T_{E}$ not only gives (2.8) but also extends (2.9) as

For further details, we refer to ^[23].

Challenges to the assumption of equilibrium conditions Even though numerous theoretical findings are based on the assumption of equilibrium conditions, it is improbable that these conditions hold true for data from the initial stages of a major epidemic outbreak. Consider the onset of an outbreak, where nearly every person infected is capable of spreading the disease, resembling the dynamics of an SI model. Here, the likelihood of an individual being sampled, $\Pr(\Delta_{i}(t) = 1)$, is close to $1$, and the rate of infection, $r$, is approximately equal to $\beta$. However, as time goes on and some of the infected individuals recover by time $t$, the probability $\Pr(\Delta_{i}(t) = 1)$ drops below $1$. Yet, it remains unrealistic to expect that the distributions of $U_{i}$ and $V_{i}$ align as they would under equilibrium conditions. When the system is not in equilibrium, the sampling probability not only varies with time $t$ but also the likelihood of each infected individual being included in the prevalence cohort changes based on how recently they were infected relative to time $t$. In essence, time-nonhomogenity induces one aspect of individual heterogeneity.

2.4.2.   Homogeneity assumed in the agent-host-environment interface

The $\beta$ parameter, often expressed as $\beta = \lambda p$, where $\lambda$ represents the rate of contact within the social network and $p$ is the probability of transmission during an interaction between an infected and a susceptible person, implies several assumptions about uniformity across the interactions between agents, hosts, and the environment.

Regarding the infectious agent, like a virus in viral infections, it is assumed that there are no mutations that alters its infectiousness during the study, ensuring that infections acquired at various times confer the same level of infectiousness at the point of infection.

Concerning the hosts, it is assumed that:

$1)$ All susceptible individuals have the same level of susceptibility;

$2)$ An infected individuals maintain a consistent level of infectiousness throughout their infectious period;

$3)$ All infected individuals have equal infectiousness.

In terms of the environment, it is assumed that there's uniform mixing, meaning each person is equally likely to interact with every other person in the population. In large populations, this social network can be simplified as a Bernoulli random graph like that described in Erdös and Rényi ^[26]. Viewing from a stochastic processes angle, the social contact network expands in a manner where the typical individual's number of contacts adheres to a stationary Poisson process with rate $\lambda$.

With these homogeneity assumptions, at any given moment during an epidemic, all infected individuals are equally infectious, regardless of when or how long ago they were infected. Each interaction between susceptible and infected individuals carries the same transmission probability $p$. Together, these assumptions of homogeneity in the agent-host-environment dynamic lead to the formulation $\beta = \lambda p$. This equation defines the rate of new infections in the population as $\beta\frac{S(t)I(t)}{N}$ in deterministic models for disease transmission.

2.5.   Frailty as a model to account for non-identifiable and unobservable heterogeneities

The growth rate can be expressed as $r = \lambda p\delta$, where $\delta = \Pr(\Delta_{i}(t) = 1)$. Data gathered at the group level cannot distinguish $\lambda, p$ and $\delta$ as distinct factors. Each factor reflects specific assumptions about homogeneity. We use a frailty model to describe variability among individuals in a phenomenological manner.

In the field of survival analysis, the proportional hazard model, denoted as $\rho(t|z) = z\rho_{0}(t)$, accounts for variability in case individuals exhibit varying hazard rates. Here, $\rho_{0}(t) > 0$ represents a fundamental hazard rate, and $z > 0$ is the frailty variable. This can also be depicted via survival functions, $G(t|z) = e^{-zH_{0}(t)}$, $H_{0}(t) = \int_{0}^{t}\rho _{0}(u)du$.. The widely utilized Cox proportional hazard model formulates $z$ as a log-linear regression based on observed variables.

The frailty model is a type of random effect model designed to account for unobserved heterogeneity among individuals. It posits that the variable $z$ is random, with an expected mean of $E(z) = 1$ and a probability density function represented by $\xi(z)$. In scenarios where there is no variation among individuals, $\xi(z)$ becomes a singular point where $z$ is exactly $1$, leading to no variability. Under these conditions, the survival function is defined as $\overline{G}_{0}(t) = e^{-H_{0}(t)}$. Consequently, when the population includes a diverse mix of individuals, the resulting survival function is derived from a combined distribution:

where $L\mathcal{[\xi]}(s) = \int_{0}^{\infty}e^{-zs}\xi(z)dz$ is the Laplace transform with respect to $\xi(z)$ and $L\mathcal{[\xi]}(H_{0}(t))$ is $L\mathcal{[\xi]}(s)$ evaluated at $s = H_{0}(t).$ The hazard function becomes

We use the frailty framework with a fixed baseline hazard rate, denoted as $\rho_{0}(t) = r$, and cumulative hazard rate $H_{0}(t) = rt$. In the frailty model, the hazard rate given frailty $z$ is $\rho(t|z) = zr$, leading to a mixed hazard rate expression $\rho^{(mixed)}(t) = -\frac{d}{dt}\log L\mathcal{[\xi ]}(rt)$. The following general observations are made:

$1)$ For any non-degenerate probability density function $\xi(z)$, the Laplace transform of $\xi$ at $rt$, satisfies $L\mathcal{[\xi]}(rt) = \int _{0}^{\infty}e^{-zrt}\xi(z)dz\geq e^{-rt}$. Referring to (2.6) when $z$ is constantly $1$, $r = \frac{d}{dt}\log C(t) = -\frac{d}{dt} \log\overline{G}(t),$ leading to the relationship $C(t) = \overline{G}(t)^{-1}$. Introducing the random effect (frailty) $\xi(z)$ and assuming $C(0) = 1$, we have $\frac{d}{dt}\log C(t) = -\frac{d}{dt}\log L\mathcal{[\xi]}(rt)$, resulting in the inequality

indicating that the growth of $C(t)$ is limited by the exponential function $e^{rt}$ from above.

$2)$ It has been demonstrated ^[27] that $L\mathcal{[\xi ]}(rt)$ is log-convex. As a survival function, its hazard rate $\rho (t) = -\frac{d}{dt}\log L\mathcal{[\xi]}(rt)$ decreases monotonically from an initial value of $r$ at time $0$. This reflects a decrease in the effective reproduction number, $R_{t}$, during the initial phase of epidemic growth, as noted by ^[21].

$3)$ When comparing the variability of $z$ using convex order (as detailed in the Appendix), it is found that higher variability in $z$ results in a higher value of the Laplace transform $L\mathcal{[\xi]}(s)$ for all $s > 0$. This, in turn, means that the more variable $z$ is, the higher the value of $L\mathcal{[\xi]}(rt)$ and thus, the smaller the value of $C(t)$ for any $t > 0$.

2.6.   Special cases

2.6.1.   Gamma distributed heterogeneity

Let $\xi(z)$ be the p.d.f. of the Gamma distribution with $E(z) = 1$ and variance $var[z] = v > 0.$ The Laplace transform is $L\mathcal{[\xi]}(s) = \left(1+s\nu\right) ^{-1/v}$. Letting $C(0) = 1,$ the expression $C(t) = \left[ L\mathcal{[\xi]}(rt)\right] ^{-1}$ becomes

Meanwhile the instantaneous growth function is

Figure 2 illustrate $r(t)$ and $C(t)$ at different values of $\nu > 0$. The parameter $r$ is a time-scale parameter and the x-axes in these plots are scaled to $\tau = rt$. We focus on the range $0 < v < 1$. During this range, $\left(1+rvt\right) ^{1/v}$ is the same as the GGM (2.2) by letting $v = 1-p$ and $C^{\prime}(t) = [r/\{1+r(1-p)t\}]C(t)$ can be re-written as $C^{\prime}(t) = r\left(1+r(1-p)t\right) ^{\frac{p}{1-p} } = rC(t)^{p},$ corresponding to the phenomenological model (2.2)$.$ When $\nu = 1,$ the linear growth $C(t) = 1+rt$ corresponds to $\overline {G}^{(mixed)}(t) = \left[ L\mathcal{[\xi]}(rt)\right] = \left(1+rt\right) ^{-1}$ recognized as the survival function of a (continuous time) power-law distribution, analogous to the discrete power-law distribution in social network analysis. $C(t)$ is a convex increasing function bounded by the exponential function from above and the linear function from below. This is the sub-exponential growth feature ^[6].

We also draw a connection to previous studies ^[28,29], which demonstrate that an initial Gamma distribution of agent susceptibility and infectivity in the SIR model leads to power-law sub-exponential growth. These results, although derived through a different theoretical pathways, aligns with our findings on frailty-induced heterogeneity using a Gamma distribution. This similarity suggests a deeper link between different approaches to modeling heterogeneity in epidemics. it is also worth noting the recent work by

2.6.2.   Inverse-Gaussian distributed heterogeneity

If $\xi(z)$ be the p.d.f. of an inverse-Gaussian distribution with $E(z) = 1$ and variance $var[z] = v,$ the Laplace transform is $L\mathcal{[\xi ]}(s) = e^{\left(1-\sqrt{1+2vs}\right) /v}.$ Letting $C(0) = 1,$ the expression (2.12) becomes

and $C(t)\rightarrow e^{rt}$ as $v\rightarrow0.$ Meanwhile,

However, the difference between $C(t)$ given in (2.14) and $C(t)$ given in (2.13) is small:

It can be shown that, if $E(z) = 1$ and variance $var[z] = v > 0,$ $e^{\left(\sqrt{1+2vrt}-1\right) /v}\geq\left(1+vs\right) ^{1/v}$ for all $s > 0.$

It is a special case of a more general result ^[30]:

where $v$ on the left hand is the variance of an inverse-Gaussian distribution and $v_{G}$ on the right hand is the variance of a gamma distribution satisfying $v\leq v_{G}$. In both distributions, it is assumed the mean value $E[z] = 1.$ Relating (2.15) to $C(t) = \left[ L\mathcal{[\xi ]}(rt)\right] ^{-1},$ it implies that if one compares two growth curves $C_{IG}(t)$ and $C_{GAM}(t),$ generated by a frailty model with random effect distribution $\xi(z)$ chosen from an inverse-Gaussian distribution with variance $v$ and from a gamma distribution with variance $v_{G},$ respectively, if $v\leq v_{G},$ $C_{IG}(t)\geq C_{GAM}(t).$

2.6.3.   Two distributions for $\xi(z)$ where $L\mathcal{[\xi]}(s)$ lacks explicit expressions and discussion of sub-exponential growth

If we replace the left hand side in the inequality (2.15) by the Laplace transform of a log-normal distribution with variance $v$, denoted by $L_{LogN}[\xi](s),$ the following statement is still true.

To show this, we briefly introduce $L\kappa$-class distributions.

Definition 1. The distribution for $X > 0$ with p.d.f. $f(x)$ and finite mean $\mu$ is said to belong to the $L_{\kappa}$-class of distributions, if

The p.d.f. of a log-normal distribution of a random variable $X$ can be written as $f_{X}(x; \gamma, \varsigma) = \frac{\varsigma}{\sqrt{2\pi}x} e^{-\frac{\left(\varsigma\log\left(\gamma x\right) \right) ^{2}}{2}}$ ^[31] with mean $E[X] = \frac{1}{\gamma}\exp\left(\frac {1}{2\varsigma^{2}}\right)$ and variance $var[X] = \frac{1}{\gamma^{2} }e^{\frac{1}{\varsigma^{2}}}\left(e^{\frac{1}{\varsigma^{2}}}-1\right).$ The gamma distribution with mean value $\mu_{G}$ and shape parameter $\kappa$ has the variance expression $v_{G} = \mu_{G}^{2}/\kappa.$ We compare a log-normal distribution with a gamma distribution with equal mean values and let $\mu_{G} = \frac{1}{\gamma}\exp\left(\frac{1}{2\varsigma^{2}}\right).$ The condition $var[X]\leq\mu_{G}^{2}/\kappa$ becomes $e^{\frac{1} {\varsigma^{2}}}\leq1+\frac{1}{\kappa_{G}}.$ Immediately following Proposition 3.4 in ^[30], if $\varsigma^{-2}\leq\log\left(1+\frac{1}{\kappa }\right),$ the log-normal distribution belongs to the $L\kappa$-class.

Applying the above to $\xi(z)$ with mean $\mu = 1,$ then $v_{G} = 1/\kappa,$ the condition $\varsigma^{-2}\leq\log\left(1+\frac{1}{\kappa}\right)$ is $v\leq v_{G}$ and the inequality (2.16) holds even though $L_{LogN}[\xi](s)$ does not have explicit expression. Since $C(t) = \left[ L\mathcal{[\xi]}(rt)\right] ^{-1},$ it implies that if one compares two growth curves $C_{IogN}(t)$ and $C_{GAM}(t),$ generated by a frailty model with random effect distribution $\xi(z)$ chosen from a log-normal distribution with variance $v$ and from a gamma distribution with variance $v_{G},$ respectively, if $v\leq v_{G},$ $C_{IogN}(t)\geq C_{GAM}(t).$

With respect to the benchmark $v_{G} = 1,$ of the three distributions for $\xi(z)$, when $E(z) = 1$ and variance $var[z] = v\leq1,$ they all belong to the $L_{1}-$class satisfying $L[\xi](s)\leq\left(1+s\right) ^{-1}.$ This gives $C(t) = \left[ L\mathcal{[\xi]}(rt)\right] ^{-1}\geq1+rt.$ We can add the Weibull distribution for $\xi(z)$ into this class because ^[30] has also shown that the Weibull distribution with $E(z) = 1$ and variance $var[z] = v\leq1$ belongs to the $L_{1}-$class. Therefore, for four distributions for $\xi(z)$: the gamma, the inverse-Gaussian, the log-normal and the Weibull distribution, when $E(z) = 1$ and variance $var[z] = v\leq1,$ $C(t) = \left[ L\mathcal{[\xi ]}(rt)\right] ^{-1}$ satisfies $1+rt\leq C(t) < e^{rt}.$

It is worth noting that a power-law sub-exponential diffusion growth has been demonstrated on the recent COVID-19 epidemic (for example in ^{[32,33,34,35,36]}).

3.   Discussion

Most popular compartmental epidemic models treat populations as homogeneous, with only a few key differences considered, such as age, gender, and location ^[7]. However, capturing all the factors that influence the transmission dynamics of infectious diseases in these models is practically impossible. This is partly because some relevant factors might be unknown, a challenge frequently encountered in epidemiology. Our analysis employs the frailty model from survival analysis, offering a robust framework for addressing this issue. This model is suitable for incorporating unobserved heterogeneity into the study of epidemic growth dynamics, without the need to define the specific mechanisms. By doing so, it allows us to explore the dynamic impacts of these unaccounted-for variables on disease transmission, providing insights that are not obtained from conventional models.

The frailty model has found increasing use across various applications, notably in epidemiology. In our study, we leverage the concept of frailty to provide a straightforward explanation for the observed diversity in the early growth dynamics of epidemics across different infectious diseases ^[15]. This approach allows us to delve into the phenomenon of sub-exponential epidemic growth by acknowledging the presence of unseen internal and external factors. These factors contribute to individual-level heterogeneity, a complexity that traditional mechanistic frameworks struggle to fully encompass.

The essence of the frailty framework lies in its ability to capture information about elusive individual heterogeneity through the structure of conditional hazards and the distribution of frailty (or heterogeneity) among individuals. However, a significant challenge we face is that of identifiability. Specifically, the hazard function (2.11) when the population arises from a random mixture of heterogeneous individuals involves a baseline hazard $h_{0}(t)$ and a random effect distribution $\xi(z).$ Thus, different combinations of $h_{0}(t)$ and $\xi(z)$ may produce the same marginal hazard rate $h^{(mixed)}(t).$ To overcome this issue, we propose limiting the baseline hazard rate to a constant value, $h_{0}(t) = r$ and also confining $\xi(z)$ to a specific family of distributions, such as the gamma distribution.

By incorporating individual unobserved heterogeneity into standard mechanistic compartmental epidemic models through the frailty framework from survival analysis offers a significant opportunity to understand the consequences of heterogeneity on infectious disease dynamics and control. Moreover, the use of different statistical distributions to represent frailty highlights the flexibility and robustness of our approach. The connection to previous work ^[28,29] further underscores the relevance and potential of incorporating heterogeneity in epidemic modeling. For instance, Novozhilov (2008) examined the spread of epidemics in a closed heterogeneous population by introducing a model that incorporates individual variability in susceptibility and infectivity through an initial Gamma distribution. He demonstrated that this variability leads to power-law sub-exponential growth, highlighting how heterogeneity can significantly alter epidemic trajectories.

While incorporating individual heterogeneity into epidemic models enhances their accuracy and predictive power, it also increases their complexity. This complexity can pose challenges in terms of model availability and computational feasibility. Therefore, a careful balance must be struck between the level of detail included and the practical utility of the model. Our approach aims to provide a nuanced understanding of epidemic dynamics through the concept of frailty, while maintaining a level of simplicity that ensures the model's applicability and ease of use.

In this paper we have distinguished the growth rate $r$ as a scale parameter of time in the growth curve function from the instantaneous growth function $r(t)$ given in (2.1). Importantly, our work here and a recent study ^[37] suggests that mathematical analyses that rely on the integration of mechanistic epidemic models and the frailty model are often feasible. For instance, in ^[37] authors incorporated the frailty framework into mechanistic models to explore the impact of individual variability in the level of compliance on a control intervention. Immediate research in this direction include studies that extend specific mechanistic models with the frailty framework to better understand the impact of individual unobserved heterogeneity on estimates of the reproduction number, epidemic size, stochastic extinction as well as the accuracy of epidemic forecasts.

The proposed frailty model provides a comprehensive mathematical framework to incorporate unobserved heterogeneity into epidemic growth models. By doing so, it bridges the gap between phenomenological observations and mechanistic theories, offering a more nuanced understanding of epidemic spread. This approach aligns with related work in the field, such as Novozhilov (2008) ^[28], and suggests promising directions for future research in modeling epidemic heterogeneity.

In summary, we demonstrate how the variability in early epidemic growth, ranging from sub-exponential to exponential, can be mathematically derived by incorporating individual unobserved heterogeneity into simple mechanistic models. This is achieved using the frailty model from survival analysis. By adopting the standard assumption that heterogeneity is gamma-distributed, we have derived a generalized growth equation based on the variance of the frailty distribution. These findings are in line with previous work ^[6], emphasizing the importance of including detailed individual-level characteristics, including dynamic contact networks, in transmission models to aim to capture sub-exponential epidemic growth dynamics.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Conflict of interest

The authors declare there is no conflict of interest.

Appendix.   Order of variability of non-negative random variables

If $T_{1}$ and $T_{2}$ have equal mean value $\mu_{1} = \mu_{2}$ (should they exist)$,$ a general description of variability is based on "marjorization". Let $f_{1}$ and $f_{2}$ are the corresponding p.d.f. for $T_{1}$ and $T_{2}$, the verbal description for $T_{2}$ being more dispersed (spread out) than $T_{1}$ is reflected in Figure A1 about the change of signs between. $f_{1}$ and $f_{2}$ and their corresponding survival functions $\overline{F}_{1}$ and $\overline{F}_{2}.$

The verbal and graphic description gives the following definition:

Definition 2. $T_{1}\underset{cv}{\leq}$ $T_{2}$ if and only if $\mu_{1} = \mu_{2}$ plus the following two statements:

1) $f_{2}(x)-f_{1}(x)$ has two sign changes and the sign sequence is: $+, -, +$ (see Figure A1).

2) $\overline{F}_{1}(x)-\overline{F}_{2}(x)$ has one sign change and the sign sequence is: $+, -$ (see Figure A1).

The above definition is verbal and has been proven to be equivalent to the following definition in mathematical terms:

Definition 3. $T_{1}\underset{cv}{\leq}$ $T_{2}$ if $E[\Psi(T_{1})]\leq E[\Psi(T_{2})]$ for all convex functions $\Psi(x)$ for which these expectations exist.

Therefore, the order $T_{1}\underset{cv}{\leq}$ $T_{2}$ is called the convex order.

$1)$ The convex order implies the ordering according to variance $var(T) = E[(T-\mu)^{2}]$ because $x^{2}$ is a convex function.

$2)$ It also implies the ordering according to the Laplace transform $E[e^{-sT}]$ for all $s > 0$, because $e^{-sx}$ is a convex function.