Conditions for a second wave of COVID-19 due to interactions between disease dynamics and social processes

Mechanisms causing a second wave

At our baseline parameter values, time series of infection prevalence I(t) and support for closure x(t) exhibit nontrivial time evolution, including a second wave of COVID-19 infections (Figure 1a, b). These results illustrate the basic mechanisms underlying the model dynamics. As infection prevalence grows, support for closure rises and eventually crosses the 50% threshold by t = 80 days. After this, infection prevalence peaks and begins to decline. Support remains at a high plateau for a period of two months, after which support for closure wanes, causing restrictions to be lifted by t = 160 days. Shortly thereafter, prevalence begins to rise again. Support for closure correspondingly rises again, but not quickly enough to prevent a second wave of COVID-19 with a peak at t = 240 days that is higher than the peak of the first wave.

• Download figure
• Open in new tab

Figure 1:

(left) Time series plot of the number of infected individuals in the face of social showing appearance of second wave (left) due to effects of social dynamics (right). The results are obtained for C₀ = 0.6, R₀ = 2.4,1/σ = 2.5, γ = 1/5, κ = 5, ϵ = 0.008, α = 1.0 * 24/365 and δ = 0.05 * 24.0/365 and initial conditions S(0) = (100000α−1)/1000000, E(0) = 0, I(0) = 1/1000000, R(0) = 0, x(0) = 0.25, L(0) = 0 and H(0) = 0.

Epidemiological conditions for a second wave

Our goal was to gain insight into the conditions that generate a second wave of COVID-19, and to test the robustness of the predicted result at our baseline parameter values illustrate in Figure 1. Hence we explored the model dynamics in the neighbourhood of our baseline parameter values (Table 1) using parameter planes that show how the dynamical regimes of the model vary with changes in two model parameters (one on each axis) around the baseline values. We explored two dynamical outcomes: the number of waves in the course of the entire pandemic, and the ratio of the peak height of the second wave to the peak height of the first wave.

We start by exploring the effects of variation in the epidemiological parameters β (transmission rate) and γ (inverse of the average duration of infectiousness) in Figure 2, while keeping the rest of the parameters at baseline values. The results show that one, two, or three waves are possible under variation in these parameter values. A second wave characterizes most of the β − γ plane, however. A third wave appears when 1/γ is between 5 to 7 days and β is up to 0.4/day, corresponding to R₀ < 2. (We note that most estimates place R₀ > 2 for COVID-19 [19, 29].) The second peak may be higher or lower than the first peak, depending on the β and γ parameter combinations. For R₀ > 2.4, the second peak tends to be lower than the first, while for R₀ < 2.4 is is higher.

• Download figure
• Open in new tab

Figure 2:

Parameter plane showing model dynamics as they vary with changes in the transmission rate β and the rate γ at which the infectious period ends, for the number of COVID-19 waves (a) and the ratio of peak height of the second wave to the first wave peak. Other parameters are at baseline values (Table 1).

Changes in the duration of infectiousness 1/γ and the duration of the latent stage 1/σ around baseline values do not change the number of peaks: a second wave is still observed across the range we explored (SI Appendix: Figure 2). However, the second peak is higher than the first when 1/y is between 3 to 5 days, while out of this range the second peak is lower. The lack of dependence of dynamics on σ is expected. When 1/γ < 3 days, the second peak is less severe because R₀ drops below levels that are feasible for continued transmission in the population. In contrast, when 1/γ > 5 days the second peak is less severe because a heightened R₀ causes rapid build-up of herd immunity in the first wave of infection.

Socio-economic and intervention conditions for a second wave

Next we explored the effect of variation in intervention, economic and social parameters. The parameter plane for α, the rate at which economic losses due to closure accumulate, and δ, the discounting rate for losses, shows little variation in these values across the ranges we explored (SI Appendix: Figure 3). Two waves are predicted and the peak of the second wave is higher than the first wave for almost all parameter combinations The only exception is that when α is very small, only a single wave occurs because the population is willing to tolerate economic losses indefinitely. As a result, x remains high over the entire time horizon of the simulation and COVID-19 is effectively controlled throughout this period.

The behavioural parameter κ is a measure of how quickly novel social behaviour spreads through a population as disease cases are reported. It has a large influence on the model dynamics, as represented in the κ − α, κ − ϵ and κ − C₀ parameter planes (Figures 3-5). Higher values of κ indicate that individuals imitate more quickly. At our baseline value κ = 5/day, we observe a second wave. As the value of κ increases from this baseline value, the number of waves increases from two to six or seven in all three parameter planes, unless the effectiveness of closure (C₀) is so low that the population experiences a single large epidemic that rapidly confers herd immunity to everyone (Figures 3-5). As κ is reduced sufficiently from its baseline value, the second wave is lost as expected, since we enter a parameter regime where the population responds with an unrealistic slowness to the presence of COVID-19, and it experiences a single, rapid pandemic wave that rapidly confers herd immunity. The second peak is higher than the first peak in all three parameter planes, except again when C₀ is too low to effectively flatten the curve, and a single large outbreak results. Some examples of model outcomes for three or more waves are shown in SI Appendix: Figure 4. In these extreme scenarios, the second wave can either dominate the first and third waves, or it is also possible that the peaks of successive epidemic waves increase over time until it reaches a maximum peak in the fourth wave.

• Download figure
• Open in new tab

Figure 3:

Parameter plane showing model dynamics as they vary with changes in the rate α at which socio-economic losses accrue and the rate κ that controls the social learning rate, for the number of COVID-19 waves (a) and the ratio of peak height of the second wave to the first wave peak. Other parameters are at baseline values (Table 1).

• Download figure
• Open in new tab

Figure 4:

Parameter plane showing model dynamics as they vary with changes in the rate κ that controls the social learning rate and the parameter ϵ which controls how sensitive the population is to economic losses relative to infection prevalence, for the number of COVID-19 waves (a) and the ratio of peak height of the second wave to the first wave peak. Other parameters are at baseline values (Table 1).

• Download figure
• Open in new tab

Figure 5:

Parameter plane showing model dynamics as they vary with changes in the rate κ that controls the social learning rate and the parameter C₀ which controls how effective closure is, for the number of COVID-19 waves (a) and the ratio of peak height of the second wave to the first wave peak. Other parameters are at baseline values (Table 1).

Model equations

Transmission dynamics are given by an SEIR model, modified to take physical distancing into account, where S is the proportion of susceptible individuals (“susceptible”), E is the proportion of individuals who have been infected but are not yet infectious (“exposed”), I is the proportion of individuals who are both infected and infectious (“infectious”), and R is the proportion of individuals who are no longer infectious (“removed”). The time-varying parameter C(t) captures the impact of school and workplace closure on the transmission of COVID-19. β is the baseline transmission rate in the absence of school/workplace closure, σ is the time rate at which an exposed person becomes infections, and γ is the time rate at which an infectious person recovers.

The decision-maker decides to “turn on” closure at some time t_close, and then decides to “turn off” closure when population support for closure, x(t) drops below 50%. Hence C(t) is given by: where C₀ is a combined measure of how many workplaces are closed (the remainder being essential workplaces such as hospitals) as well as the effectiveness of physical distancing in those workplaces that remain open. The percentage of the population that supports school and workplace closure, x, evolves according to an imitation dynamic that represents social learning processes, where κ is the social learning rate, ω is sensitivity to infection prevalence, and e is sensitivity to accumulated socio-economic losses L. In this equation, support for closure goes up when the prevalence of infection goes up, but it declines when the accumulated socio-economic losses, L, become too large. The quadratic term x(1 − x) represents a social learning dynamics where individuals sample others at some rate, and they change opinion based on the utility difference ωI − ϵL. (A full derivation of this type of differential equation appears in Ref. [43]). We can absorb ω into κ, yielding: and then, setting κ′ = κω and ϵ′ = ϵ/ω and dropping the primes for simplicity we obtain

Finally, the variable L is a phenomenological representation of accumulated socio-economic losses obeying where α controls the rate at which school and workplace closure impacts socioeconomic health of the population, and δ is a decay rate that represents adjustment to baseline losses.

Parameterization

A full list of parameter definitions, baseline values, and literature sources appears in Table 1. The basic reproduction number R₀ and the transition rates σ and γ, were set based on COVID-19 epidemiological literature [16, 25, 27, 28, 34]. We note that the last compartment of the model, R, does not correspond to a stage of illness preceding recovery but rather a stage of infectiousness [2], which wanes quickly after the imposition of case isolation, in addition to the decline in viral shedding after the first 5 days [15]. Moreover, the infectious stages is preceded by a latent stage in which the virus is still replicating inside its new host until it can reach a level where the host can transmit the infection to others. These features of COVID-19 disease history guided our choice of γ and α.

The social parameters κ and ϵ were calibrated. While κ dictates how quickly novel social behaviour spreads through a population, ϵ dictates how sensitive the population is to changes in case reports relative to socio-economic losses. κ can be estimated from behaviour early in the epidemic when socio-economic losses are small as dx/dt ≈ kx(1 − x)I. For I(t) we used reports of confirmed positive cases from the early stages of the United States epidemic, adjusted by a case under-ascertainment factor of 8.7 in the United States [26] in order to estimate the true prevalence of infections I(t) in the model. We used 21 January 2020 as the initial date of the epidemic, when the first case of COVID-19 was reported in the United States. Most populations rapidly adopted physical distancing measures against COVID-19. Gallup polls indicate that 59%/79%/92% of the United States public avoided going to events with large crowds, as of 13-15 March/16-18 March/20-22 March 2020 respectively [40]. Similarly, 30%/54%/72% avoided public places, and 23%/46%/68% avoided small gatherings [40]. Taking the average of these responses across the three question types, we obtain that x(52) = 0.373, x(59) = 0.597, x(66) = 0.773 where time is measured in days since January 21. Finally, we shifted these points forward 14 days as physical distancing at time t will not be reflected in infection data until t + 14 due to the delays in testing and symptom recognition. We assumed x(0) = 0.25 when fitting to these three data points using least-squares minimization for the κ calibration. Accordingly, we fit the model to the number of cases in the United States (SI Appendix: Figure 1b). The fitted infection trajectory is in good agreement with reported cases in the US during the initial phase of the epidemic leading up to 4 April 2020. For the special case where there is an absence of any control measures, the model predicts that about 80% of the population becomes infected by the end of the outbreak (SI Appendix: Figure 1a).

The remaining two parameters, ϵ and α, were calibrated to obtain the result that x remains high after the initial surge in support for closure, but begins to drop after two months. This period of time was based on the observation that two months that have elapsed since the declaration of the national emergency in the United States on 13 March 2020, and the process of re-opening state economies that has unfolded over the month of May 2020. These two parameters control the timescale of lifting school and workplace closure based on its socio-economic impacts. Finally, the parameter δ was set such that δ = 0.05α on the basis of commonly used discounting rates in economics and assuming that economic losses accumulated through the αC(t) term would be discounted at a rate of 5% per year [20].

In order to illustrate curve flattening and show that the model has the expected response to reduction in the transmission rate due to closure, we generated model timeseries of I(t) for the special case where closure is applied throughout the entire outbreak. The epidemic curve for different values of the closure efficacy C₀ is shown, ranging from C₀ = 0 (no intervention) to C₀ = 0.6 (SI Appendix: Figure 1c). The timeseries show that the epidemic curve is flattened and delayed as closure becomes more efficacious, which reduces peak demand for intensive care beds and buys time for developing pharmaceutical interventions like vaccines and antiviral drugs, improving testing capacity, and establishing novel approaches to patient care. For the remainder of our analysis, to determine C₀ it was assumed that C₀ should be large enough to bring the effective reproduction number R_eff below 1, reflecting the observed success in multiple jurisdictions where physical distancing and closure have maintained R_eff < 1 [3, 8, 32]. Hence we chose C₀ = 1 − 1/R₀ based on the elimination threshold for the SEIR model [18]. We also assumed t_close = 20 days but in practice, our second requirement that x ≥ 50% was not reached until after 20 days in all of the model simulations.