As mentioned above, the high dimensionality of the RSV transmission model with two levels of social mixing was a limiting factor on the possible complexity of the compartmental framework representing the possible combinations of age and disease state (see appendix 2). In order to both capture the structure of the population in households and incorporate finer-grained information about the ages of the modelled individuals, we calculated empirical joint distributions for the proportion of individuals of different ages in various household sizes, and whether that household contained an under-one year old. We did not restrict the age categories of this joint age-and-household distribution to just under-one or over-one, instead preferring finer-grained age categories: (i) each month of first year of life, (ii) each year of life aged 1–18 and (iii) 18+ years old. We used the Kilifi health and demographic surveillance system (KHDSS; Scott et al., 2012) to construct the joint distributions, which records for each individual a unique person ID, a birth date, immigration into the KDHSS date(s), out-migration from the KHDSS date(s), and a unique building ID for where they live during their time in the KHDSS. By combining this data we could calculate,

where $N_t(a,n,U)$ was the number of individuals on day $t$ who were jointly in age category $a$, lived in a household of size $n$, which either contained at least one under one year old ($U=1$) or not ($U=0$), and $N_t$ was the total population size on day $t$. The joint distribution changed over time, we calculated ${\mathbb{P}}_t(a,n,U)$ for a series of year-start days t = 1 st Jan 2000, 2001,…, 2016. We then used ${\mathbb{P}}_t$ as representative for the rest of the year. Because the exact birth dates where missing for a large number of people, and for model simplicity, we assumed that all U1 individuals aged to become O1 individuals at a constant rate 1 per year, which was equivalent to assuming that given that the exact age of an U1 individual was uniformly distributed between 0 and 1 years old, independently of the U1’s household configuration.

The dynamic model of transmission tracks whether individuals are under-one or over-one years old; however, for estimating the risk of disease per infection it was useful to use the conditional age distribution for the finer-grained age category of an individual based on her dynamic model age category $a<1\text{year}$ or $a>1\text{year}$, her household size and whether the household contained an U1 or not, for example,

The conditional distributions for an individual’s household size and whether they lived in a household containing an U1 based on their age were constructed similarly. The reason we included a variable indicating whether the household of the individual contained an under one or not was because it was important to capture the pathway to transmission to the under-one year olds most at risk of disease due to contracting RSV.

The fundamental unit of the RSV transmission model developed for this paper was the household. Each household was described by the number of each type of individual inhabiting it, which we call the household configuration. The type of individual within each household was identified by her RSV disease state and age category. The RSV transmission model described the dynamics of the number of households that were in each possible household configuration using an approach introduced by House and Keeling, 2008. Mathematically, the number of households in a given household configuration at time $t$ was denoted $H_{s_1,i_1,r_1,s_2,i_2,r_2}(t)$, referring to the household configuration with exactly $s_1$ U1 susceptibles, $i_1$ U1 infecteds, $r_1$ U1 recovered, $s_2$ O1 susceptibles, $i_2$ O1 infecteds, and $r_2$ O1 recovereds. In order to limit the number of possible household states, we included only households of total size ten or less with two or fewer under ones. We chose these limits on the household size based on capturing ≈99% of the U1s in the population, and therefore the pathway to them catching RSV (appendix 2). There were 1926 possible household configurations in the RSV transmission model. The vector $\mathit{𝑯}(t)$ of number of households in each possible household configuration evolved according to the semi-linear ODE:

Each term describing the vector field of Equation (3) corresponded to a dynamic component of the model:

1. RSV transmission within households, recovery of infected individuals, loss of immunity of recovered individuals, aging from U1 to O1 and turnover in household occupancy due to births and individuals leaving the household ($A_t\mathit{𝑯}(t)$).

2. RSV transmission between households due to age-group specific mixing (${\mathit{𝒇}}_t(\mathit{𝑯}(t))$).

3. Change in household numbers due to population flux, (${\mathit{𝝆}}_t(\mathit{𝑯}(t))$).

See appendix 2 for further details. The force of infection due to transmission within a household of generic configuration ($s_1,i_1,r_1,s_2,i_2,r_2$) was density dependent; that is the person-to-person infection rate in the household did not depend on household size,

where $\tau$ is the basic within-household transmission rate, ${\iota }_2$ is the relative decrease in infectiousness of O1s compared to U1s, and $\beta (t)$ is the seasonal variation in the transmission rate of RSV (see appendix 1). Transmission outside of the household within the wider community was assumed to be based on the finer-grained age categories introduced above. The conditional age distributions of the individuals allowed us to construct matrices ($P_{H\to A,t}$) to convert between the household configuration vector into a vector of number of infected individuals in each age category, weighted by their relative infectiousness, for any time $t$ during the simulation: $\mathit{𝑰}(t)=P_{H\to A,t}\mathit{𝑯}(t)$ (appendix 2). The force of infection on each individual due to age-based mixing in the community was,

where $T$ was the community infection rate matrix and $N(t)$ was the total population size at time $t$. In this formulation, the rate at which an infected in age group $b$ creates infectious contacts in the community with individuals of age group $a$ is $T_{ab}N(a,t)/N(t)$ where $N(a,t)$ is the number of individuals in age group $a$ at time $t$(Keeling and Rohani, 2008). The force of infection on an individual within a given household was calculated using matrices constructed from the conditional distribution of an individual’s household type given her age, ${\lambda }_{com}=P_{A\to H,t}{\mathit{𝝀}}_{age}$. The total force of infection on each individual was the sum of her infectious contact rates within the household and within the community, $\lambda ={\lambda }_{hh}+{\lambda }_{com}+{\lambda }_{ext}$. Where ${\lambda }_{ext}=ϵ\beta (t)/N(t)$ was the force of infection from outside KHDSS.

The actual infection rate for each individual was the force of infection ‘felt’ by the individual times the susceptibility of the individual. The susceptibility of under-one year olds (${\sigma }_{U1}$) depended on whether or not the U1 individual was still protected from RSV by maternally acquired antibodies, which we modelled as giving a random $M$ days of protection; that is for an individual of age $A$ days, ${\sigma }_{U1}=0$ if $M>A$ and ${\sigma }_{U1}=1$ otherwise. In general, the infection status of an individual correlates with her age. However, because RSV is strongly seasonal we do not treat the age of an U1 as correlated with her susceptibility arguing that every U1 is facing her first RSV season irrespective of whether she is 1-month old or 11 months old. Therefore, the mean susceptibility for under-ones was ${\overline{\sigma }}_{U1}=\mathbb{P}(M\le A)$. The susceptibility of over-one year olds was chosen as if the individual had definitely received at least one RSV infection in the past, and definitely had no chance of being maternally protected. We modelled the duration of maternal protection $M$ as a truncated exponential distribution conditioned on being less than 1 year in duration; that is ${\displaystyle M\sim \exp (\alpha )|(M\le 1\text{ year})}$ (appendix 2).

The chance of an infected individual becoming severely diseased after contracting RSV, and then seeking care at hospital, depended on that person’s age and number of infections (Nokes et al., 2008; Ohuma et al., 2012). When an U1 was infected in the model her age at infection was given by conditioning on the age of the U1 being greater than her maternal protection period,

Which was calculated exactly (see appendices 2 and 4). This took into account that increasing the duration of maternal protection would increase the age at infection and therefore reduce the risk of disease. O1s were assumed to have no maternal protection but their conditional age depended on their household type [Equation (2)]. We used these conditional distributions to convert the incidence rate of U1s and O1s in each household type into dynamic incidence rates in each age category, ${ℐ}_a(t)$. By assuming that all O1s had been infected at least once we could use previously published age-dependent hospitalisation odds per infection $h_a$ (Kinyanjui et al., 2015 and appendix 3) to determine the cumulative hospitalisations predicted by the model for each age category $a$ and week interval $w_i=(t_{i,1},t_{i,2})$,

The majority of the parameters for the RSV transmission model were drawn from the RSV literature (see Table 2 and appendix 3) leaving four parameters, and the five hyperparameters of a normal distribution describing the random yearly variation in log-seasonality, to be inferred from hospitalisation data (see Table 3 for parameter estimates and appendix 1 for further details on seasonality model). The free parameters and distribution of the RSV transmission model were:

1. Community infection rate outside the household between U1s and all others in the community accessing KCH ($b_{U1}$).

2. Community infection rate outside the household among all O1s in community ($b_{O1}$).

3. Infectious contact rate within the household to all other household members ($\tau$).

4. Mean duration of maternally derived immunity to RSV ($M$).

5. The joint normal distribution of the yearly log-seasonality amplitude and phase ($[\xi ,\varphi ]\sim 𝒩(\mathit{𝝁},\mathbf{𝚺}$)).

We also included an infectious contact rate for children of schooling age (5–18 years old; $b_S$) which acted additionally to $b_{O1}$; that is children of schooling age were at additional risk of contracting RSV on top of the risk due to mixing in the community. This meant that the mixing matrix in Equation (5) was in block form,

where the blocks represented respectively under-one age categories, over-ones at school age categories and over-ones above school age categories. Unfortunately, we were unable to reliably identify $b_S$ parameter jointly with the other parameters. Investigating a range of $b_S$ values gave similar results for model fit and predictions for vaccine efficacy, the results in the main paper were for the highest value of $b_S$ considered which was mildly pessimistic compared to $b_S=0$ (see appendix 3).

The data for parameter inference was RSV-confirmed, age-specific weekly admissions to Kilifi county hospital (KCH) hospitalisation data from September 2001 until September 2016 (see Nokes et al., 2009 for study details). KCH serves as the primary care facility for the KHDSS population, and we assumed that all KHDSS members who accessed urgent hospital treatment due to RSV disease accessed their treatment at KCH. However, a significant number of admissions were from people not within the KHDSS survey leading to data re-scaling (see appendix 3). The log-likelihood for a particular simulation corresponded to Poisson errors,

where ${𝒟}_{i,a}$ was the cumulative number of hospitalisation observed at KCH in age category $a$ on week $w_i$ and $f_{poi}(x|\mu )$ is the probability mass function for a Poisson distribution with mean µ.

If the yearly realisations of the random seasonality (see appendix 1) were known, then the entire model would be deterministic and $\ln ℒ$ would be a function of the unknown parameters. Therefore, we treated the yearly variation in seasonality as missing data and used the Expectation-maximisation (EM) algorithm (Dempster et al., 1977) to converge onto maximum likelihood estimates for the four free parameters, and the two hyperparameters of the log-seasonality model, 95% confidence intervals were constructed using the likelihood profile technique (e.g. King et al., 2008 and appendix 3).

There were two vaccines used in this modelling study, which were deployed as part of the prenatal contact between pregnant women and skilled health professionals. We assumed that the maternal vaccine was delivered as one injection to the pregnant women in her third trimester. This achieved some unknown additional period of maternal protection, $P$ days, on top of the random period $M$, that is after maternally vaccinating the period of protection became $M_{vac}=M+P$. Achieving an effective maternal vaccination coverage of $V_{cov}$ shifted the mean susceptibility of U1s to ${\overline{\sigma }}_{U1}=\mathbb{P}(M_{vac}<A)V_{cov}+\mathbb{P}(M<A)(1-V_{cov})$, a linear increase in $V_{cov}$. The change in distribution of age at infection was non-linear in $V_{cov}$ because, conditional on an U1 being infected, it was more likely that the U1’s mother had not been vaccinated than the unconditional probability of non-vaccination, $1-V_{cov}$ (see appendix 4). We also assumed that there was a vaccine available that provoked an immune response in the vaccinated individuals similar to a natural infection; that is a susceptible $O1$ who is vaccinated immediately becomes ‘recovered’ and immune to RSV infection until her immunity waned. Immune response provoking vaccination was offered to all O1s in households when a birth occurred, as an addendum to the prenatal contact between mothers-to-be and health professionals. In principle, there were three dimensions to the coverage of the immunity provoking vaccine: (i) coverage of households, (ii) coverage within households, and (iii) vaccine effectiveness. For simplicity, we bundled these dimensions together, and vaccinated whole households at an effective vaccination coverage (the product of the three dimensions of coverage). Over 10 years of forecasted RSV epidemics if a MAB vaccine was available, and given to every pregnant mother, 8601 MAB vaccines were deployed each year. 0–24,095 IRP vaccines were deployed each year depending on household coverage. It should be noted that by 2016 the KHDSS population was around 240,000 people, hence 100% effective coverage of the households where births occurred corresponded to ~10% effective coverage of the total population.

We simulated the model by numerically solving the high dimensional ODE [Equation (3)] simultaneously with the ongoing cumulative hospitalisations in each age category, ${\dot{\mathscr{H}}}_a=h_a{ℐ}_a(t)$, which allowed us to solve for the model predicted weekly hospitalisations [Equation 7]. The initial state of the model was unknown. We initialised the model by starting with a completely susceptible population with the population demography set to mimic that of the KHDSS on 1st Jan 2000. We then simulated RSV transmission for 10 years, with demographic rates (e.g. birth rates) chosen to match those of KHDSS in year 2000 and the seasonal amplitude and phase of $\ln \beta$ set to their latest mean estimate, in order to provide an initial state of the household model. Finally, we ran the model from 1st Jan 2000 until 1st September 2001. This provided the initial point for comparison to hospitalisation data. Numerical solutions were provided using the Sundials CVODE solver (Cohen et al., 1996) implemented within the DifferentialEquations package for Julia 0.6 (Rackauckas and Nie, 2017). For retrospective simulations comparing model predictions to data (Figure 2), we used the most probable values of the yearly seasonality. For forecast simulations, we generated 500 realisations of yearly seasonality over 10 years from the distribution inferred in model inference, this gave 500 predictions for the time series of future hospitalisations. We typically presented medians of these predictions (e.g. Figure 4). The code for the RSV household model used in this paper, and the data used for parameter inference, is available from https://github.com/SamuelBrand1/RSVHouseholdModel (Brand, 2020; copy archived at https://github.com/elifesciences-publications/RSVHouseholdModel).

RSV is a seasonal virus, in temperate climates the peak month for RSV incidence tends to be consistent year-on-year. Therefore, modelling approaches aimed at understanding RSV transmission in temperate climates have used an annually periodic deterministic function, with the timing of peak infectiousness of RSV being either a model parameter (Yamin et al., 2016) or itself a function of climatic variable to be fitted using regression methods (Pitzer et al., 2015).

The seasonal drivers of RSV transmission in the tropics are less clear (Paynter, 2015). At KCH the most common trough month for RSV hospitalisations was September, which lead us to define the RSV ‘year’ as September - September. The most common month for peak hospitalisation in each RSV year was January, however there was significant variation in peak month between RSV seasons with peaks occurring in each month November - April between 2002 and 2016 (Appendix 1—figure 1).

Appendix 1—figure 1

Download asset Open asset

The year-on-year variation in peak month for RSV hospitalisation means that naively inferring a single fixed peak infectiousness parameter would not be a successful inference strategy. However, determining the precise mechanistic reason for shifting seasonality was challenging for the KHDSS population. RSV has been positively associated with the rainy season in some tropical settings (Paynter et al., 2013; Paynter, 2015); however, this is not obviously the case in Kilifi county where the rainy season is April to June with short rains October to December. There have been many proposed mechanisms for erratic periodicity in transmission (for a wide variety of infectious pathogens) which could be relevant to RSV transmission in Kilifi, for example, dynamical attractor switching (Keeling et al., 2001), or the effect of species/strain interaction (Bhattacharyya et al., 2018). In particular, strain competition between RSV A and RSV B has been identified a mechanism for generating complex seasonal dynamics (White et al., 2005).

In this paper, we took an agnostic view and rather than choosing a mechanistic hypothesis for erratic seasonality from the many possible, we assume that the time-varying infectiousness of RSV alters randomly (but from a common distribution) year to year:

where the RSV infectiousness (${\xi }_n$) and seasonal peak timing (${\varphi }_n$) for each RSV year $n$ are drawn jointly from a normal distribution common to each year $({\xi }_n,{\varphi }_n)\sim 𝒩(\mathit{𝝁},\mathbf{𝚺})$. During model inference the yearly ${\xi }_n$ and ${\varphi }_n$ realisations are treated as latent variables; their mean and covariance matrix are imputed along with other model parameters.

As described briefly in the main text, we developed a dynamic model for simulating the spread of RSV through the KHDSS population. The model was a hybrid between a mechanistic ODE approach, this included detailed household structure but only a simplified set of age-and-disease states for individuals within the households, and a data-driven empirical model, this used the observed joint distributions of KHDSS individuals’ household occupancy and ages to generate conditional predications of individual detail beyond that of the mechanistic part of the model.

A commonly used conceptual framework for modelling epidemic transmission with a population is the compartmental model (Anderson and May, 1991; Keeling and Rohani, 2008); each person’s disease state is described as being one of a finite number of possibilities, for example susceptible, infectious, recovered, which define that person’s risk of contracting the infectious pathogen or transmissibility whilst infected with the pathogen. Additionally, it is usually important to capture the heterogeneity of the population, also called the population structure, in contrast to unstructured populations where every individual is treated as interchangeable. Therefore, each person will be described by their position in the population with sufficient detail that a rate of contact can be modelled between any pairs of individuals, see Diekmann and Heesterbeek for a more detailed discussion on modelling population structure (Diekmann and Heesterbeek, 2000). RSV transmission models have most commonly used age structure to describe heterogeneity in the population; each individual is described jointly by their disease state and which age interval (from some predetermined set of intervals) they occupy (Pitzer et al., 2015; Kinyanjui et al., 2015; Yamin et al., 2016). For age-structured RSV transmission models, there are two dynamical elements: the transmission of disease and the demographic turnover of the population (births, deaths and ageing). At the level of the individual these are modelled as discrete random events occurring at some per-capita rate (Rock et al., 2014). However, for large populations, there will be a very large number of individuals in each age-and-disease state, and the flux of population density in each age-and-disease state converges in probability onto the solution of a set of ordinary differential equations (ODEs) as the population size is treated as converging to infinite size (Kurtz, 1970; Kurtz, 1971; Diekmann and Heesterbeek, 2000). The limiting ODE model has as many degrees of freedom as there are age-and-disease state combinations in the epidemic model. In most epidemic modelling studies, it is the deterministic evolution of the solution to these ODEs that is usually given as the transmission model description.

In this paper, the essential modelling concept was to shift the focus away from numbers of individuals in each age-and-disease state and towards the number of households in each possible household configuration. A household configuration describes the number of individuals in each age-and-disease state who cohabit within a single household. Including households within the model adds a potentially relevant layer of realism; the social contacts within a household are persistent, therefore pairs of individuals that cohabit will repeatedly have the opportunity to infect one another if RSV enters the household but be relatively cocooned from infection if RSV has not entered the household. Age-structured transmission models implicitly assume that no two individuals contact one another more than once. To see this consider a population size of $N$; the rate of any individual contacting another single individual is $𝒪(1/N)$ therefore the probability that an individual selects the same other individual twice for contact over any finite time horizon goes to zero as $N\to \mathrm{\infty }$ (which is also the limit at which the ODE model is valid). For household models the discrete random events that change the state of individuals (infection, death etc.) also change the household configuration. When the number of households is very large, there will be a large number of households in each possible household configuration and, as with age-structured models, there is convergence onto a set of ODEs with as many degrees of freedom as the number of possible household configurations.

The possible household configurations, or state space, of a household- and age-structured RSV transmission model is considerably larger than it would be for the equivalent age-structured model. If there are $m$ possible age-and-disease states then the number of possible household configurations for a household of size $n$ is given by a standard combinatorial identity, $\left(\genfrac{}{}{0pt}{}{n+m-1}{n}\right)$. In this paper, we consider a range of household sizes up to a maximum size $n_{max}$, therefore the number of household configurations was,

The number of possible household configurations grows very rapidly (Appendix 2—figure 1). Therefore, having a sufficiently large $n_{max}$ to capture the target population required using a relatively simple compartmental age-and-disease state model for RSV infection.

Appendix 2—figure 1

Download asset Open asset

The age-and-household structured model we used in the main paper to make predictions of potential vaccine effectiveness in a population with persistent social structure. However, it can be useful to compare comparatively complex simulation studies to simpler models which are at least partially analytically tractable; this comparison identifies which features of a model are generic as opposed to emerging from more complicated factors (like seasonality or social structure).

A simple unstructured compartmental model of RSV transmission with two types of vaccine in a population of size $N$ was presented in the main paper (Box 1). Individuals are born into the population at rate $B$ and are initially protected against RSV by maternal antibodies ($M$). All individuals die at rate µ. They lose maternal protection at rate ${\alpha }_{vac}$ (the rate associated with the maternal vaccine) and become susceptible to RSV infection ($S$). Each susceptible is infected at a rate $\beta I/N$ where $\beta$ is the product of the contact rate and the probability of transmission per contact and $I$ is the number of infected individuals in the population. Infected individuals clear their infection and become recovered and are temporarily immune to reinfection ($R$) at rate $\gamma$. Recovered individuals lose their temporary immunity to reinfection at rate $\nu$. A vaccine aimed at provoking an immune response akin to a natural infection (IRP vaccine) is also used to control RSV. This is given to individuals in the population at effective rate $V$ (rate of delivery times probability the vaccine dose is successful). For simplicity, we assume that the IRP vaccine is not given to children so young they are likely to be in the $M$-compartment, but their isn’t memory of which individuals have been vaccinated recently, therefore the chance that an individual selected for vaccination is actually susceptible is $S/(S+I+R)$. If a susceptible individual is vaccinated she transitions to becoming temporarily immune to RSV, this temporary immunity being lost at rate $\nu$.

The ODE equations for the dynamics of the basic unstructured model are:

We solve for the equilibrium state of this simple model, denoted $(M^{*},S^{*},I^{*},R^{*})$, assuming that the population has reached a steady size of $N$, with replacement birth rate $B=\mu N$. For the simple RSV model, we use a mortality rate µ that corresponds to a life expectancy of 65 years, the Kenyan average. The reproductive ratio for the model is $R_0=\beta /(\gamma +\mu )$.

Since, the rate of loss of maternal immunity is fast compared to the mortality $({\alpha }_{vac}\gg \mu )$ nearly all the population survive their $M$ period and become available for infection,

We use $S^{*}+I^{*}+R^{*}=N$ below to simplify the notation, but $N$ could be replaced with $N_{eff}=\frac{{\alpha }_{vac}}{{\alpha }_{vac}+\mu }N$. Note that the maternal vaccine does not alter the incidence rate for the simple RSV model at equilibrium, it simply delays the typical infection time. Equation (13) implies that either $I^{*}=0$ (disease free state), or,

Therefore,

Combining Equations (12), (15), (16), (17) gives that if RSV is endemic then,

Equation (18) implies that for the simple RSV model the critical rate at which an IRP vaccine eliminates RSV is $V_c=(\mu +\nu )N(R_0-1)$.

At an endemic equilibrium, the RSV incidence rate with vaccination rate $V$, denoted ${\iota }_V^{*}$, is therefore,

Equation (19) implies the two results which are presented in Box 1 of the main text:

• The relative reduction in incidence due to IRP vaccination compared to no vaccination is,

  • (20) $\frac{{\iota }_0^{*}-{\iota }_V^{*}}{{\iota }_0^{*}}=\min \{\frac{V}{N(\mu +\nu )(R_0-1)},1\}.$

  • In this paper, we model a scenario where co-habitants of newborn children each receive an IRP vaccine. This fixes $V$ to be proportional to the birth rate, $V=\mu N⟨H⟩V_{cov}$, where $⟨H⟩$ is the average number of co-habitants that a newborn has and $V_{cov}$ is the effective IRP coverage of households. This gives,

  • (21) ${\displaystyle \mathrm{R}\mathrm{e}\mathrm{l}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{v}\mathrm{e}\mathrm{r}\mathrm{e}\mathrm{d}\mathrm{u}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}\mathrm{m}\mathrm{i}\mathrm{s}\mathrm{s}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{d}\mathrm{u}\mathrm{e}\mathrm{t}\mathrm{o}\mathrm{v}\mathrm{a}\mathrm{c}\mathrm{c}\mathrm{i}\mathrm{n}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}=min\{\frac{\mu ⟨H⟩V_{cov}}{(\mu +\nu )(R_0-1)},1\}.}$

• Whilst RSV is not eliminated the reduction in incidence rate due to IRP vaccination is linear in $V$, with the improvement per extra vaccine used being a constant

  • (22) ${\displaystyle \mathrm{R}\mathrm{e}\mathrm{d}\mathrm{u}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}\mathrm{m}\mathrm{i}\mathrm{s}\mathrm{s}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{I}\mathrm{R}\mathrm{P}\mathrm{v}\mathrm{a}\mathrm{c}\mathrm{c}\mathrm{i}\mathrm{n}\mathrm{e}=\frac{(\gamma +\mu )}{(\gamma +\mu +\nu )R_0}.}$

The mean number of over-one year olds living in households with at least one under-one year old in the KHDSS (see below) fluctuated yearly, but was never greater than five ($⟨H⟩<5$). Therefore, using a reversion to susceptibility rate $\nu =2$ per year (see Table 2) with Equation (21) suggests that if, say, $R_0=2$ then the maximum achievable relative reduction in RSV incidence using this strategy with a Kilifi like population implied by the simple RSV model is 3.8%.

A literature review of mechanistic RSV transmission models revealed a number of critical common features:

• At birth newborns are born protected against RSV infection due to antibodies gained from their mother via trans-placental transfer. This is typically modelled as a maternally protected disease state $M$ e.g. (Yamin et al., 2016).

• The probability of developing severe disease and being hospitalised depends on a person’s age, and number of times infected in the past, e.g. (Kinyanjui et al., 2015).

• The susceptibility to RSV infection per infectious contact, their infectiousness after infection, and the expected time taken to become recovered from RSV depend on number of times previously infected, e.g. (Kinyanjui et al., 2015).

The high dimensionality of household- and age-structured models necessitated using the most minimal age-and-disease state model possible for RSV (see above). To do this we use an extremely parsimonious approach. The possible age-and-disease state for individuals are: susceptible or maternally protected and under the age of one ($S_1$), infectious and under the age of one ($I_1$), recovered and under the age of one ($R_1$), susceptible and over the age of one ($S_2$), infectious and over the age of one ($I_2$) and recovered and over the age of one ($R_2$). An under-one year old (U1) experiencing some force of infection $\lambda$ becomes infected ($S_1\to I_1$) and infectious to RSV at a rate ${\sigma }_{U1}\lambda$ where ${\sigma }_{U1}$ is the average susceptibility of an U1 year old to RSV. After becoming infected the U1 ceases to become infectious at a rate ${\gamma }_1$ ($I_1\to R_1$) and then is immune to reinfection to RSV for a period of time. The immunity derived from natural infection is lost at a rate $\nu$, and the U1 revert to susceptibility but in the $S_2$ category ($R_1\to S_2$). The reason we transition recovered U1s to a susceptible over-one year old (O1) is that due to the seasonality of RSV it is very rare for a person to be infected more than once in one epidemic season, therefore functionally by the time an individual is facing the risk of their second RSV lifetime infection they will very likely be over one. All U1s age at the rate $\eta =1/365.25$ days_-1 becoming individuals in the same disease state but over-one ($S_1\to S_2$, $I_1\to I_2$, $R_1\to R_2$). An O1 individual experiencing a force of infection $\lambda$ becomes infected and infectious ($S_2\to I_2$) with RSV at a rate ${\sigma }_{O1}\lambda$ where ${\sigma }_{O1}$ is the relative susceptibility of O1s compared to an U1 no longer protected by maternal antibodies. Infectious O1s cease being infectious ($I_2\to R_2$) at a faster rate than U1s, ${\gamma }_2>{\gamma }_1$, but revert to susceptibility ($R_2\to S_2$) at the same rate $\nu$ (Appendix 2—figure 2).

As mentioned in the main document we relate this simple age-and-disease state model to more complicated RSV models by (i) using the conditional age distribution of individuals to address questions that required a more complicated age structure than a simple under/over-one binary choice, for example whether susceptible under ones were still protected by maternal antibodies, and (ii) by assuming that all over-ones have been infected at least once and all susceptible U1s have never been infected and might still be protected by maternal antibodies.

Appendix 2—figure 2

Download asset Open asset

A household configuration is a tuple of the number of individuals in each age-and-disease state who cohabit a household. The generic household configuration is denoted $h=(s_1,i_1,r_1,s_2,i_2,r_2)$, indicating that the household has precisely $s_1$ individuals in state $S_1$, $i_1$ individuals in state $I_1$ etc. The household size is the number of people living in the household (i.e. $s_1+i_1+r_1+s_2+i_2+r_2$). We denote the space of possible household configurations $\mathrm{\Sigma }$ and number of households in the state $h$ at time $t$ as $H_h(t)$. It is useful to consider a vector quantity over all possible household configurations such as $\mathit{𝑯}(t)=(H_h(t)|h\in \mathrm{\Sigma })$ where we have generated some ordering for elements $h\in \mathrm{\Sigma }$. It is clear that the knowledge of $(\mathit{𝑯}(t),t\ge 0)$ would allow us to reconstruct the dynamics of individuals. For example, using the function $f(h)=s_1$ for each $h\in \mathrm{\Sigma }$ in a vectorised form $\mathit{𝒇}=(f(h)|h\in \mathrm{\Sigma })$ allows us to track the dynamics of numbers of individuals: ${\displaystyle (f\cdot \mathbf{H}(t),t\ge 0)}$.

As mentioned above, age-structured models are constructed by considering the per capita rate of events affecting the state of individuals. Household- and age-structured models are constructed by considering the per household rate of events that affect the household configuration (see House and Keeling, 2009 for further mathematical details). In the following we list the events that change the household model divided into three groups: events due to transmission within the household, events due to transmission between households and events due to demographic turnover.

The model above could in principle have an infinite number of states if the household size was not limited (see above). We chose limits on the household size based on capturing ≈99% of the U1s in the population, and therefore the pathway to them catching RSV. The limits were: (i) no household is bigger than size 10, and (ii) no household has more than 2 U1s. This also covers the big majority of the total numbers of households (see Appendix 2—figure 3). The $n_{max}=10$ limit was imposed by initialising the model without households of size $>10$, and setting $r(n,t)=0$ for all $n>10$. The ≤2 U1 limit was imposed by setting the birth/turnover rate to zero for all households with 2 U1s. Putting the limits in reduces the dimensionality of the system to 1926 different household configurations.

Appendix 2—figure 3

Download asset Open asset

Note that the events that either change a household’s configuration or change the number of households described above can be divided into two categories: (Nair et al., 2010) those with rates that only depended on the household’s configuration, e.g. infection within the household, or ageing of U1s, and, (Glezen et al., 1986) those with rates that depended on the configurations of other households, e.g. transmission between households or the rate of change of household numbers. The events in category (Nair et al., 2010) translate to linear dynamics for $\mathit{𝑯}(t)$, events in category (Glezen et al., 1986) translate to non-linear dynamics (House and Keeling, 2008). Overall, the dynamics of $\mathit{𝑯}(t)$ obey the semi-linear dynamical system,

A_t is a matrix which encodes the dynamics of events in category (Nair et al., 2010), ${\mathit{𝒇}}_t(\mathit{𝑯}(t))$ encodes the transmission between households, and ${\mathit{𝝆}}_t(\mathit{𝑯}(t))$ encodes the rate of change of numbers of households in each configuration. We initialised the dynamics of Equation (51) by starting with a completely susceptible population on 1 st Jan 1990, allowing RSV to be introduced via the external force of infection and running for 10 years (see main text).

Equation (51) has two properties that are important to note:

• The change rate in households of size $n$ is independent of the transmission dynamics:

  • (52) ${\partial }_t\left({\displaystyle \sum _{h\in {\mathrm{\Sigma }}_n}}{H}_h(t)\right)=r(n,t),n=1,\mathrm{\dots },10.$

• The dynamics of the proportion of households in a given state $P_h(t)=H_h(t)/{\sum }_{h^{\prime }}{H}_{h^{\prime }}(t)$ is not directly affected by the change rates (${\mathit{𝝆}}_t$) in households:

  • (53) ${\partial }_t{\mathit{𝑷}}_t=A_t{\mathit{𝑷}}_t+{\displaystyle \frac{{\mathit{𝒇}}_t\left({\mathit{𝑯}}_t\right)}{{\sum }_{h^{\prime }}{H}_{h^{\prime }}(t)}}$

Equations (52) and (53) guarantee the desired modelling features discussed above. Equation (52) gives that the change in the number of households of each size matches the empirical rate of change for each year, we also verified this by numerical solution of Equation (51) (Appendix 2—figure 4). Equation (53) shows that the rate of change of household numbers doesn’t directly effect the proportion of households in any given configuration. We also verified that the number of U1s and O1s was close to their empirical values (Appendix 2—figure 5).

Equation (51) was difficult to solve efficiently because it is both numerically stiff and high dimensional. We numerically solved Equation (51) using the Julia DifferentialEquations package implementation of the CVODE solver, with an efficient Krylov method (GMRES) to solve the implicit timestepping (see main text). We also used the DifferentialEquations efficient event handling which allowed us to change parameters (like the household change rate) at specific times without damaging the performance of the solver, or having to restart simulations.

Appendix 2—figure 4

Download asset Open asset

Appendix 2—figure 5

Download asset Open asset

The parameters for the household- and age-structured transmission model were drawn from four sources:

• A literature review of infectiousness duration and other epidemiological quantities; main Table 2.

• Calculated from the empirical joint distributions (see above Appendix 3—table 1).

• Age-dependent hospitalisation probability per RSV infection derived from Kinyanjui et al., 2015; Appendix 3—table 2. Hospitalisation probability was the probability that an infected individual would develop severe disease, multiplied by the probability that severely diseased individuals would require hospitalisation. The probability that an infected individual became diseased depended on whether it was the individual’s primary infection episode or not. The underlying data for estimating these probabilities was drawn from cohort studies on RSV disease rates (Ohuma et al., 2012; Nokes et al., 2008). We adapted these probabilities for our model using our assumption that all infected under-ones were experiencing their first RSV episode, and all over-ones were experiencing their second or subsequent infection.

• Inferred from the KCH hospitalisation data set (see below).

As mentioned in the main text we used the EM algorithm (Dempster et al., 1977) to estimate parameters for the model. Again, as described in the main text the parameters we chose for inference were:

• Infectious contact rate outside the household between U1s and all others in the community accessing KCH ($b_{U1}$).

• Infectious contact rate outside the household among all O1s in community ($b_{O1}$).

• Infectious contact rate within the household ($\tau$).

• Rate of loss of maternally derived immunity to RSV ($\alpha$).

• The joint normal distribution of the yearly log-seasonality amplitude and phase ($[\xi ,\varphi ]\sim 𝒩(\mathit{𝝁},\mathbf{𝚺}$)).

where the community age mixing matrix $T(a,b)$ was in block form:

The log-likelihood for our model [Equation (8) main text] was defined using the incidence rates ${ℐ}_a(t)$ predicted by solving the model. The incidence rate for all the households in the generic household configuration was,

where the household force of infection for the generic household configuration was ${\lambda }_{hh}=\tau (i_1+{\iota }_2{i}_2)$. We converted the household incidence rate into an age structured incidence rate by using conditional age distributions, and this allowed us to calculate the cumulative hospitalisations in age category $a$, predicted by a given set of parameters and yearly seasonality realisations, in weekly intervals $w_i=(t_{i,1},t_{i,2})$ using the age-dependent hospitalisation rates per infection $h_a$ (see Table 2)