The basic reproduction number as a predictor for epidemic outbreaks in temporal networks.

The basic reproduction number R0--the number of individuals directly infected by an infectious person in an otherwise susceptible population--is arguably the most widely used estimator of how severe an epidemic outbreak can be. This severity can be more directly measured as the fraction of people infected once the outbreak is over, Ω. In traditional mathematical epidemiology and common formulations of static network epidemiology, there is a deterministic relationship between R0 and Ω. However, if one considers disease spreading on a temporal contact network--where one knows when contacts happen, not only between whom--then larger R0 does not necessarily imply larger Ω. In this paper, we numerically investigate the relationship between R0 and Ω for a set of empirical temporal networks of human contacts. Among 31 explanatory descriptors of temporal network structure, we identify those that make R0 an imperfect predictor of Ω. We find that descriptors related to both temporal and topological aspects affect the relationship between R0 and Ω, but in different ways.

Introduction

The interaction between medical and theoretical epidemiology of infectious diseases is probably not as strong as it should. Many results in the respective fields fail to migrate to the other. There are of course exceptions. Perhaps the most important are the ideas of epidemic thresholds and the parameter R 0—the basic reproduction number—as a key predictor of the epidemiological severity of a disease [1,2]. R 0 is defined as the expected number of other individuals that an infected individual will infect if he or she enters a population entirely composed of susceptible individuals. It is thus a combined property of the process of contagion and the contact patterns of the population. In classic mathematical models of infectious disease spreading, R 0 = 1 marks an epidemic threshold. If R 0 < 1, the expected fraction of infected people in an outbreak, denoted by Ω, will not depend on the total population size N. If R 0 > 1, the expected value of Ω is proportional to N. In other words, in the limit of large populations, a finite fraction of the population can be infected. The focus on R 0 in the literature has sometimes been so strong that researchers rather calculate R 0 than quantities directly related to the outbreak, such as prevalence, incidence, and time to the peak prevalence.

The use of R 0 is not entirely unproblematic. First, it is hard to estimate both in models [3-5] and from outbreak data [6-8]. Second, the result that R 0 = 1 defines an epidemic threshold rests on very coarse assumptions [3,9,10]. For example, one needs to assume that every pair of individuals has the same chance of interacting at any given time. In fact, interaction rates depend on pairs of individuals—people living in the same city are more likely to interact than those living in different cities. The derivation of R 0 has been extended to the case in which information about contact networks (describing who can spread the disease to whom) is available [11-15]. In this case, the derivation is usually restricted to the case of regular networks, where all individuals have the same degree (number of neighbors in the contact network) [14,15]. Sometimes people use definitions of R 0 that differs from the original [11-13,16] in a strict sense (but typically captures some similar property relevant for the modeling framework in question). The assumption that a pair of individuals interacts at the same rate over time does not hold true in reality either. For example, interaction is more likely to take place when most people are awake. This point is a reason for the increasing interest in temporal networks (showing who is in contact with whom, at what time) as a representation for the interactions underlying epidemic spreading, which focus on time dependence of networks [17-19]. Another reason is the increasing availability of data sets of temporal networks—typically lists of anonymized id numbers of two individuals and the times when these two individuals have been in contact (close enough for a disease to spread). The temporal network literature has focused on spreading processes (not only epidemic spreading) and how these are affected by the structure. Structure, in this case, refers to the way the network differs from a random temporal network (where the contact can happen with any pair of nodes with equal probability, at any time). Studies of epidemic models on temporal networks have found that e.g. a broad distribution of interevent times slows down the spreading [18,19].

There have been a few attempts to examine R 0 for temporal networks. Ref. [16,20], for example, derives R 0 for a specific model of temporal networks. Ref. [21] measures R 0 in empirical temporal networks, but does not relate it to prevalence, final outbreak size or other direct measures of outbreak severity.

One possible approach in this line of research is to find more accurate estimators than R 0 of disease severity. However, R 0 is routinely estimated for different infectious diseases by public health organizations worldwide. These estimates constitute an important resource for monitoring and comparing disease outbreaks. Rather than discarding this data by proposing another quantity, we will investigate what R 0 really tells us about disease spreading in empirical temporal networks of human contacts. Including the temporal information can make a big impact on the outbreak dynamics compared to modeling epidemics on a static network, let alone a fully mixed model [17-19].

We use the Susceptible–Infectious–Recovered (SIR) model with constant disease duration [22]. This model has two control parameters—the probability of disease transmission (upon a contact between an infectious and susceptible individual), denoted by λ, and the duration of the infectious stage, denoted by δ. We numerically simulate the SIR model on various temporal networks. First, we observe that in this case Ω is not uniquely determined from an R 0 value. A combination of λ and δ can give a larger R 0 but a smaller Ω than another combination does. Then, we investigate how the structure of the temporal contact network explains the relationship between R 0 and Ω. Instead of building a theory that bridges the microscopic structure of temporal network data and the emergent properties of the outbreak, we screen many potentially interesting descriptors of the temporal network structure by identifying those that are strongly correlated with the descriptors of the shape of scatter plots of Ω vs. R 0.

Results

Empirical data

We analyze empirical sequences of contacts between people. These data sets can be divided into physical proximity and electronic communication data. The former type could be interesting for studying information and disease spreading mediated by human contacts. The latter type is primarily of interest in the context of information spreading (bearing in mind that information spreading not necessarily follows the same dynamics as infectious diseases). In all data sets, nodes are human individuals. We list some basic statistics of the data sets in Table 1.

Table 1

The basic statistics of the data sets.

	Number of individuals	Number of contacts	Sampling time	Time resolution
Conference	113	20,818	2.5d	20s
Dating	28,972	529,890	512.0d	1s
E-mail 1	57,189	444,160	112.0d	1s
E-mail 2	3,188	115,684	81.6d	1s
Facebook	293,878	876,993	4.36y	1s
Forum	7,084	1,412,401	8.61y	1s
Gallery	159(8)	6,027(350)	7.32(11)h	20s
Hospital	293,878	64,625,283	9.77y	1d
Online community	35,624	472,496	8.27y	1s
Prostitution	16,730	50,632	6.00y	1d
Reality mining	64	26,260	8.63h	5s
School	237(1)	62,886(2,263)	8.61(3)h	20s

The numbers in parenthesis indicate the standard deviation in order of the last digit for the two composite data sets (Gallery and School).

One data set belonging to the physical proximity class comes from the Reality mining study [23], where contacts between university students were recorded when their smartphones were within Bluetooth range (10–15 m). We use the same subset of this data as in Ref. [24]. Another class of proximity data was collected from groups of people wearing radio-frequency identification sensors. One such dataset comes from the attendees of a conference [25] (Conference), another from a school (School) [26], another from a hospital (Hospital) [27] and yet another from visitors to a gallery (Gallery) [25]. School and Gallery are collected for two and 69 days, respectively. We analyze the days separately and average the results over the days. In these data sets, a contact between people closer than 1~1.5m was recorded every 20 seconds. Finally, we use a data set of sexual contacts between sex sellers and buyers collected from a Brazilian web forum (Prostitution) [28].

The class of electronic communication data includes two e-mail networks. These data sets are described in detail in Refs. [29] (E-mail 1) and [30] (E-mail 2). E-mails have a natural direction from the sender to the recipient. However, to analyze all the data sets in the same way, we treat them as undirected temporal networks. We furthermore study two Internet communities: a dating community (Dating) [31] and a film community (Online community) [32]. The contacts in these data sets represent messages from person to person like e-mails do. In Dating there are also “flirts” with which one user expresses interest in another (but does not send text, images or other information). A slightly different form of online pair-wise interaction is posting to public web pages. We study one data set of posts to the home page (“wall”) of Facebook [33] and a data set from the aforementioned film community where a contact represents a reply to a post at a public forum (Forum) [32]. One contact in these data sets is thus a publically accessible message from one user to another.

Final outbreak size as a function of R 0

In Fig. 1, we show scatter plots of Ω vs. R 0 for our data sets (see S1 Fig. for the results for the Gallery data). One scatter plot corresponds to one data set. More precisely, we measure R 0 directly from the simulations according to the definition—the average number of others infected by the infection source. Ω is the fraction of recovered individuals once the outbreak has subsided, i.e., when there no longer are infectious individuals. A point in a scatter plot represents an average over 106 runs for given parameter values (λ,δ). Each run starts with one infected node that is selected from all nodes with the equal probability. We assume the source of the infection is infected at the time of the first contact. In total, we sample 20×20 points in the (λ,δ) parameter space, where each parameter varies from 0.001 to 1 with exponentially increasing intervals. δ is defined as a fraction of the total sampling time.

Fig 1

The average outbreak size plotted against the basic reproduction number for 12 data sets (indicated in the Fig.) of human interaction.

Each point of the scatter plots corresponds to one pair (λ,δ), where λ is the infection probability and δ is the duration of infection. In the upper left corner there is a legend for the color-coding of these points. In the other panels, a data point is an average over 104 runs of the SIR model as described in the Methods section. The vertical lines mark R 0 = 1—the epidemic threshold for the canonical, fully mixed SIR model.

For all the data sets, there is a significant deviation from a deterministic relationship between R 0 and Ω. Here, a deterministic relationship is operationally defined as the situation in which the Ω value is uniquely determined by the value of R 0 (as it would be in most fully mixed and network models we are aware of). Interestingly, the way these scatter plots deviate from a deterministic relationship depends on data sets. For example, for the Hospital data red points are typically on top of the green ones—i.e. points with higher λ and lower δ give larger outbreaks than points with similar R 0 but lower λ and higher δ. For the Facebook data the situation is reversed.

Characterizing the shape of the Ω vs. R 0 point cloud

To explore the causes of the imperfectness of R 0 as a predictor of Ω, we define six so-called shape descriptors, which measure the shape of the point clouds shown in Fig. 1. The shape descriptors are listed in Table 2, their definitions are illustrated in Fig. 2, and their values for each data set are shown in S2 Fig.

Table 2

Shape descriptors for the point clouds shown in Fig. 1.

Symbol	Definition	Explained in
τR0Ω	Kendall’s τ of R 0 vs. Ω (number of concordant pairs of parameter values—number of discordant pairs) / total number of pairs	Fig. 2A
ταΩ	λδ-balance, the Kendall’s τ of α vs. Ω averaged over ten equal sized bins of R 0. α is the angle to the λ = δ line of a point in λ,δ-space	Fig. 2D–F
ΔR0	Largest difference in R 0 among discordant pairs of points in R 0-Ω space, where one point corresponds to one combination of λ and δ.	Fig. 2B
ΔΩ	Largest difference in Ω among discordant pairs of points in R 0-Ω space	Fig. 2C
ρR0	Midpoint of the R 0 values of the pair defining ΔR0	Fig. 2B
ρΩ	Midpoint of the R 0 values of the pair defining ΔΩ	Fig. 2C

Fig 2

Explanation of shape descriptors to characterize the point clouds shown in Fig. 1.

All examples come from the Conference data set. Panel A describes Kendall’s τ—a correlation coefficient based on the counting of discordant pairs (pairs of points connected by a line of negative slope). Panels B and C show the maximal separation of discordant pairs. In B, the measures focus on the pair with the largest separation in the R 0 direction. Δ denotes the maximum separation; ρ is the mean R 0 value for the maximally discordant pair. Panel C shows the similar quantities, ΔΩ and ρΩ, defined along the Ω direction. Panels D, E, and F illustrate the measurement of λδ-balance via ταΩ. This descriptor captures the tendency of some data sets to have high-λ, low-δ points above high-δ, low-λ points, while for other data sets, the situation is reversed. Panel D illustrates how the R 0 axis is segmented into bins. Panel E shows how we assign a (λ,δ)-plane angle, α, to all points in the bin. Panel F shows how we measure the correlation between α and Ω, which is very weak in this particular case.

The first shape descriptors is the Kendall’s τ (Fig. 2A), which captures how good R 0 is as a predictor of Ω. We chose Kendall’s τ because the Ω vs. R 0 curve is highly non-linear such that the Pearson’s correlation coefficient would underestimate how good a predictor R 0 is. Among non-linear correlation measures, Kendall’s τ, is the most principled and easiest to understand. It counts the number of point pairs that are connected by a line with a positive slope (concordant pairs) and a negative slope (discordant pairs). Kendall’s τ is then the number of concordant pairs minus the number of discordant pairs divided by the total number of pairs. In the context of measuring the R 0-Ω correlation, we denote Kendall’s τ by τ.

Next four shape descriptors focus on the region in the (R 0,Ω) space where the spread of the points is the largest (Fig. 2B, C). We look for the discordant (λ,δ) pair with the largest difference between the its R 0 values. This difference defines Δ. Similarly, the largest difference in Ω among discordant pairs defines ΔΩ. We also measure the average R 0 value, ρ, of the two R 0 values derived from the discordant pair maximally separated in R 0. Similarly ρΩ is the average R 0 value of the discordant pair maximally separated in Ω. The shape descriptors ρ and ρΩ thus show the locations on the R 0 axis of the maximally separated discordant pairs. They may be related to the location of the epidemic threshold, where Ω takes off from zero in an infinite population.

As mentioned above, for some data sets, given a value of R 0, higher δ implies higher Ω (Hospital), whereas the relationship is reversed for other data sets. To quantify this observation, we define the sixth shape descriptor ταΩ that we call λδ-balance for short. To define ταΩ, we start by dividing the range of R 0 into ten equidistant bins between the smallest and largest observed values (Fig. 2D). Within a bin, the points have fairly similar R 0 values, but their λ and δ values can be diverse. To measure the effect of the balance between λ and δ on Ω, we calculate the angle α that a (λ,δ) pair relative to the origin makes to the diagonal in the (λ,δ)-plane, i.e., the λ = δ line (Fig. 2E). Then, we measure the correlation between α and Ω by Kendall’s τ (Fig. 2F). Finally, we average the values for the different bins. To avoid confusion, we denote the calculated Kendall’s τ by ταΩ.

Temporal and static network descriptors

To characterize the structure of the contact structures modeled as temporal networks, we use 31 different quantities, which we call network descriptors. They are listed in Table 3. We have chosen quantities that are relatively simple and intuitive.

Table 3

Descriptors of temporal network structure.

Type	Symbol	Definition
TE	f NC	Fraction of nodes present (i.e. having had at least one contact) when half of the contacts happened.
TE	f NT	Fraction of nodes present at half the sampling time.
TE	f LC	Fraction of links present when half of the contacts happened. This is illustrated in Fig. 3A and B.
TE	f LT	Fraction of links present at half the sampling time.
TE	F NC	Fraction of nodes present at both the first and last 5% of the contacts.
TE	F NT	Fraction of nodes present at both the first and last 5% of the sampling time. This is illustrated in Fig. 3C and D.
TE	F LC	Fraction of links present at both the first and last 5% of the contacts.
TE	F LT	Fraction of links present at both the first and last 5% of the sampling time.
LA	μLt	Mean of interevent times over links.
LA	σLt	Standard deviation of interevent times over links.
LA	cLt	Coefficient of variation of interevent times over links. In the terminology of Ref. [31], this is the burstiness of link activity.
LA	γLt	Skewness of interevent times over links.
LA	μLτ	Mean of the number of other contacts between two consecutive contacts of a link.
LA	σLτ	Standard deviation of the distribution of the number of other contacts in the data between two consecutive contacts of a link.
LA	cLτ	Coefficient variation of the distribution of the number of other contacts in the data between two consecutive contacts of a link.
LA	γLτ	Skewness of the distribution of the number of other contacts in the data between two consecutive contacts of a link.
NA	μNt	Like μLt but for nodes.
NA	σNt	Like σLt but for nodes.
NA	c Nt	Like c Lt but for nodes, i.e., the burstiness of node activity.
NA	γNt	Like γLt but for nodes.
NA	μNτ	Like μLτ but for nodes.
NA	σNτ	Like σLτ but for nodes.
NA	c Nτ	Like c Lτ but for nodes.
NA	γNτ	Like γLτ but for nodes.
DD	μk	Average degree of the network of accumulated contacts.
DD	σk	Standard deviation of the degree distribution of the network of accumulated contacts.
DD	c k	Coefficient of variation of the degree distribution of the network of accumulated contacts.
DD	γk	Skewness of the degree distribution of the network of accumulated contacts.
NS	N	Number of nodes.
NS	C	Clustering coefficient of the network of accumulated contacts.
NS	r	Degree assortativity of the network of accumulated contacts.

The types are: time evolution (TE), link activity (LA), node activity (NA), degree distribution (DD) and network structure (NS).

Time evolution

We calculate eight network descriptors that characterize the long-term behavior of the contact dynamics—basically, how the contacts process differs from a stationary process. The background is that some of these data sets (e.g. Prostitution, Dating, Forum and Online community) are growing throughout the sampling period. A fast-spreading outbreak would thus, effectively, spread in a larger population (defined as the set of individuals possible to be infected) in the end than in the beginning. The Gallery data is also special in that the individuals in the beginning of the sampling are not present in the end. Ref. [34] argues, in more general terms, that when the first and last contacts of a link (pairs of nodes that are in contact at least once) happen is important for the behavior of outbreaks.

The first such set of quantities focuses on the time when nodes and links appear for the first time. For example, Ref. [34] points at the growth of the Prostitution data set as a factor behind the observation [35] that the order of events speeds up disease spreading in this data. We use f to symbolize this class of network descriptors. We measure the fraction of links present at half the sampling time relative to the final number of links. Because several studies in temporal networks address the role of the order of events [35,36], rather than the time itself, we also measure the corresponding quantities if time is replaced by the contact index (the index of the contact number—1 for the first contact, 2 for the second, etc.). These have subscript ‘C’ as opposed to ‘T’ for time. Furthermore, the descriptors concerning nodes and links have the subscripts ‘N’ and ‘L’, respectively.

Another class of network descriptors, denoted F, focuses on persistent nodes or links. F is the fraction of nodes (subscript N) or links (subscript L) present in the first and last 5% of time (T) or contact index (C). Fig. 3 illustrates f and F.

Fig 3

Illustration of two descriptors of temporal network structure, f and F .

The measure illustrated in A and B, f , uses the order of the contact to separate the contacts; the measure in C and D, F , uses the real time. Panels A and C are time-line representations of a temporal network data set. Each horizontal line represents an individual. A contact between two individuals is indicated by a vertical arc. In A and B, we focus on the first contact between a pair of nodes. We measure the fraction of the number of node pairs that have been in direct contact when a fraction ν of the total number of contacts has been observed. This fraction is plotted against ν in B. The value at ν = 1/2 defines f . In the timeline (A) we highlight the first half contacts, which contribute to the calculation of f , in color and the first contact between each node pair by black contours. In panels C and D, we illustrate the calculation of F , which looks at nodes (rather than links) present in both the first and last time interval of width ϕ (measured as a fraction of the sampling time), shown in color in the timeline (C). The fraction of such nodes as a function of ϕ is graphed in D. F is defined as the value at ϕ = 0.05.

These network descriptors calculated across the different data sets span a relatively wide range. For example, f , the fraction of links present at half the sampling time, takes values from 0.17 (Facebook) to 0.98 (School).

Node and link activity

The node activity descriptors relate to the bursty nature of human activity as characterized by intense periods of activity separated by long periods of quiescence [37]. To characterize burstiness, one usually starts from interevent times, i.e., the times between consecutive contacts for a node or link. For simplicity, we ignore correlations between consecutive interevent times and focus on the probability distribution of interevent times. The distribution is often right-skewed—a structure that has been shown to slow down epidemic spreading [38-41]. To characterize the distribution, we measure four descriptors, i.e., the mean μ, standard deviation σ, coefficient of variation c (i.e. the standard deviation divided by the mean) [37], and the sample skewness given by where μ₂ and μ₃ are the second and third moments of the distribution, respectively.

Some studies have pointed out that the duration of presence of a node or link in the data can be more important for spreading dynamics than interevent times [34,42]. For this reason, we also study the distribution of node and link durations and use the same four descriptors. In sum, we use 16 network descriptors in this category—μ, σ, c and γ for interevent times and duration of activity, for both nodes and links.

Degree distribution

In the following, we define static network descriptors, i.e., those for aggregate contact networks. Among them, the degree distribution is arguably the most important for disease spreading. A right-skewed degree distribution, which is observed in many empirical networks, is known to facilitate disease spreading [43]. For simplicity, we use the network of accumulated contacts (even though one may be able to find network representations of temporal network data that better captures the important structures for disease spreading [44]). To summarize the shape of the degree distribution, we use the same four descriptors as for the interevent time and duration distributions—μ, σ, c and γ.

Other static network descriptors

We also measure other static network descriptors. First, we count the number of nodes, N. Because the number of links is equal to the half of the mean degree times N, we do not include it in the analysis.

We also measure the degree assortativity r (essentially, the Pearson correlation coefficient of the degrees at either side of a link). This network descriptor measures the tendency for assortative mixing by degree, i.e., whether high-degree nodes tend to connect to high-degree nodes and low-degree nodes to low-degree nodes. It has been shown that assortativity affects disease spreading (exactly how depends on the specific epidemic model and other structures of the contacts) [45-48].

Finally, we measure the clustering coefficient—the number of triangles in the network divided by the number of connected triples (not necessarily a full triangle) normalized to the interval [0,1]. Similar to assortativity, the relative number of triangles (clustering) is also a contact-structural factor influencing disease dynamics [46-51]. As an example, if we compare SI disease spreading on a clustered network with a random network with the same number of nodes and links, the early stage of the spreading would be faster in the less clustered network [49,50]. Intuitively, if a disease spreads from one individual to two neighbors, and the three individuals are connected as a triangle, then the third link of the triangle is useless for the spreading process. If the third link were connected elsewhere, the disease would spread faster.

Structural determinants of the Ω vs. R 0 point cloud

Ultimately, one would like to explain how the relations between R 0, Ω, λ and δ emerge from the contact structure. In this work, as mentioned, we take a different approach and look at the Pearson correlation coefficient between the shape descriptors (Table 2) and network descriptors (Table 3). In this way, we search for network descriptors that contribute to the deviation from a deterministic relationship between Ω and R 0. A temporal network data set defines a data point that is fed to the calculation of the correlation coefficient; there are 12 data points available for regression analysis. We decided to use the Pearson correlation coefficient and not multivariate regression methods because there are 31 dependent variables, i.e., network descriptors (and 6 independent variables, i.e., shape descriptors), whereas we have only 12 data points.

In Fig. 4, we plot the results from our correlation analysis. In each panel, we plot the coefficient of determination R 2 (square of the Pearson correlation coefficient) between a shape descriptor and each network descriptor. The network descriptors are grouped in accordance with the subsections of the previous section. Scatter plots of all pairs of network descriptor and shape descriptors are shown in S2 Fig.

Fig 4

The coefficient of determination R 2 between the shape descriptors of the R 0 vs. Ω point cloud and network descriptors.

The error bars are standard errors estimated by the jackknife resampling method. *: p < 0.05, **: p < 0.01, ***: p < 0.001.

The predictability of R 0 with respect to Ω, as measured by τ (Fig. 4A), is to some extent (p < 0.05) explained by the coefficients of variation of the interevent time for the node and link interevent time distribution, c and c . This correlation is positive (see S2 Fig.), so broader interevent time distributions (burstier contact patterns) imply worse predictability. Furthermore, the α dependence of Ω is most strongly correlated with the burstiness of the nodes c . In this case the correlation is negative (S2 Fig.). This means that if we compare two points with the same R 0 value, where the first parameter set has a comparatively large transmission probability and short disease duration than the second, then the first parameter set tends to trigger a larger outbreak size than the second. These quantities are strongly affected by burstiness. The remaining four shape descriptors concern the location (in the (R 0,Ω) space) of the biggest deviation from a deterministic relationship and the size of the deviation.

Fig. 4C shows the correlation coefficient with the location along the R 0 axis of the mid-point of the discordant pair with the largest separation in R 0, i.e., ρ. Also in this case, network descriptors derived from the interevent-time distributions are relatively strongly correlated with ρ. The mean μ and standard deviation σ as well as the skewness γ show strong correlations. Furthermore, the fraction of links present in both the first and last 5% of the contacts (F ) shows an R² = 0.4 correlation with ρ (p = 0.06). Furthermore, even though they do not reach the p < 0.05 significance criterion, other link-related quantities of the time evolution (μ, c , γ, μ, σ, c and γ) show R² values over 0.3. Fig. 4D indicates that the largest width of a discordant pair, Δ, is strongly correlated with a number of temporal network descriptors. First, Δ is correlated with both those relating to the node and link activity when the real time, not the contact index, is used (μ, c , γ, μ, σ and γ). Second, Δ is correlated with the time evolution, especially with the F quantities—measuring the fraction of links and nodes present both in the beginning and end of the sampling period (f , f , F , F , F , F ); p < 0.01). Fig. 4E shows the correlation with the R 0-location with the discordant pair with the largest separation in Ω, ρΩ. Just like ρ (Fig. 4C), much of the variance in ρΩ is explained by the time-related descriptors in real time (f , f , F , F and σ). More interestingly, the largest Ω-separation of discordant pairs, ΔΩ (Fig. 4F) is strongly and positively correlated with some static network descriptors, i.e., the coefficient of variation and the skewness of the degree distribution (c and γ).

Discussion

In this work, we have shown that temporal network structure of human contacts can change the interpretation of the basic reproduction number R 0. We have found pairs of SIR parameter values (λ1,δ1) and (λ2,δ2) such that R 0 (λ1,δ1) < R 0 (λ2,δ2) and Ω(λ1,δ1) > Ω(λ2,δ2). In other words, the expected number of secondary infections of the outbreak’s source is smaller for (λ1,δ1) than (λ2,δ2), but the expected final fraction of individuals that had the infection is larger for (λ1,δ1) than (λ2,δ2). It is hard to give a succinct explanation for this phenomenon, and we do not attempt that in the present paper. It relates to many aspects of the contact patterns—static network structures, dynamic network structures, and the fact that empirical data is finite-sized, non-equilibrium and inhomogeneous [18,19,52]. On the other hand, it is easy to imagine scenarios where this happens. Assume, for simplicity, that λ1≪ λ2, δ1 ≫ δ2 and the nodes split in two halves—one half active throughout the sampling time, the other half entering after some time. Then, in the (λ2,δ2) scenario, the larger λ (i.e., λ2) could cause a burnout outbreak that ends before the second group of nodes enters the system. Therefore, R 0 would be high, whereas Ω does not exceed 1/2. In the (λ1,δ1) scenario, R 0 would be smaller. However, the duration of infection would be long enough for the second half of the nodes to be infected, so Ω could be larger than 1/2. Therefore, a larger value of R 0 does not necessarily mean that the disease spreads more easily. At the same time, the correlation between R 0 and Ω is often strong, especially if one accepts a non-linear relationship. For most practical purposes, it probably suffices to assume that R 0 is a good predictor of Ω.

Looking closer at the deviation of the Ω vs. R 0 scatter plots from a deterministic relationship and structural correlates of the amount of the deviation, we notice that a combination of seemingly unrelated descriptors of temporal network structure often shows a significant correlation. This result suggests that—although a better achievement may be obtained through identification of microscopic factors contributing to these phenomena—such factors could be interdependent and hard to fully disentangle. Probably a fruitful path would be to vary the structure in models of contact patterns and look at responses in the Ω vs. R 0 plots. However, already based on the current numerical results, we can draw some conclusions. One of them is that the temporal network factors often seem important. In particular, the quantities relating to the interevent-time distributions are significant predictors of e.g. the overall correlation between Ω and R 0. This is a bit surprising in the light of Refs. [35] and [41] that have found that the birth and death of links and nodes influence (some other quantities relating to) spreading phenomena (probably also the importance of the “loyalty” metrics in Ref. [52]). Only one aspect of the Ω vs. R 0 plots—ΔΩ (see Table 2 and Fig. 2C for definition)—is primarily explained by the static network properties, specifically the coefficient of variation and skewness of the degree distribution. This result is accompanied by the largest confidence level (p < 0.001) of the correlation. In contrast to ΔΩ, a similar shape descriptor Δ (see Table 2 and Fig. 2B for definition) is strongly correlated with several of the temporal network properties and not with the static ones. Especially the former observation is interesting—even though temporal structure is needed to see any spread in ΔΩ at all, it is the degree distribution that is the most strongly correlated with the actual value of ΔΩ.

Needless to say, this work opens more questions than it answers. In particular, it calls for mechanistic modeling connecting R 0 and Ω. Another direction would be to develop improved estimators of disease severity.

Methods

In this section, we will go through technicalities of the methods that are not fully explained in the Results section.

SIR simulations

In this work we use the constant duration SIR model (that defines a Monte Carlo simulation of the SIR model). We initialize all individuals to susceptible and pick one random individual i to be the source of the infection. We assume that i becomes infected at the same time as its first appearance in the data. In a contact between an infectious and susceptible, the susceptible will (instantaneously) become infectious with a probability λ. Infectious individuals stay infectious for δ time steps after which they become recovered. If many contacts happen during the same time step, we go through them in a random order.

A more common version of the SIR model is to let infectious individuals recover with a constant rate. Qualitatively, both versions give the same results [21]. We use the constant duration version because it is a bit more realistic [53,54] and makes the code a bit faster than the exponentially distributed durations.

Measuring the λδ-balance

A combination of a large λ and small δ can give the same R 0 value as a combination of a small λ and large δ. At the same time, Ω may depend on one of these parameters more strongly than on the other. The result is a vertical trend in the colors of the points as seen in Fig. 1 (most clearly for the Forum, Dating and Online community data). We measure this tendency—the λδ-balance—as illustrated in Fig. 2D, E, and F. First, we segment the R 0 axis into ten bins. The number of bins is determined based on a trade-off between minimizing the spread of the points along the R 0 axis, and maximizing the number of points per bin. After the division into bins, we capture the λδ-balance via the angle α between the line from the origin to the parameter value (λ,δ) and the λ = δ line. Finally, we calculate Kendall’s τ for the relationship between α and Ω and average the τ values over all bins.

Data availability

The Conference, Gallery, Hospital and School data sets are available from http://www.sociopatterns.org/datasets/, the Prostitution data set is available as the Supporting Information of Ref. 35, and the Facebook data is available from http://konect.uni-koblenz.de/. Other data is available from the authors of the papers where they were first analyzed (as cited above).

R 0 vs. Ω plots for the Gallery data.

This file contains plots corresponding to Fig. 1 for all 69 days of data for the Gallery data set.

(PDF)

Click here for additional data file.

Scatterplots between the R 0-Ω cloud shape descriptors and network structural descriptors.

The numbers identify the data sets as follows: Conference (1), Dating (2), E-mail 1 (3), E-mail 2 (4), Facebook (5), Forum (6), Gallery (7), Hospital (8), Online community (9), Prostitution (10), Reality mining (11), School (12).

(PDF)

Click here for additional data file.