Fabienne Krauer1, Tor Erlend Fjelde2, Mihaly Koltai1, David Hodgson1, Marina Treskova-Schwarzbach3, Christine Harvey4, Mark Jit1, Ole Wichmann3, Thomas Harder3, Stefan Flasche1
1 Centre for Mathematical Modelling of Infectious Diseases, London School of Hygiene & Tropical Medicine, London, UK
2 Computational and Biological Learning, University of Cambridge, UK
3 Robert Koch Institut, Berlin, Germany
4 Health Protection NSW, NSW Ministry of Health
Background:
Respiratory syncytial virus (RSV) is a leading cause of respiratory tract infections and bronchiolitis in young children. The seasonal pattern of RSV is shaped by short-lived immunity, seasonally varying contact rates and pathogen viability. The magnitude of each of these parameters is not fully clear. The disruption of the regular seasonality of RSV during the COVID pandemic in 2020 due to control measures, and the ensuing delayed surge in RSV cases provides an opportunity to disentangle these factors and to understand the implication for vaccination strategies. A better understanding of the drivers of RSV seasonality is key for developing future vaccination strategies.
Methods:
We developed a mathematical model of RSV transmission, which simulates the sequential re-infection (SEIRRS4) and uses a flexible Von Mises function to model the seasonal forcing. Using MCMC we fit the model to laboratory confirmed RSV data from 2010-2022 from NSW while accounting for the reduced contact rates during the pandemic with Google mobility data. We estimated the baseline transmission rate, its amplitude and shape during RSV season as well as the duration of immunity. The resulting parameter estimates were compared to a fit to pre-pandemic data only, and to a fit with a cosine forcing function. We then simulated the expected shifts in peak timing and amplitude under two vaccination strategies: continuous and seasonal vaccination.
Results:
We estimate that RSV dynamics in NSW can be best explained by a high effective baseline transmission rate (2.94/d, 95% CrI 2.72-3.19) and a narrow peak with a maximum 13% increase compared to the baseline transmission rate. We also estimate the duration of post infection temporary but sterilizing immunity to be 412 days (95% CrI 391-434). A cosine forcing resulted in a similar fit and posterior estimates. Excluding the data from the pandemic period in the fit increased parameter correlation and yielded less informative posterior distributions. The continuous vaccination strategy led to more extreme seasonal incidence with a delay in the peak timing and a higher amplitude whereas seasonal vaccination flattened the incidence curves.
Conclusion:
Quantifying the parameters that govern RSV seasonality is key in determining potential indirect effects from immunization strategies as those are being rolled out in the next few years.
The main model code is written in Julia and fitted with Turing. You can download Julia here.
This repository contains all Julia and R code for this analysis as well as the input data and the final posterior samples (traces). The main scripts are:
| script name | description |
|---|---|
| prep.R | preparation of input data |
| results.R | produces all final figures |
| scripts/fit.jl | fits the main model to the weekly RSV cases |
| scripts/visualize.jl | to quickly visualize the results of the selected fit(s). Produces a PDF file |
| scripts/postprocessing.jl | simulates the posterior predictive of the selected fit(s). Produces several csv files and figures |
| scripts/vaccsim.jl | simulates the theoretical disease dynamics under different vaccination scenarios, based on the posteriors from the selected fit(s). Produces several csv files |
| scripts/loo.jl | calculated the Pareto-smoothed importance sampling leave-one-out cross-validation (LOO), based on the posteriors from the selected fit(s). |
| src/models.jl | Contains the turingmodel |
| src/differential_equations/ | Contains the ODE models for the fitting (SEIRRS4) and the vaccination simulatino (SEIRRS4Vacc and SEIRRS4Vacc_pulse) |
To get started with any julia file, run the following from the commandline:
julia --project -e 'using Pkg; Pkg.instantiate()'To start the fitting process with the default settings, run the following from the julia console:
julia> using DrWatson
julia> _args = ["seed", "--betafunc", "--symb", "--verbose"]
julia> include(scriptsdir("fit.jl"))Alternatively you can run the following line in the console:
julia --project scripts/fit.jl seed --betafunc=mises --symb --verbose
seed is any seed no. you prefer, betafunc refers to the seasonal forcing function (mises or cosine) and symb makes a symbolic version of the init ODE model for faster fitting. verbose returns more detailed information of what the script is doing at any time, we recommend to turn it on. fit.jl takes multiple arguments to adjust the solver and sampling settings as well as the model arguments (see code for arg options). The fitting script stores the model, the trace, the arguments and the state after adaptation in a date-time stamped folder within the output directory.
The current fitting script does not allow multithreading because the sampling progress cannot be displayed in multithreading mode. Instead, we suggest to open multiple instances of julia to run the script several times (with different seeds). In our experience, running four instances in parallel does not impact performance (but this depends of course on the performance of your machine, please check that the number of instances does not exceed the number of available cores on your machine).
The following settings/seeds reproduce our fitting results:
| model | args | seeds |
|---|---|---|
| Mises 2010-2022 | --betafunc=mises --symb --verbose | 25,95,104,222,303,462,1886,1983,9333,68136 |
| Mises 2010-2019 | --betafunc=mises --symb --verbose --endyear=2019 --target-acceptance=0.95 | 95,104,222,303,462,1886,9333,68136 |
| Cosine 2010-2022 | --betafunc=cosine --symb --verbose | 25,644,874,7425,8319,336447 |
| Mises 2010-2022, without AR data | --betafunc=mises --symb --verbose --noar | 25,95,104,222,303,462,1886,9333,68136,83411 |
To run the diagnostics/visualization of the fitting results, run the following from the console:
julia --project scripts/visualize.jl --verbose --theme=ggplot2 output/foldernamewhere foldername is the name of the folder containing the fitting results. You can also combine the results from different fits with the same settings with
julia --project scripts/visualize.jl --verbose --theme=ggplot2 --prefix=folderprefixwhere folderprefix is the prefix of the folders with the fits to combine
To generate the posterior predictive of the fit, run the following from the console:
julia --project scripts/postprocessing.jl --verbose output/foldernamewhere foldername is the name of the folder containing the fitting results. You can combine the results from different fits as above.
To simulate the vaccination scenarios, run the following from the console:
julia --project scripts/vaccsim.jl seed output/foldername --verbose where seed is the random number see and foldername is the name of the folder containing the fitting results. You can combine the results from different fits as above.
Do you have questions about the study or the code? Contact Fabienne Krauer or Tor Fjelde by email or open a GH issue.