Multi-Season Analysis Reveals the Spatial Structure of Disease Spread

TL;DR

This study introduces a tensor-driven multi-compartment SIR model to analyze and predict the spatial and seasonal dynamics of disease spread using internet search data, revealing key transmission patterns.

Contribution

The paper develops a novel tensor-based multi-compartment SIR model and applies it to internet data to uncover spatial structure and predict disease spread for seasonal viruses.

Findings

Human movement explains 27-30% of inter-state transmission.

Model fits RSV and WNV data with R^2 of 0.70 and 0.52.

Disease load patterns can be predicted with minimal current season data.

Abstract

Understanding the dynamics of infectious disease spread in a heterogeneous population is an important factor in designing control strategies. Here, we develop a novel tensor-driven multi-compartment version of the classic Susceptible-Infected-Recovered (SIR) model and apply it to Internet data to reveal information about the complex spatial structure of disease spread. The model is used to analyze state-level Google search data from the US pertaining to two viruses, Respiratory Syncytial Virus (RSV), and West Nile Virus (WNV). We fit the data with correlations of $R^2=0.70$, and $0.52$ for RSV and WNV, respectively. Although no prior assumptions on spatial structure are made, human movement patterns in the US explain 27-30\% of the estimated inter-state transmission rates. The transmission rates within states are correlated with known demographic indicators, such as population density and average age. Finally, we show that the patterns of disease load for subsequent seasons can be predicted using the model parameters estimated for previous seasons and as few as $7$ weeks of data from the current season. Our results are applicable to other countries and similar viruses, allowing the identification of disease spread parameters and prediction of disease load for seasonal viruses earlier in season.

Equations

$$\frac{dS}{dt}=-\beta SI,\frac{dI}{dt}=\beta SI-\gamma I,\frac{dR}{dt}=\gamma I,$$

$$\begin{array}{c}\begin{array}{rl}\frac{d\mathbf{S}}{dt}& =-\mathbf{I}{\mathbf{\beta }}^{⊺}\mathbf{S}\\ \frac{d\mathbf{I}}{dt}& =\mathbf{S}{\mathbf{\beta }}^{⊺}\mathbf{I}-{\mathbf{\gamma }}^{⊺}\mathbf{I}\\ \frac{d\mathbf{R}}{dt}& ={\mathbf{\gamma }}^{⊺}\mathbf{I},\end{array}\end{array}$$

Full text

1 Introduction

The spread of infectious diseases through a population of susceptible individuals varies in time, space, and according to the characteristics of the susceptible individuals and the disease[1]. Modeling this spread can provide information on the spreading mechanism of the infection and can assist in designing strategies for control of the disease. One common method to model infections is compartmental models[2, 3], which in their basic form assume a single population divided into 3 sub-groups of susceptible, infected, or recovered (SIR) individuals. This model describes a global spreading mechanism where each infected individual can infect any susceptible individual. Another common model describes local spreading[4], where infected individuals only spread infection to a limited subset of susceptible individuals.

In reality, many epidemics spread through a combination of two or more spreading mechanisms; hybrid spreading[5]. Such spreading can be modeled by multiple compartments encompassing different demographic characteristics[4], various disease strains[6, 7, 8, 9, 10], and cross-immunity in an age structured model[11, 12, 13]. However, fitting data from a given epidemic to these models remains a challenge owing to a dearth of data of sufficiently high spatial and temporal resolution. This problem is exacerbated the higher the resolution of the model, since the number of model parameters grows with the number of compartments of the model. Thus, in many cases researchers introduce simplifying assumptions in order to reduce the dimensionality of the models[14, 15, 16]. In the absence of actual data, the spreading topology can generally be analyzed by numerical simulation of models using synthetic data, for example, of the patterns of human movements at the desired scale[16, 17]. It can also be analyzed through the use of theoretical tools to detect critical phenomena in networks[18, 19, 5].

Another approach to overcome a lack of epidemiological data is to use proxy data, such as mobile phone data[20] or Internet data[21]. Internet data, including search engine queries, social media postings[22, 23], or Wikipedia visits can serve as a proxy for such information and provide data at high spatial and temporal resolutions from very large populations. One demonstration of the value of these data is the recent evaluation of the childhood influenza vaccine campaign in England[23], which showed that vaccination of primary school-age children significantly reduces influenza rates as estimated from Twitter posts and Bing queries.

In this study, we focus on modeling the spatial and temporal spread of infectious disease using a tensor-driven multi-compartment version of the classic Susceptible-Infected-Recovered (SIR) model. We illustrate the ability of the model to capture disease parameters by fitting it to Internet data on two common viruses, Respiratory Syncytial Virus (RSV) and West Nile Virus (WNV) in the United States (US). Our results demonstrate that the spatial and temporal dynamics of these viruses can be captured by this model with almost no a-priori assumptions on the spreading topology and interaction. This allow us to characterize the spread of the virus of future seasons. The inferred model parameters are shown to be correlated with population indices such as human movement dynamics, and with demographic and environmental factors which are likely modulating transmission.

2 The model

The proposed model is a generalization of the single population SIR model[2], that describes the average evolution of an infectious disease in a single population. In the SIR model, disease spread is represented by a system of ordinary differential equations for the susceptible $S$, infected $I$, and recovered $R$ individuals. It accounts for only one pathogen and assumes equal probability of being infected by any individuals regardless of his position i.e. mean field interaction.

The model equations are:

[TABLE]

where $\beta>0$ is the infection rate and $\gamma>0$ is the recovery rate. This model, Eq. (1), can be generalized to account for multiple viruses and spatial effects due to groups of sub-populations and virus strains, by transforming to a multidimensional representation in which $S$, $I$, $R$, $\beta$, and $\gamma$ are tensors (Eq. (2)). These tensors act as time-dependent state representation of the entire system, reflecting its various attributes or dimensions. That is, each two-dimensional projection represents different sub-groups e.g. viruses, spatial regions. Here, the system describes the evolution of an infectious sub-population interacting with other sub-populations under the influence of multiple viruses. The dynamics for the system’s state tensors evolution is represented by the following system of ordinary differential equations:

[TABLE]

where boldface represents a tensor. We refer to this model as a multi-compartment SIR model (mcSIR). This model, and its reduction to account for multiple viruses and virus strains in one population, is presented by Levy et al.[6]. Thus, the multi-compartment model is represented by a multi-dimensional set of equations. This representation of sub population and viruses can also be generalized to other epidemic models, such as the SIRS and the SIS[1].

In this work, we focus on the spreading mechanism in a group of $N$ sub-populations (different spatial regions) analyzed separately for each virus. Within each sub-population the interaction between individual is well mixed i.e., of mean field type. In this case, $\bm{S}$, $\bm{I}$, and $\bm{R}$ take the form of square matrices with dimension $N$. The specific elements in these matrices, e.g. $S_{ij}(t)$, are defined as the fraction of susceptible individuals traveling from region $i$ to $j$ at time $t$, where $i=j$, represents the fraction of susceptible individuals in sub-population $i$ present at time $t$. The off-diagonal elements of the matrix $\bm{\beta}$ account for the infection rates between sub-populations and the diagonal elements represent the infection rates within each sub-population. Note that, $\bm{\beta}$ is not a symmetric matrix, as the probability to move from $i$ to $j$ is not assumed to be the same as moving from $j$ to $i$. The matrix $\bm{\gamma}$, under the assumption of distinct sub-populations, is defined as a diagonal matrix. Each element, $\gamma_{ii}>0$, is the average recovery rate for individuals belonging to the $i$th sub-populations. For our purposes, we also assume that the number of infected people in each region is the sum of all the infected people traveling to this region, and that the number of susceptible people traveling is negligible compared to the sizes of the susceptible sub-populations in the region. In this case, the model is reduced to a traditional sub-population model[4].

The advantage of using a sub-population model is that in contrast to situations where populations might be well mixed (a common assumption of models without spatial structure) local disease spread is taken into account. This means that most of the population is not exposed to the infection immediately upon its introduction. Rather, the disease may have to pass through many intermediate individuals before reaching all members of the population. This is mainly due to the fact that at a larger scale, individuals that are close together in space are likely to come into contact with each other more frequently than individuals that are more distant. In addition, this locally structured populations may exhibit clique behavior. Moreover, in many cases one can find a spatial epidemic wave pattern depending on the virus type and location[24, 21, 25]. This assumption, which to some extent holds true even today when people move longer distances over shorter times, may smooth-out the clique behavior. Understanding the spatial nature of disease transmission has important consequences for control measures aimed at a disease per location.

In order to validate the use of sub-population models presented in Eq. (2) to describe the observed Internet data, we tested the model on seasonal data derived from search engine queries about two commonly spread viruses in the US, RSV and WNV. In this case, the sub-populations are the 50 states. We chose these two viruses to fit to an SIR model since there is no vaccine or medicine for prevention or treatment. Moreover the symptoms are usually mild meaning people are less likely to visit a doctor or have the case reported to the national health care authorities. On the other hand, Internet search data in these cases may provide real time high resolution estimate for the disease activity.

3 Methods

3.1 Data sources

We used two Internet data sources in our analysis: Google Trends[26] and Twitter. Google Trends was used to determine the relative number (per state in the US) of queries made about each of the two viruses as described below. Data was extracted at a weekly resolution. Twitter was used to estimate human mobility patterns as described below.

3.1.1 Respiratory Syncytial Virus (RSV)

Data about RSV were extracted for the years 2013 to 2018 i.e., 5 seasons. A season starts at the beginning of October every year. The data comprised of the number of queries for the term “RSV”. A single term was sufficient in this case because of the observation[21] that this number is in high correlation with the number of people reported as infected in each state, as published by the Centers for Disease Control and Prevention (CDC). Our data included all US states excluding Hawaii.

3.1.2 West Nile Virus (WNV)

Data about WNV were extracted for the years 2014 to 2018 i.e., 4 seasons. A season starts at the beginning of April every year. Preliminary analysis of Google Trends data showed that the number of queries for the term ”West Nile” from each state during 2017 was highly correlated with the yearly total number of infected cases reported by the CDC, reaching $R^{2}=0.58$ (p-value=$1.4\cdot 10^{-5}$). Thus, we used the number of people querying for ”West Nile” from 2014 to 2018 as a proxy for WNV disease load. Our data included all contiguous US states, excluding Hawaii and Alaska.

3.1.3 Data on human mobility patterns in the US

We collected all messages from Twitter having a GPS location in the US from October 1, 2015 through April 31, 2016. For each message we extracted an anonymized user identifier, the time of the message, and the location from where it was made. These data comprised of approximately 50 million messages. The exact GPS location of each message was mapped as its encompassing US state. From these data we created a matrix of mobility patterns, where the $(i,j)$-th cell comprised of the number of people whose location in one tweet was state $i$ and in their following tweet was state $j$. This matrix was normalized by dividing the number of people moving from state $i$ to state $j$ by the total number of people who moved from state $i$ to any other state.

3.1.4 Demographic information

We extracted from the US census[27] variables that may affect the spreading rate of a virus within a population. These data, extracted per state, included:

Population size (number of people) 2. 2.

Density per square mile of land area. 3. 3.

Children in the following age groups: 0-4, 5-11,12-14, 14-17 (%) 4. 4.

Percentage of people in each of 6 main racial groups 5. 5.

Average income per household and per family (Dollars) 6. 6.

Poverty rate in adults and child under 18 (%) 7. 7.

Senior people (65 years and over) in the population (%) 8. 8.

Average age 9. 9.

Urban population (%) 10. 10.

Unemployment rate (%).

3.2 Parameter matching

Fitting is carried out in two steps. In the first step, a dictionary of functions is generated. Each function is the infection rate function from an instantiation of a solution to the single population SIR model (Eq. (1)) chosen from a wide range of parameter values. Specifically, we used $\gamma\in[10^{-4},0.1]$ with $50$ steps and $\beta\in[10^{-4},1]$ with $200$ steps. Each function entry is intended to be proportional to the number of infected individuals per state over time sampled at a weekly resolution. Out of the dictionary, one function is chosen, such that it has the highest correlation with the data from each state over a period of one season.

In the second step, we fit the data to the multi-compartment SIR (mcSIR) model presented in Eq. (2). The parameters found in the previous step are used as an initialization to the mcSIR fitting process. Specifically, the on-diagonal elements of the matrix $\bm{\beta}$ were initialized to the values found from the single SIR model. The off-diagonal values of $\bm{\beta}$ were initialized to uniformly distributed random values between 0 and the average $\beta$ found from the fit to the single population SIR. We perform a search in the parameter space for the solution of the mcSIR model which maximizes the highest geometric mean across states of the correlations between the mcSIR estimated infection rates and data acquired from Google Trends calculated over four seasons.

The use of multiple seasons provides us with more data to increase robustness to noise, under the assumption that the change in the infection rate from season to season and the recovery rate over the years can be accounted for by a single season-dependence multiplicative scaling factor $D$.

Another assumption made in the fitting of the mcSIR model is the use of a scalar recovery rate $\gamma$ for all seasons and all states. This assumption was made because of the observation from the single SIR fit that the recovery rates are similar for all states and seasons.

The gradient descent optimization process was limited to $15$ iterations. We repeated the optimization process $100$ times, each time with random initialization on the off-diagonal terms of the $\bm{\beta}$ matrix. The run with highest correlation between the data and the infected rate estimated by the mcSIR model was chosen.

The single population SIR and mcSIR models are solved numerically. All analyses were performed using Matlab R2017b[28].

4 Results

The algorithm reached an average correlation (across states) $R^{2}$ of 0.70 for RSV and $R^{2}$ of 0.52 for WNV between the Google Trends data and the mcSIR model.

Figure 1 shows the distribution of the model fit ($R^{2}$) between the Google Trends data and the mcSIR estimated infection rate in each state over $4$ seasons for RSV ($200$ state and season pairs), and $3$ seasons for WNV ($147$ state and season pairs). For both viruses, the majority of states show high correlation between mcSIR model and the Google Trends data. For RSV (Figure 1a) the correlation distribution is centered around $0.8$. Although, the data for WNV (Figure 1b) was noisier, there was still a significant number of states with a correlation higher than $0.6$. We hypothesize that WNV data is noisier than RSV due to the lower incidence of the former.

4.1 Infection rate within states

After fitting the mcSIR model to the data, the infection rate matrix $\bm{\beta}$ and the recovery rate $\gamma$ of each virus could be examined. We note that the best value found for the recovery rate is $\gamma=0.02$ for both RSV and WNV.

The values along the main diagonal of the matrix $\bm{\beta}$ represent the infection rate within each sub-population normalized by the sub-population size[4]. Google Trends data is normalized to the maximal number of queries in each state in a given season. This, the infection rate within each state require similar normalization. Since the number of queries in a state are correlated with its population size, in order to estimate the infection rate in each state we multiply the diagonal elements of $\bm{\beta}$ by the state’s population size.

Then, in order to explain the estimated infection rate values, we built a rank regression model where the independent variables were the known demographic variables of each state (see Methods) and the estimated infection rates in each state.

In the case of RSV, the model fit is $R^{2}=0.42$ (p-value$=9.9\cdot 10^{-7}$). Statistically significant variables are population density (slope: $0.64$) and average age (slope: $-0.35$). The model for WNV reached a fit of $R^{2}=0.49$ (p-value$=2.92\cdot 10^{-7}$). Statistically significant variables are population density (slope: $0.64$), poverty (slope: $0.31$), and average age (slope: $-0.30$).

Infection rates of viruses increase in areas of denser population[29] and so both viruses have a positive correlation with population density. For both viruses there is a negative correlation with average age meaning that younger people are more susceptible. RSV is known to be more severe in infants[30]. WNV is commonly spread by mosquitoes, which may explain the positive correlation with poverty, as mosquitoes could be more prevalent in such areas[31].

4.2 Inter-state infection rates

The off-diagonal elements of the infection rate matrix $\bm{\beta}$ represent the rate of infection between states. That is, these elements represent the rate at which susceptible people in one state are likely to be infected by people from another state. Thus, we hypothesized that these elements should be correlated with human mobility patterns among states. To estimate the correlation between the state-level mobility patterns estimated from Twitter data and the values of $\bm{\beta}$, we normalized $\bm{\beta}$ by the average number of infected people estimated by the mcSIR model in each state during a particular season i.e., the infection rate matrix $\bm{\beta}$ is multiplied by a diagonal matrix of the average number of people transmitting the disease in a season in each state. This quantity is defined as the average transmission of the infection between different states.

The Spearman correlation between the movement patterns, as estimated from Twitter, and the estimate normalized infection rates for RSV was $\rho=0.30$ (p-value $<10^{-10}$) and for WNV was $\rho=0.27$ (p-value $<10^{-10}$). Thus, approximately a significant portion of of the variance in transmission rates of the examined viruses is explained by human mobility patterns.

Embedding the transmission rate matrix on a map of the US provides a clearer understanding of disease spread. Figure 2 shows the transmission rate matrix averaged over the first $25\%$ of the season for the RSV virus in 2016. The color of each state represents the week of peak infection as estimated by the model. Thus, the model predicts Florida would be the first state where the disease peaks. This also identifies Florida as the source of the virus. This is in line with our knowledge of the virus spread in the US[30]. The color of the arrows in Figure 2 represents the strength of the interaction between each two states given by the elements of the transmission matrix $\bm{\beta}$. Since the transmission matrix is correlated with the movements of people, the matrix elements (arrows on the map) represent the beginning of spatial virus spread in the US predicted by the model.

4.3 Characterization of the disease’s spread in the following seasons

Our results suggest that the parameters of disease spread change little from year to year. This can be observed by looking at the season dependent scaling factors, $D$ (defined in the Methods). This factor approaches a value of $1$ when there is seasonal invariance. For RSV the average value of $D$, over $4$ seasons was $0.93\pm 0.06$, and for WNV the average value over $3$ seasons was $0.92\pm 0.08$. Thus, we sought to predict infection rates by applying the model with the estimated parameters from previous seasons, corrected using the first few weeks of data from the current season. This allow us to predict the infection rates of the entire season and the temporal location of the infection peak in each state. The latter is simply the extremal point of the model Eq. (2).

Here, we used the estimates of $\bm{\beta}$ and $\gamma$ from the previews seasons, together with increasing amounts of data from the 2018 season to predict the seasonal dynamics for the entire 2018 season. Specifically, we used the parameters from the previews seasons as initial parameters for the model and adjusted them using the first $N$ weeks of data from the 2018 season, by running a few iterations of the optimization algorithm.

Figures 3a and 3c show the correlation between the actual and predicted infection rates given increasing amounts of data for the current season (increasing $N$), for RSV and WNV, respectively. Figures 3b and 3d show the error in prediction of peak infection for RSV and WNV, respectively. Figure 3 shows that within 7 weeks from the beginning of the season it is possible to predict the progression of the entire season with reasonably high accuracy, both for RSV and WNV.

5 Discussion

Multi-compartment epidemic models provide a useful tool in the understanding of disease spread. However, limits on availability of epidemiological data has made the testing and validation of such models complicated. Here, we demonstrated the utility of matching proxy data from search queries to theoretical models. This matching allows us to confirm the spatial and temporal dynamics of disease spread for two sample viruses, RSV and WNV. Additionally, this matching provides information on the parameters of disease spread. It reveals and emphasizes the importance of accounting for the spatial structural complexity of the spreading mechanism. These parameters were shown to be correlated with the demographics of people in each state, and the mobility of people across states. In addition, we show that learning the parameters of the disease spread from proxy data over multiple seasons generalizes relatively well to the next season.

Our work has drawbacks, especially, in the use of proxy data, rather than epidemiological data. Still, although less accurate, the proxy data are derived at higher spatial resolution and from a much larger, and arguably a more representative, population. The correlation between the parameters we find and the demographics of people provide an indication for the validity of this data. Therefore, this approach can be used to estimate factors influencing disease intervention and control especially in cases where epidemiological data is lacking or insufficient providing more information about the spreading mechanism.

The modeling approach we proposed makes minimal assumptions on the structure of the population, interventions, or behavior. Our approach can be used to track and model other viruses and sub-populations. Future work will focus on modeling sub-populations (for example, different age groups) within a geographic area, as well as the effects of interventions such as vaccination.