Upper Bounds on Overshoot in SIR Models with Nonlinear Incidence

We first examine the effect of nonlinear incidence on overshoot for ODE models, where the computations can be made analytical. The equations of the Kermack-McKendrick SIR ODE model are given as follows with generic incidence term $\beta f(S)g(I)$, where $f(S)$ and $g(I)$ are functions of $S$ and $I$ respectively to be specified.

where $S$, $I$, $R\in[0,1]$ are the fractions of the population that are susceptible, infected, and recovered respectively, $\beta,\gamma\in\mathbb{R}_{>0}$ are positive-definite parameters for transmission and recovery rate respectively.

For the SIR model, the overshoot is given by the following equation:

where $S_{t^{*}}$ is the fraction of susceptibles at the time of the prevalence peak (i.e. when $I$ is maximal in value), $t^{*}$, and $S_{\infty}$ is the fraction of susceptibles at the end of the epidemic. To solve this equation, the easiest approach is to derive an equation for $S_{t^{*}}$ in terms of only $S_{\infty}$ and parameters. We do this by first setting (5) equal to 0 and solving for the critical susceptible fraction $S_{t^{*}}$.

By using the usual definition for the basic reproduction number, $R_{0}\equiv\frac{\beta}{\gamma}$, we obtain the following equation for $S_{t^{*}}$.

We can see from this equation that $S_{t^{*}}$ will have $I$ dependence unless $g(I_{t^{*}})=I_{t^{*}}$. Thus to make what follows analytically tractable, let us assume $g(I_{t^{*}})=I_{t^{*}}$. We will provide even stronger justification why $g(I)$ must take this form later in the results. This assumption of $g(I_{t^{*}})=I_{t^{*}}$ reduces the above equation to the following.

Taking this equation for $S_{t^{*}}$ (9) and the overshoot formula (7), we obtain:

Thus the main challenge now becomes a problem of finding an equation for $R_{0}$ and the inverse function $f^{-1}$. Based on previous results [35], the following outlines the general steps for calculating the maximal overshoot for a SIR model:

1. A.

   Take the ratio of $\frac{dI}{dt}$ and $\frac{dS}{dt}$. Integrate the resulting ratio. The indefinite integral requires a constant of integration, which is a conserved quantity that applies at every time point along the system’s trajectory in time.

2. B.

   Evaluate the equation for the conserved quantity at the beginning of the epidemic ($t=0$) and the end of the epidemic ($t=\infty$) using initial conditions and asymptotic values. Then, rearrange the resulting equation for $\frac{1}{R_{0}}$.

3. C.

   Find the form for the inverse function, $f^{-1}$.

4. D.

   Combine the equations for $\frac{1}{R_{0}}$ and $f^{-1}$ with the overshoot equation.

5. E.

   Maximize the resulting overshoot equation by taking the derivative of the equation with respect to $S_{\infty}$ and setting the equation to 0 to find the extremal point $S_{\infty}^{*}$. This step usually leads to a transcendental equation for $S_{\infty}^{*}$, which can be solved numerically.

6. F.

   Use the maximizing $S_{\infty}^{*}$ value in the overshoot equation to calculate the corresponding maximal overshoot.

7. G.

   Calculate the corresponding $R_{0}^{*}$ using $S_{\infty}^{*}$ and the $\frac{1}{R_{0}}$ equation.

Thus, the analytical exploration of nonlinear incidence terms of the type $\beta f(S)g(I)$ is reduced to exploring different forms of $f(S)$.

The first step is to rule out what forms for the incidence term will not work with the procedure outlined above.

We now show the principle reason why we require $g(I)=I$. We can see from calculating Step A. in (11) that any incidence term that does not take the form $g(I)=aI,a\in\mathbb{R}$, where $a$ is a real scalar, retains $I$ dependence upon simplification.

Any deviation from the form $g(I)=aI$ results in $I$ in the numerator and the denominator not completely cancelling out which will result in having to integrate $I$ with respect to $S$, which we will not be able to do analytically. Therefore, $I$ must enter linearly into the incidence term. Since $a$ can be absorbed into the $\beta$ parameter, all possible incidence terms for the purpose of calculating overshoot analytically will take the form $\beta\cdot f(S)\cdot I$.

We now turn to what restrictions there are on the form of $f(S)$. We start first with two conditions. First, we must enforce that when there are no susceptibles ($S=0$), then the incidence rate must go to zero ($\beta f(S)I=0$). Otherwise, since $I$ does not have such a restriction, violating this condition would leave open the unrealistic possibility that the model can generate infected people when there are no susceptibles available. To ensure this condition is met, we need the function on S to output zero if the input is zero.

For the second condition, for an outbreak to occur in an SIR model, we must have a minimum value for the basic reproduction number, $R_{0}$. Another way to view this condition is by inspecting (2) and recognizing that that the incidence term ($\beta f(S)I$) needs to be greater than the recovery term ($\gamma I$). Otherwise the epidemic cannot grow in size. Comparing the two terms leaves an inequality ($\beta f(S)I>\gamma I$) which can be rewritten as:

Beyond these conditions, an obvious requirement is that $f(S)$ should be a continuous function. In order to be able to calculate the maximal overshoot analytically, the function should be integratable with respect to $S$ and should also have a closed-form inverse $f^{-1}$. As we will demonstrate, non-monotonic functions for $f(S)$ are possible.

For $f(S)$, the following examples are constructed using basic functions that satisfy the above criteria:

1. 1.

   $\text{exp}(S)-1$

2. 2.

   Invertible polynomials of $S$

3. 3.

   $\text{sin}(aS)$

Conversely, there are many examples of functions that would not work. An example that satisfies the boundary conditions but that does not have a closed form inverse is $f(S)=\text{log}(S+1)$. While similar to examples listed, the following violate the boundary conditions: $\text{exp}(S),\text{log}(S)$, $\text{cos}(S)$. Examples that violate conditions of continuity include step functions of $S$ or $f(S)$ with cusps.

As an example, we will now look at a model with nonlinear incidence and apply the whole procedure previously outlined for finding the maximal overshoot.

In the previous sections, we focused on nonlinear incidence that was generated from the introduction of a nonlinearity in an ordinary differential equation model. Importantly, this fundamentally assumes homogeneity in the transmission, in that all infected individuals are identical in their ability to spread the disease further. In contrast, network models allow for heterogeneous spreading which depends on the local connectivity of each infected individual. This provides an entirely different mechanism through which nonlinear incidence can be generated compared to the ordinary differential equation models. As many real-world complexities and details can be more easily captured and explored in a network model, the consequences of those heterogeneities on the overshoot becomes more transparent through the targeted and fully-fleshed out explorations that can be done through numerical experimentation in network simulations [53, 54, 55, 9, 56]. However, a trade-off for the increased realism is that it becomes more difficult to perform analytical calculations for a general network model, so here we conduct a numerical exploration of the behavior of the overshoot in network models across a range of network structure.

The space of all possible network configurations is immense, so we must restrict ourselves to analyzing a particular subset of possible networks. Here we explore what happens to the overshoot when the contact structure of the population is given by a network graph that is roughly one giant component. While it is possible to construct pathological graphs that produce very complex dynamics, we consider more classical graphs here. Using a parameterization of heterogeneity ($\sigma$) given by Ozbay et al. [57] and the configuration model of Newman [58] to randomly construct networks (see Methods for details), we simulated epidemics on a spectrum of networks with structure ranging from the homogeneous limit (well-mixed, complete graph) to a heterogeneous limit (heavy-tailed degree distributions). The spectrum of distributions shapes across the space of $\sigma$ used to generate the degree distributions can be seen in the supplement. As the dynamics of the epidemic on a network are stochastic and occur in discrete-time here, they are not parameterized by $R_{0}$ as in the ODE models. Instead, we use an analogous parameter that we will denote as a basic reproduction number for networks ($R_{0,network}$) and is defined as follows:

The transmission probability ($\tau$) and recovery probability ($\rho$) are directly analogous to their counterparts ($\beta$ and $\gamma$ respectively) in $R_{0}$ in that that correspond to transmission and recovery parameters, albeit for a stochastic model. The inclusion of the mean network degree ($\langle k\rangle$) as a scalar makes intuitive sense as that represents the average number of potential neighbors an infected node could potentially infect at the beginning of the epidemic, which is analogous to the interpretation of $R_{0}$ as the number of secondary infections. We observed what happens to the overshoot on these different graphs as we changed $R_{0,network}$.

On a network model, where contact structure is made explicit, the homogeneous limit is a complete graph, which recapitulates the well-mixed assumption of the Kermack-McKendrick model. It is not surprising then that the overshoot in the homogeneous graphs ($\sigma\approx 0)$ peaks also around 0.3 (Figure 3, $purple$ and $blue$), which coincides with the analytical upper bound previously found in the ordinary differential equation Kermack-McKendrick model [35]. However, while the homogeneous graphs shown here are highly regular (i.e. symmetric in the average number of contacts each node has), the network’s connectivity is not close to complete as the mean degree is significantly less than the network size. Thus, the average overshoot is lower than would be the case for a complete graph.

We also see that increasing contact heterogeneity (i.e. $\sigma\rightarrow 1$) qualitatively suppresses the overshoot peak both in terms of the overshoot value and the corresponding $R_{0,network}$. Furthermore, increased heterogeneity also flattens out the overshoot curve as a function of $R_{0,network}$. We see that for more heterogeneous graphs, the overshoot is larger at very low $R_{0,network}$. For some intermediate values of contact heterogeneity (Figure 3, $yellow$ and $orange$), the overshoot shows larger overshoot than the homogeneous case for large $R_{0,network}$. This would suggest a larger public health hazard for a larger range of communicable diseases for networks with structure in this regime. The intuition for the larger overshoot at high $R_{0,network}$ in more heterogeneous networks is that both the peak occurs earlier and that a significant number of cases can occur in the periphery of a network (Supplemental Materials). The probability of being connected to a high-degree node increases as the heterogeneity of the network increases, which drives the peak of infections to occur sooner. In the second phase of the epidemic as the epidemic expands out to the periphery the epidemic burns more slowly as individuals in the periphery have fewer neighbors.

However, this trend of a overshoot distribution with a big right tail does not monotonically increase with the contact heterogeneity. We see that at very high amounts of contact heterogeneity, the overall area under this overshoot curve decreases. This can be partially explained by the fact that for very heterogeneous networks, a significant fraction of the population have no neighbors at all. This limits the amount of people that can potentially be infected and caps the outbreak size and subsequently overshoot size.

In Figure 3, we presented the results of SIR simulations of epidemics run on graphs of size N = 200 and the mean degree is 7, where the parameter of interest is $\sigma$. Each curve represents a different value of the graph heterogeneity. We implemented the following procedure from [57] for generating graphs as a function of a continuous parameter ($\sigma$):

The following simple procedure generates a graph that has the desired heterogeneity:

1. 1.

   Choose values of $\sigma$ (heterogeneity), the mean node degree (which can be set through the relationship $\text{Mean}=\lambda\Gamma\left(1-\frac{1}{\text{log}(\sigma)}\right)$ via $\sigma$ and an appropriate scale parameter ($\lambda$), and $N$ (number of nodes).

2. 2.

   Draw N random samples from the following distribution using $\sigma$ and $\lambda$, rounding these samples to the nearest integer, since the degree of a node can only take on integer values.

   $$f(x;\lambda,\sigma)=\frac{-\text{\text{ln}}(\sigma)}{\lambda}\left(\frac{x}{\lambda}\right)^{-\text{\text{ln}}(\sigma)-1}e^{-(x/\lambda)^{-\text{\text{ln}}(\sigma)}};x\geq 0;\sigma\in(0,1];\lambda>\mathbb{R}^{+}$$

3. 3.

   With the sampled degree distribution from the previous step, now use the configuration model method [58] (which samples over the space of all possible graphs corresponding to a particular degree distribution) to generate a corresponding graph.

This yields a valid graph with the desired amount of heterogeneity as specified by $\sigma$.

We implemented the following simulation procedure from [57] for implementing an SIR epidemic on a graph:

Given a graph $G(\sigma)$ of heterogeneity $\sigma$, fix a transmission probability $\tau$ and recovery probability $\rho$:

1. 1.

   At time $t_{0}$, fix a small fraction $f$ of nodes to be chosen uniformly on the graph and assign them to the Infected state. The remaining $(1-f)$ fraction of nodes start as Susceptible.

2. 2.

   For each $i\in[1,T]$, for each pair of adjacent S and I nodes, the susceptible node becomes infected with probability $\tau$.

3. 3.

   For each $i\in[1,T]$, each infected node recovers with probability $\rho$.

4. 4.

   For each simulation, record the overshoot as the difference between the number of susceptibles at the peak of infection prevalence and the end of the time dynamics.

5. 5.

   Repeat steps (1-4) $n$ times for each value of $\tau$.

6. 6.

   Repeat steps (1-5) for each value of $\sigma$.

Code for the network simulations is provided in the Supplemental Information.

Where needed, equations were solved numerically using the $ode45$ numerical solver in $MATLAB$. Code for the network simulations is provided in the Supplemental Information.