## Abstract

While seasonal variation has a known influence on the transmission of several respiratory viral infections, its role in SARS-CoV-2 transmission remains unclear. As previous analyses have not accounted for the implementation of non-pharmaceutical interventions (NPIs) in the first year of the pandemic, they may yield biased estimates of seasonal effects. Building on two state-of-the-art observational models and datasets, we adapt a fully Bayesian method for estimating the association between seasonality and transmission in 143 temperate European regions. We find strong seasonal patterns, consistent with a reduction in the time-variable *R*_{t} of 42.1% (95% CI: 24.7% – 53.4%) from the peak of winter to the peak of summer. These results imply that the seasonality of SARS-CoV-2 transmission is comparable in magnitude to the most effective individual NPIs but less than the combined effect of multiple interventions.

## 1. Introduction

Since the onset of the COVID-19 pandemic, the role of seasonal variation in SARS-CoV-2 transmission has received significant scientific and political attention [1]. Understanding seasonal patterns is vital, as it enables more accurate inferences about current trends in transmission and how they may change over the longer term. For example, a proper understanding of seasonality can help policymakers avoid attributing declining incidence over the summer to population immunity alone, when in fact seasonality may be playing a meaningful role.

While seasonal variation is well-established for many respiratory viral infections [2], and some studies have suggested associations between temperature, humidity, and COVID-19 incidence [3, 4, 5], other analyses have failed to show a robust role of climate and weather [6, 7, 8]. A recent review found that the evidence remains inconclusive [9].

A further complication is that temperature, humidity, and UV radiation plausibly affect transmission and incidence through a range of biological and epidemiological mechanisms [2, 10]. These include virus stability and viability [11, 12], host susceptibility and immune response [13, 14], human behaviour [15, 16], and social factors such as holidays and school calendars [17, 18]. This multitude of plausible causal pathways makes it exceedingly difficult to disentangle the influence of various seasonal factors, particularly given the extensive multi-collinearities and interactions between environmental, biological, and behavioural elements [19]; see Appendix A for an overview of the various causal pathways. As Lofgren and colleagues note in the context of influenza, “the myriad theories accounting for seasonality (…) suggest that the elegant and predictable periodicity of nonpandemic influenza is caused by a less-than-straightforward interaction of many different factors,” meaning that “recognition of this complexity, as well as the likelihood that seasonality arises from many different factors, is essential for continued examination and elucidation of seasonality” [15]. Given the severe methodological challenge of disentangling these interrelated factors, a more tractable solution is to approach seasonality holistically with the purpose of understanding its overall effects. In this study, we infer a single seasonality parameter, describing the amplitude of the yearly variation in the time-variable reproduction number, *R*_{t}, for one climate region. While this single parameter does not disentangle the individual effects that comprise a seasonal profile, it accounts for the overall magnitude of the seasonal effect on SARS-CoV-2 transmission and thereby provides valuable insights for long-term policy planning.

Since both COVID-19 incidence and the presence of governmental non-pharmaceutical interventions (NPIs) have waxed and waned in consecutive waves since early 2020, adjusting for NPI effects is crucial for any effort to infer the influence of seasonality on transmission, yet previous analyses of environmental drivers have largely not done so [9, 5, 3]. We accomplish this by extending two hierarchical Bayesian models of NPI effects from Brauner et al. [20] and Sharma et al. [21] to include a term representing the multiplicative seasonal influence on the effective reproduction number.

Employing a common technique in infectious disease modelling [22, 23], we assume the seasonal variation itself is described by a sinusoidal modulation. We re-analyse data from existing studies [20, 21] while restricting the scope to European regions in the temperate climate zone, where we assume the seasonality effect to be comparable both in its environmental and behavioural causal components.

## 2. Methods

Published hierarchical Bayesian models of NPI effects [20, 21, 24] typically assume the time-variable reproduction number *R*_{t} to be a product of *R*_{0}, the “natural” reproduction rate without mitigations, multiple terms representing the effects of interventions, and noise terms modelling growth speed variance, of the form:
where exp(−*α*_{i}) is the estimated effect of intervention *i* and *x*_{i}, *t, l* are indicator variables of interventions implemented on a given date in a given location. Daily reproduction rates are then typically connected to observed data on cases or deaths by growth rates, compartmental epidemiological models, renewal processes, etc. The noise term *N*_{t,l} then varies with the model, e.g. a log-normal multiplicative factor in Brauner et al. [20] and a random walkbased multiplicative factor in Sharma et al. [21]. Mechanistically, the noise term can be intuitively thought of as a random effect that accounts for residual variation not captured by the NPI effects.

To account for seasonality, we substitute *R*_{t,l} with (adjusted for seasonality) and let the model infer a single seasonality amplitude parameter *γ* along with its other parameters. This minimal modification aims to preserve the demonstrated robustness of the original models [20, 21, 25]. The final seasonality amplitude estimate is then pooled from the two models, equally sampling from their posterior distributions. Note the seasonal adjustment to *R*_{t} is *shared* across all locations *l* and therefore captures common dynamics between locations not explained by the location specific noise terms or NPIs.

We model seasonality as a sinusoidal multiplicative factor Γ(*t*) to *R*:
where *γ* is the intensity (amplitude) of the seasonal effect, *d*_{γ} is the day of the year of the highest seasonal effect on *R*, and *d*_{0} is the first day of the respective dataset. We assume a single, common seasonal effect for countries in similar climates and relative proximity along dimensions such as income, political structure, and culture.

For both models, the time and location-specific *R*_{t,l} is replaced with seasonal :
Note that we divide Γ(*t*) by Γ(0) to have Γ(*t*)/Γ(0) = 1 at *t* = 0 and , i.e. the seasonality multiplier is normalised to 1 at the start of the window of analysis. This means that the priors over need not be adjusted in either model. For both models, we assume an uniform prior *γ* ∼ *U* (−0.95, 0.95)^{1}.

The amplitude of the cyclical seasonal variation (*γ*) can be converted to the reduction in transmission associated with going from the peak of winter to the peak of summer (i.e., peak-to-trough) as . Similarly, the amplitude can be directly converted to the expected reduction between adjacent seasons such as peak winter to midspring or mid-spring to peak summer (i.e., peak-to-mid).

Our analysis utilises January 1 as the seasonal peak day *d*_{γ}, as this date is both close to the center of a stable range of *d*_{γ} in the sensitivity analysis of Appendix E.2, as well as close to January 3, the median peak date inferred by a model with variable *d*_{γ} in Appendix E.1. Note that while we show January 1 to be a robust choice of *d*_{γ}, we are not trying to infer its exact value.

## 3. Results

Using two model structures and datasets on non-pharmaceutical interventions covering 72% of the 2020-2021 period in Europe, we estimate the seasonality parameter *γ* and the time-variable seasonal multiplier Γ(*t*) (Figure 1).

Our combined estimates from the two models are consistent with a reduction in *R* of 24.7% to 53.4% (95% CI) from January 1, the peak of winter to July 1, the peak of summer, with a median reduction of 42.1% (Figure 2, Figure C.5, Table C.2).

Modelling seasonality alongside non-pharmaceutical interventions allows us to gain a sense of the epidemiological importance of environmental factors. We find that the transition from winter to summer is associated with a reduction in transmission that is comparable to or greater than the effects of individual interventions, but less than the total effect of combined interventions (Figure 3).

Figure 3 compares our seasonality estimates to the original effect estimates from Brauner et al., as their robustness is well-established [25, 20]. Although these estimates are based on analysis that included countries outside temperate Europe, we find that restricting our analysis to temperate regions has little effect on the inferred total effect of NPIs and thus should not invalidate the comparison (Appendix E.4).

Incorporating seasonality into models of NPI effectiveness may also improve their estimates by explaining residual variation in the inferred reproduction rate. A key advance in the model proposed by Sharma et al. was the incorporation of a stochastic random walk process on the basic reproduction number to flexibly account for trends in transmission due to unobserved factors [21]. We find that including the seasonality term reduces the magnitude and asymmetry of the random walk considerably, thereby reducing the internal model variation (Appendix D). Specifically, we find that the mean square displacement (MSD) of the random walk in log-space is 0.131 for the non-seasonal model and 0.072 for the seasonal model. These results suggest a considerable amount of the residual variation can be explained by a common seasonality profile.

Estimates of seasonality and NPI effects are sensitive to modelling choices [20, 25, 21]. It is therefore vital to include a sensitivity analyses of free parameters and inputs to ensure consistent results. Relying on the demonstrated robustness of the original models, we focus primarily on the parameter that we introduce in the form of peak seasonality day. We find that the inferred mean peak-to-trough reduction in *R* varies by less than 5% across all the analysed peak seasonality dates in December and January (Appendix E.2). Although the seasonality magnitude is somewhat sensitive to setting the winter peak to different dates in February, these dates are considerably later in the year than the median peak date inferred in our sensitivity analysis, January 3 (see Appendix E.1, Figure E.9).

Since the seasonality term we introduce is directly related to *R*_{t,l} through Equation (3), we also examine the sensitivity of our results the mean initial *R*_{0} prior. We find that our results are robust to univariate variation in this parameter, with the seasonal Sharma et al. model being the most sensitive (Appendix E.3).

## 4. Discussion

The clear seasonal patterns of other respiratory viruses give us strong prior reasons to expect seasonal variation in SARS-CoV-2 transmission [2], and the strong associations we observe in temperate Europe match this expectation. While reductions in reproduction rates and case numbers are not directly comparable, another recent analysis by Chen et al. [5] infers a 64% reduction in cases from one season to the next based on a cross-sectional regression at a single point in time, similarly suggesting a significant role of environmental factors. The general magnitude of our results is also in line with previous assumptions about the reduction of SARS-CoV-2 *R*_{0} between winter and summer peaks, which range from 10% to 40% depending on the degree of seasonal forcing [1]. Moreover, recent analyses have suggested a role of environmental factors in the B.1.1.7 lineage transmission intensity and that such factors may differentially affect the transmission of different variants of concern [4].

It is important to note that our results are not inconsistent with widespread outbreaks in warmer regions, nor do they imply that temperate regions cannot face surges in transmissions during summer periods [8]. Despite moderate seasonal forcing, the time-variable reproduction number can remain well above 1, as was in fact the case in parts of Europe during the studied period and is clearly still the case in several warmer regions across the world.

Moreover, this study utilised variation in environmental and behavioural factors across time while holding the climate zone constant, and the observed results may not directly translate to comparisons across regions holding the season constant. In other words, the relationship between cooler periods and transmission within the temperate zone does not necessarily imply an exactly similar association between regional climate and transmission rates at any given point in time. This is because latitude is correlated with a wide range of epidemiological, demographic, and societal factors, each of which may affect transmission.

One major limitation of our analysis is that it relies on data from only one period of seasonality. We present the inferred seasonality estimates as the best estimate given the available data. Moreover, since our analysis focused exclusively on European regions in the temperate climate zone, the findings may not generalise to other climates, particularly as we have not identified the relative contributions of different causal mechanisms. Other respiratory infections show less seasonality in tropical regions relative to temperate regions as well as seasonal patterns with different peak timings, for example, during the monsoon season [2, 26]. Further research can shed light on the extent to which this is the case for SARS-CoV-2, and on the interaction between seasonality and latitude within climate regions.

More generally, this observational study demands caution when drawing conclusions about causality. Our analysis did not attempt to disentangle the various plausible causal pathways through which seasonality may affect transmission, and both environmental and behavioural factors can vary over the years. For example, behavioural patterns throughout the first year of the pandemic were likely exceptional, and while some behavioural changes are closely tied to modelled NPIs and thus do not bias our analysis, other relevant behavioural aspects may differ in subsequent years. Consequently, a granular focus on specific factors such as temperature, humidity, and behaviour is required for short-term prediction to inform policy.

Notwithstanding these limitations, the parsimonious seasonality form may be adequate to understand variations over time and aid long-term policy planning. Even without disentangling the underlying factors, incorporating seasonality can augment modelling efforts to more reliable anticipate changes in transmission patterns, particularly when adjusting for important factors such as non-pharmaceutical interventions.

For such forward-looking analyses of SARS-CoV-2 seasonality, it should be noted that our inferred seasonal associations do not include two factors that play significant roles in the seasonality of other respiratory viruses. First, we treat school closures, including for holidays, as NPIs in our model due to the role of closing educational institutions in the epidemic responses of many countries. This means that any effects of closing schools are attributed to the school NPI, rather than to seasonality. This is noteworthy considering that school calendars are considered an important driver of seasonality for other respiratory viruses [17, 27]. Consequently, the full extent of seasonality would likely be greater if it is construed to include school calendars. Second, the seasonal variation of some respiratory viruses, such as influenza, owes to a combination of both the direct seasonal forcing from biological and behavioural factors as well as the indirect influence of waning population immunity [28]. Given what is known about the robustness of acquired immunity within the first year of SARS-CoV-2 infection [29], the patterns we observe likely owe almost entirely to seasonal forcing. Going forward, the long-term seasonality of SARS-CoV-2 will depend in part on developments in population immunity as well as on the emergence of variants.

## 5. Conclusion

Failing to account for seasonality may lead to grave policy errors or Panglossian outlooks. For instance, a reduction in transmission over the summer may be misinterpreted as the result of herd immunity [30], and so lead to inadequate preparation for a resurgence during the colder months. Overestimating the role of environmental factors may be equally perilous if policymakers anticipate a greater reduction due to seasonality than will actually occur.

## Data Availability

Spring 2020 dataset has been published in Brauner et al.: "The effectiveness of eight nonpharmaceutical interventions against COVID-19 in 41 countries". Fall-winter 2020 dataset is currently available upon request. The entire implementation is publicly available on GitHub.

## Appendix A. Causal pathways for SARS-CoV-2 seasonality

### Appendix A.1. Diagram of potential causal pathways

A complex web of environmental, biological, and behavioural factors contribute to the seasonality of respiratory viruses. Recent reviews by Moriyama et al. [2] and Tamerius et al. [31] provide discussions of the role that each of these factors play, with a particular focus on influenza. Figure A.4 illustrates some of the factors and potential causal pathways.

### Appendix A.2. Existing evidence on causal pathways

For each of these pathways, theory and evidence have been presented in support of a causal relationship. However, extensive multi-collinearities and interactions complicate any effort to tease apart the exact contributions of different factors, particularly when considering population-level transmission dynamics where experimental approaches are intractable.

As Lipsitch and Viboud succinctly put it: *“Unfortunately, this potpourri of possible mechanisms places us in a kind of Popperian purgatory, in which data in support of every hypothesis exist, yet none of the hypotheses has been subjected to tests that are rigorous enough to reject it”* [19].

Table A.1 presents a (non-comprehensive) selection of evidence relating to some of the important causal pathways for viral seasonality.

## Appendix B. Datasets and implementation

*Our Brauner et al. seasonal model implementation is based on the Brauner et al. code-base [20] and can be obtained at https://github.com/gavento/covid_seasonal_Brauner together with the datasets used. We restrict the dataset of Brauner et al. [20] to the following 29 regions (out of 41 total):*

*Albania, Andorra, Austria, Belgium, Bosnia and Herzegovina, Bulgaria, Croatia, Czech Republic, Denmark, Estonia, France, Germany, Greece, Hungary, Ireland, Italy, Latvia, Lithuania, Malta, Netherlands, Poland, Portugal, Romania, Serbia, Slovakia, Slovenia, Spain, Switzerland, United Kingdom*.

Our Sharma et al. seasonal model implementation is based on the Sharma et al. code-base [21] and can be obtained at https://github.com/gavento/covid_seasonal_Sharma/. We use the Sharma et al. dataset without any modifications and with the same preprocessing, in particular, we also exclude datapoints with non-negligible novel SARS-CoV-2 strain prevalence.

## Appendix C. Detailed results

## Appendix D. Random walk noise comparison

The Sharma et al. model contains a random walk process on log *N*_{t,l}, the logarithm of a multiplicative factor in *R*_{t,l}, in order to account for continuous slow changes of *R*_{t} through unobserved external factors such as unobserved NPIs or environmental transmission factors [21].

To show the effect of modelling seasonality on the inferred random walk noise, we compare *N*_{t,l} for the two models in Figure D.6. While the random walk trajectories are comparable in width, the median trajectory of the non-seasonal model follows an increasing trend. Note that the random walk multiplier is modelled as a symmetrical random walk in log-space. Also note that the Sharma et al. model allows only weekly changes in *N*_{t,l}.

To quantify the improvement, we compute the mean squared deviation (MSD) of log *N*_{l,t}, i.e. the random-walk multiplier in log-space, across the sampled random walks. We find this MSD to be 0.131 for the non-seasonal model, and 0.072 for the seasonal model, a 45% decrease. Note that here we compute MSD as
where is the *i* -th sample of *N*_{l,t}.

We compare for the two models – the reproduction factor derived from region-specific *R*_{0} by the random walk process and by seasonality effect (in the seasonal model) but *before* applying transmission reduction of the active NPIs:

where *N*_{t,l} is the random walk noise.

Figure D.7 illustrates how the inferred are comparable for the non-seasonal and seasonal model. However, the non-seasonal model random walk is of a larger overall amplitude and has an asymmetric trend compared to the seasonal model, as shown on Figure D.6. We interpret this as an indirect evidence towards seasonality improving the quality of model fit on Sharma et al. data.

Note that the noise terms used in Brauner et al. are of different type and the model does not contain a comparable random walk noise term.

## Appendix E. Sensitivity analysis

### Appendix E.1. Inferring the peak seasonality day

To examine the fitness of our seasonality peak estimation, we place a prior of 𝒩 (Jan 1, 45^{2}) on *d*_{γ} instead of a fixed date. Figure E.8 shows the distribution of the seasonality multiplier cosine curves Γ(*t*) inferred with prior on *d*_{γ}. Figure E.9 shows both the inferred seasonality peak day *d*_{γ} and the seasonality amplitude *γ*.

Note that the estimated *d*_{γ} are shown as a model validation, illustrating the range of seasonality peak the models and data are consistent with – we do not claim the models and data can infer the peak with accuracy. Note that the inferred *γ* in the model inferring the seasonality peak is virtually unchanged relative to a fixed seasonal peak day model (Figure E.9 bottom).

### Appendix E.2. Sensitivity to peak seasonality day

We test model sensitivity to the choice of peak seasonality day *d*_{γ} for *d*_{γ} ∈ {Dec 4, Dec 18, Jan 1, Jan 15, Jan 29, Feb 12, Feb 26}. We observe that the inferred combined effect of the NPIs and the inferred seasonality are stable for *d*_{γ} in December and January with the exception of the combined NPI effect in Brauner et al. Note that Sharma et al. is particularly robust in this range.

### Appendix E.3. Sensitivity to initial *R*_{0} prior

We test our model sensitivity to the choice of the mean of the initial prior (i.e. location-specific *R*_{0} on the first day of the dataset). We analyse the prior mean in ranges similar to the sensitivity analyses in Sharma et al. [21] and Brauner et al. [20].

This analysis is motivated by the seasonal amplitude parameter *γ* being closely connected with via Equations (1) and (3). Mis-specifying the initial could be compensated for by the model e.g. by a different amplitude *γ* and therefore also the slope of Γ(*t*) in the seasonality multiplier sine curve.

In Figures E.12 and E.13 we observe the inferred combined effect of the NPIs and the inferred seasonality to be mostly stable in the Brauner et al model. However, in the Sharma et al. seasonal model the inferred seasonality amplitude and peak-to-trough reduction are mildly sensitive to the prior. (Note that both those parameters are very closely tied together.)

Note that both original models also exhibit some sensitivity of the effect of NPIs to *R*_{0} prior mean; see Figure S11 in Supplementary material of Brauner et al. [20] and Figure S13 in Supplement of Sharma et al. [21] (v1).

### Appendix E.4. Inferred NPI effects in various models

We compare the inferred total NPI effect in different models and data subsets to verify its stability: seasonal vs non-seasonal (original) models, and the original full dataset vs the dataset restricted to temperate Europe countries (“TE” in plot for Brauner et al. model).

We observe that restricting the dataset regions in Brauner et al. model has very little effect on the inferred combined NPI effect.

Switching to a seasonal model produces a small decrease (resp. increase) in combined NPI effect in Brauner et al. (resp. Sharma et al.) model. While this may be coincidental, this effect is consistent with a hypothesis that a part of the seasonality-related change in *R*_{0} (i.e. the proposed spring decrease of *R*_{0} in Brauner et al., fall increase in Sharma et al.) is in part attributed to NPI activations in both models. Recall that Brauner et al. only considers NPI activations in their model, and Sharma et al. dataset is dominated by NPI activations compared to deactivations. However note that both models do contain noise terms for growth rate and other mechanisms to model small or slow changes in *R* due to unobserved factors, so the extent of this effect remains unclear.

## Footnotes

↵

^{1}Sinusoidal seasonality is well-defined only for amplitudes −1 ≤*γ*≤1. We restrict it to 0.95 −≤*γ*≤0.95 to improve model stability.