Inferring high-resolution human mixing patterns for disease modeling

We report here the results for populations of 277 sub-national administrative regions of Australia, Canada, China, India, Israel, Japan, Russia, South Africa, and the United States of America, characterizing contact patterns for about 3.5 billion individuals. We also include data at the national level for 26 European countries [\cite[cite]{[\@@bibref{}{missing}{}{}]}]FUM12. The inferred age-specific contact matrices reveal striking patterns, of which many are common to the diverse locations under study. Figure  2 shows the age mixing patterns $F_{ij}^{k}$ defined as the per capita frequency of contact of an individual of age $i$ with an individual of age $j$ in setting $k$.

Starting with the households setting in Fig. 2A, we observe that contacts between individuals can largely be characterized as that between couples living together, and parents and their children in the same household [11, 20]. The increased frequency of contact between adults of similar ages along the main diagonal of the household contact matrix represents couples of similar ages living together, while the bands of high frequency above and below the main diagonal indicate contact between parents and children. While most locations share these overall features, the contact matrices show different age-mixing patterns. For instance, in China (Fig. 2A), the lower frequency of contact between children within households is the reflection of the country’s so-called “One-child policy”. The policy, enacted in 1979, has resulted in over a generation of many Chinese youth growing up without siblings, and hence having less contact on average with other children in this setting. This is in stark contrast with the United States and India (Fig. 2A), where the presence of multiple children born to a family results in an increased frequency of contact between this age group in the household matrix. The presence of multigenerational families in countries like India is also evident from the increased frequency of contact between all age groups, notably between the elderly (60 years and older) and young children. Even within the same country, contact patterns may be markedly different. Fig 2A shows the age-mixing patterns within households for two different provinces of China: Beijing and Guizhou. While the household contact patterns in Beijing show a clear signal of the “One-child policy”, Guizhou shows the presence of multigenerational families, as well as an increased presence of multiple children living in the same household. This can be traced back to the fact that the Guizhou Province is characterized by a large frequency of minority groups and the “One-child policy” was less strictly applied for minorities.

Fig 2B,C show the inferred contact matrices in the school and workplace settings for China, the United States, and India. In both settings, the age mixing patterns vary strongly, reflecting differences in the educational systems, and economic conditions unique to each location. For all locations in our study, the school setting consistently exhibits the highest frequencies of contact between children and young adults attending school together. Interaction with older adults in this setting reflects the contact students have with instructors and other staff members in school. The variability of age mixing patterns between children in India (Fig 2B) also reflects the many different kinds of schools that children can attend throughout the country and the different age groups found in those schools. In the workplace environment most interaction takes place between individuals in the range of 20 to 65 years of age, with the age range depending on local retirement, employment regulations, and culture. For instance, in many parts of the world it is common for teenagers to be fully or partially employed (see the work contact matrix for the US - Fig. 2C); in India, census records for employment list even children among the population of workers.

Statistical validation of the contact matrices against summary statistics of a large set of socio-demographic indicators has been performed to validate our results (see SM).

The contact matrices obtained in each settings acquire epidemiological relevance when combined together to generate the effective descriptions of human mixing patterns relevant to the spreading of a specific disease. Here we define the matrix of effective contacts relevant to influenza transmission based on the relative contribution of the household, school, and workplace. In addition to these three social settings, we consider also the contribution of less structured causal encounters in the population [30], by considering a community contact matrix that assumes individuals as potentially fully mixed [20]. To combine the different matrices, we propose a weighted linear combination of the derived matrices for the four considered social settings, and compute the overall matrix of contacts between individuals of age $i$ and individuals of age $j$, $M$ (whose elements are denoted as $M_{ij}$), as a weighted linear combination of setting specific contact matrices:

where the element $M_{ij}$ represents the average number of contacts with individuals of age $j$ for an individual of age $i$ per day, and each $\omega_{k}\geq 0$ is indicating the number of contacts in each setting $k$.

Generally the $\omega_{k}$ are unknown disease-specific weights accounting for the relative importance of the different social settings in the transmission of a specific infectious disease. In the case of airborne infectious diseases we leverage data from diary-based survey contact matrices reported in [11, 15] for Finland, Germany, Italy, Luxembourg, The Netherlands, the United Kingdom, and the Tomsk Oblast region of Russia. Namely we perform a multiple linear regression analysis to find the values of $\omega_{k}$ such that the resulting $M_{ij}$ best fits the empirical data. Note that the empirical matrices derived in [11, 15] describe the average number of contacts of age $j$ for an individual of age $i$, and in the Methods section we show how $\omega_{k}$ is related to average number of contacts $<c>$ per individual. The regression yields 4.11 contacts (standard error, SE 0.41) in the household setting, 11.41 contacts (SE 0.27) in schools, 8.07 contacts (SE 0.52) in workplaces, and 2.79 contacts (SE 0.48) for the general community setting. It is worth remarking that this approach provides overall best matching $\omega_{k}$ and that, in principle, some of the difference of social behavior of specific countries may not be captured by this approach. For this reason, as a validation of this calibration method, in Fig. 3A we report the correlation between the resulting synthetic matrices for France, Japan, and the Shanghai Province of China and the available empirical matrices for these additional locations [26, 27, 28]. We find significant ($p$-value $<0.001$) Pearson correlations of 0.92, 0.9, and 0.8 for France, Japan, and Shanghai Province, respectively. Moreover, we use the Canberra distance as a measure of the similarity between two contact matrices [20] (see Methods section for the definition of the Canberra distance). We estimate the distance between the seven survey-based matrices used in the calibration phase and their respective synthetic matrices to be 0.21 on average (range: 0.17-0.28). (Note that the resulting Canberra distance is normalized by the square of the number of elements of the contact matrix to account for the different number of age groups considered by the different diary-based contact surveys). When considering the three locations used as out-of-sample validation, we estimate a slightly larger average distance of 0.29 (range: 0.21-0.37), suggesting the adequacy of the employed methodology. Finally, Fig. 3B shows a visual comparison between the synthetic and survey matrices, which highlights that the synthetic contact matrices are able to capture the specific features of each location such as contact patterns at school and the relative intensity of the main diagonals.

Figure  4A shows the synthetic overall contact matrices for China, the United States, and India at the country level. The contact matrices for all locations share many similarities: bands of increased contact along the main and off diagonals reflect the familiar household contact patterns, increased contact between adults age 20 to approximately 65 years old account for the interactions between the population’s workforce, and the dominant contact patterns in the lower left of the contact matrices reflect the high number of interactions between school aged individuals. Depending on the age structure of the population the intensity of interactions occurring in the school setting can vary, however this feature consistently dominates the contact matrix for all locations in our study.

To quantify the similarity between the overall contact matrices in different locations, we use a hierarchical clustering algorithm based on the Canberra distance to identify clusters of locations (dis)similar to each other [20]. We find that locations tend to cluster together by country (Fig. 5A), indicating that overall the contact patterns within a single country are more similar to each other than to the patterns observed in other countries. Strikingly, though not surprisingly, locations within developed countries such as Australia, Canada, and the United States are similar to each other and are clustered together, while at the same time locations throughout India, South Africa, and the North Caucasus region of Russia also cluster together, indicating a similarity in patterns between locations in the developing and transition world. Interestingly, a few territories of Canada, Russia, and India are outliers, indicating that the contact patterns in these locations are different from what is observed in all other locations (including their respective countries). A more detailed discussion is reported in the SM. If we consider the US state of New York as reference and compute the distance from all other locations to it, a geographical pattern clearly emerges (Fig. 5B). Indeed, the contact patterns in most states of the US, and the urbanized areas of Canada and Australia appear to be very closely related to the one inferred for New York. In contrast, most of India, South Africa, and of the territories in Canada, Russia, and Australia have contact patterns noticeably different from those obtained for the state of New York.

To investigate the effect of the contact matrices on the infection transmission dynamics, we develop an age-structured SIR model to describe influenza transmission dynamics. The SIR model describes the spread of influenza in terms of the transition of individuals between different epidemiological compartments. Susceptible individuals (i.e., those at risk of acquiring the infection - S) can become infectious (i.e., capable to transmit the infection - I) after coming into contact with infectious individuals. Subsequently infectious individuals recover from the infection and become removed (R) after a certain amount of time (the infectious period). In an age structured implementation of the model, individuals are now identified also by their age, and the contact matrix is introduced to describe the number of contacts between susceptible individuals of age $i$ and all of their possible infectious contacts of age $j$ [2, 10] (see Sec. Methods for details). More specifically we considered a transmission model with identical disease parameters. The contact matrices are thus the only factor driving the difference in dynamics and attack rate (total number of infected individuals) of the simulated epidemic.

Compared to the case of homogeneous mixing, where all individuals are assumed to be in contact with each other in equal proportions, the inclusion of the contact matrices in the epidemic model consistently yields a lower overall attack rate for all locations (Fig. 6A). This difference is also reflected in the strong variability of the basic reproduction number $R_{0}$, representing the number of cases generated by a typical index case in a fully susceptible population, which depends on the spectral radius of the matrix $M$ as well as population structure (see SM). To provide further validation of the adequacy of the matrices in characterizing the specific dynamics of influenza transmission in the SM we report the simulations of the age-structured SIR model calibrated on real data from the H1N1 influenza pandemic in multiple locations. The model adequately reproduces the the age-specific seroprevalence profiles in Israel, Italy, Japan, the UK and the US [31, 32, 33, 34, 35].

To understand the underlying factors of the observed heterogeneities across geographical locations, we use a linear regression model to compare the attack rates and various socio-demographic features of each location (see SM). We identified two socio-demographic features that correlate strongly with the attack rate: the average age of the population (Fig. 6B) and the fraction of the population in the educational system including instructors (Fig. 6C). Indeed, if we examine the attack rates by age and setting (see SM), we observe that the greatest proportion of infections occur as a result of contact due to the school setting, and that attack rates in general are highest for school-aged individuals. Going further, inspection of the incidence profile by age (see SM) also clearly shows that individuals with high contact frequencies with others in the school setting are infected earlier in higher proportions. These results mirror well known influenza spreading trends/patterns observed in the real world [11, 20]. The observed results are robust (although with quantitative differences) to changes in transmissibility patterns and susceptibility to infection by age (see SM). Taken together, our results suggest that developing countries with younger populations, and thus more school-aged individuals, are likely to experience higher overall attack rates when compared to older, developed countries.

We can also investigate how the attack rate and $R_{0}$ for each location would differ if we only had knowledge of the contact patterns at the national level. In this scenario, we use the country-level influenza transmission contact matrices in each location (note that each location is still characterized by its own specific age structure) and compare the results with those obtained by using the location-specific contact matrices everything else kept identical for the disease transmission model (Fig. 7A-C). By using the country-level matrix we observe a much lower variability than by using location-specific mixing patterns. Moreover, location-specific attack rates and $R_{0}$ show a nonlinear relation with the results obtained using country-level contact patterns. Interestingly, we can observe clear geographical trends in the percent difference in attack rate using location-specific contact patterns in comparison to the corresponding country-level ones. For instance in much of the western area of China where most of the nation’s ethnic minorities live, using the average matrix would lead to underestimating the final impact of an epidemic, while we would overestimate it in the more traditionally urbanized/industrialized areas in the north-east of the country, like Beijing, and Shanghai (Fig. 7A).

To construct synthetic populations in different countries, we made use of a wide array of data sources (see SM). These data provide distributions of key socio-demographic characteristics such as the age structure, household size, age of the head of the household, age gaps between household members, household composition, employment rates, the educational system and enrollment rates, etc. Distributions such as these are typically available either as macro-level data from census databases and other governmental sources, or as micro-level data coming from surveys conducted on a sample of the population. Census databases routinely provide information at a broader scope such as the age structure of a population, or the fertility rates, however they often lack more detailed information related to the household composition and age relationships between household members. For this, we rely on micro-level surveys which collect data at the household and individual level and ask participants for information in regards to their health, household condition and composition, economic conditions, and more. The kind of data available also varies by country and even at the subnational level, thus necessitating the development of adaptive algorithms that can take in the available data and accommodate for variability in data organization to produce a faithful reconstruction of each population. With this in mind, the procedure implemented can be summarized as follows.

The first step in the reconstruction of a real-world population is the generation of households. In this process we use two types of multinomial sampling. The first is based on the probability distribution $\mathcal{M}(y)$ of an independent socio-demographic characteristic $y$. For instance, such characteristic $y$ can be the household size or composition, depending on the data available. The second type of multinomial sampling is based on the probability $\mathcal{M}(x|y_{1}=i_{1},y_{2}=i_{2},...,y_{n}=i_{n})$ of characteristic $x$ conditional on the value $i$ of a previously determined variable(s) $y$. In this case $x$ and $y$ are assumed (when supported by available data) to have bivariate or multivariate joint distributions. Typically, the larger the number of joint distributions incorporated, the more precise the reconstruction of the real-world population. The precision of such a reconstruction is however often limited by the scope of the data (such as the survey sample size for each characteristic $y$) and its availability. For example, of the multinomial joint distributions used here one is the distribution of the age of the head of the household by the size of the household and the household composition (whether a couple, a single parent with children, siblings, multigenerational families, etc.). The bivariate joint distributions incorporated is considerably long and includes (but is not limited to) distributions of the age of household members by the age of the head of the household, the age gap between couples living together by the age of one in the pair, the mother’s age at child birth by the age of the child, the number of household members by their relation to the age of the household head (such as a spouse, parent, child, grandchild, sibling, in-law, etc.) by the age of the household head and the household composition. These joint distributions were either found in the macro data or estimated from the micro survey data. Characteristics of the resulting synthetic households are compared to the distributions of the summary statistics available from the macro-level data using a goodness of fit test at the desired level of significance (generally 5%).

A similar procedure is used to assign those individuals to their respective schools and workplaces based on enrollment and employment records. These records detail the enrollment and employment rates by age, institutional sizes and their age structures, as well as the student-to-teacher ratios in the case of schools. A more detailed explanation of the construction of the synthetic population can be found in the SM together with the results of the comparison between the synthetic and actual population statistics.

We use the synthetic contact networks to infer average age-based contact patterns within each social setting. For each location, these age-based contact patterns are encoded in a contact matrix $F^{k}$, whose elements $F_{ij}^{k}$ describes the average frequency of contact between a given individual of age $i$ and individuals of age $j$ in setting $k$. We focus on 4 social settings: the households ($H$), schools ($S$), workplaces ($W$), and the general community ($C$). Specifically, here we adopt the frequency dependent (mass action) transmission model, with the implicit assumption that an increased population density has no effect on the per capita contact rate between individuals [36]. This choice of modeling mechanism was already proved to represent a good approximation for the description of the transmission patterns of several infectious diseases [2]. Moreover, it allows us to readily compare epidemiological parameters between social settings and locations with disparate population density, and thus makes for an appropriate framework when modeling the transmission dynamics of heterogeneous populations around the world. The calculation of the contact matrices can be described as follows.

First, we compute the relative abundance of contacts between individuals of age $i$ and individuals of age $j$ in each instance $s$ of setting $k$, $\Gamma_{ij}^{k(s)}$.

where $\phi_{i}^{k(s)}$ is the number of individuals of age $i$ in the instance $s$ (i.e., a specific household, school, or workplace) of setting $k$; $\delta_{ij}$ is the Kronecker delta function, which we use to omit the individual $i$ from their own set of contacts; $\nu^{k(s)}$ is the number of individuals (of all ages) in instance $s$ of setting $k$. Note that to compute $\Gamma_{ij}^{k(s)}$ we assume homogeneous mixing within each instance of the setting, and as a result the matrix $\Gamma_{ij}^{k(s)}$ has the expected symmetric property $\Gamma_{ij}^{k(s)}=\Gamma_{ji}^{k(s)}$.

Second, we compute the per capita probability of contact of an individual of age $i$ with an individual of age $j$ in setting $k$ as $F_{ij}^{k}$.

where $N_{i}$ is the total number of individuals of age $i$. Note that matrix $F^{k}$ (i.e., the matrix of elements $F_{ij}^{k}$) is not symmetric.

Third, we combine the setting-specific contact matrices by age $F^{k}$ to derive a matrix of the overall contacts by age $M$. We propose a weighted linear combination of the derived matrices in the four focus settings, calibrated to match the empirically estimated contact matrices from two contact diary survey studies in 7 locations throughout Western Europe and Russia [11, 15]. We perform a multiple linear regression to calibrate the weights of the synthetic setting contact matrices such that their linear combination matches the overall contact matrix for all 7 locations coming from the survey studies (see SM for details and for a comparison between the empirical and synthetic contact matrices).

The average number of contacts $<c>$ can be computed as

where $N=\sum_{i}N_{i}$ is the total number of individuals in the population. Therefore,

where $Z_{k}$ is the number of individuals having at least one contact in setting $k$.

To make side-by-side comparisons of the inferred contact matrices by age, we use the Canberra distance [20]. Specifically, each matrix is treated as a vector on which the Canberra distance is defined as

This yields a distance value of 0 for two locations with identical contact matrices, and increasingly larger distance values for two locations with increasingly different contact matrices.

For each location $l$ the transmission dynamics of influenza are modeled through an age-structured SIR model, where the mixing patterns are defined by the contact matrix previously introduced, $M_{ij}$.

The model is defined by the following set of equations:

where $S_{i}$ is the number of susceptible individuals of age $i$, $I_{i}$ is the number of infected individuals of age $i$, $R_{i}$ is the number of recovered or removed individuals of age $i$; $\gamma^{-1}$ is the infectious periods (which corresponds to the generation time in the simple SIR model [37, 38]), which is set to 2.6 days [39]; and $\lambda_{i}$ represents the force of infection to which an individual of age $i$ is exposed to other infected individuals and expressed as

where $\beta$ is the transmissibility of the infection, $N_{i}$ is the total number of individuals of age $i$, and $M_{ij}$ measures the average number of contacts for an individual of age $i$ with all of their contacts of age $j$.

The basic reproduction number $R_{0}$, representing number of cases generated by a typical index case in a fully susceptible population, can be defined for this model as

where $\rho(M)$ is the dominant eigenvalue of the matrix $M$ [40].

A database containing the inferred setting-specific matrices as well as the contact matrices for influenza transmission for all locations (and countries) is publicly available on the dedicated website: https://github.com/mobs-lab/mixing-patterns. Python routines to work with the contact matrices and examples of how to use them in age-structured compartmental models are also available.