Empirical evidence on mixed densification scaling
We focus our analysis on temporal contact networks taken from the following four datasets:

WS16: Contacts between participants of the Computational Social Science Winter Symposium 2016 at GESIS in Cologne on November 30, 2016^{19}.

IC2S217: Contacts between participants of the International Conference on Computational Social Science 2017 at GESIS in Cologne on July 12, 2017^{19}.

Hospital: Contacts among patients, nurses, doctors and staffs in a Hospital in Lyon on December 8, 2010^{20}.

Workplace: Contacts between workers in a office building in France on June 27, 2015^{21}.
These data consist of contacts between individuals collected every 20 seconds using RFID sensors^{18,22}. A “contact” is here defined as a physical, facetoface proximity event. The datasets thus give us temporal networks in which nodes are individuals and edges encode the contacts occurring between them. All datasets exhibit large and abrupt fluctuations of the number of edges that are typical in these nonstationary systems (see Fig. 1, lower panels). In these particular contexts of social interactions, these transitions between high an low activity periods are often related to specified schedules: from talk sessions to coffee breaks in the conferences, changes in shifts in the hospital, from desk work to meetings in the workplace.
In many social and economic dynamical networks, the numbers of aggregate edges and nodes have a superlinear scaling relationship called the “densification power law”^{14,23,24}, in which the average degree is increasing with the number of nodes, i.e., “densification”. For temporal networks, where there is a sequence of network snapshots, a similar type of scaling emerges from the dynamics of the population, in which nodes enter and leave the system, keeping the chance of two nodes being connected constant^{15,25}. However, another type of scaling emerges in realworld systems for which the population is fixed. In such systems densification is “explosive”, with the scaling exponent increasing with N^{15}. While these two classes of scaling could be differentiated and identified from data if we observe a specific type of scaling^{15}, in general there may exist a mixture of them that cannot be easily classified as one of the two classes. Indeed, in the four datasets we study, no clear scaling relationship appears (Fig. 1, upper panels). In the following we show that the mixed shape of empirical densification behavior reflects a mixing of both classes of scaling.
Two dynamical regimes in the dynamic hiddenvariable model
To explore the temporal dynamics of densification and sparsification, we consider a dynamic version of the hidden variable model. The probability that two nodes i and j are in contact within a given time window t is:
$$beginaligned mathscr P_ij,t = kappa _t a_i a_j, ;;; i,j=1,ldots , N_mathrmp,t, ;; t = 1,ldots , T. endaligned$$
(1)
where (a_i) is the “fitness” that represents the intrinsic activity level of node i^{26,27,28}, and T denotes the last time window in the data. There are two timevarying parameters in the model. The first one is (kappa _t>0), which modulates the overall activity rhythm of nodes. A variation in (kappa) would reflect the timeschedule of a conference or a school, working hours in an office or a hospital, or the circadian rhythm of individuals^{9,22,29,30}. The second timevarying parameter (N_mathrmp,t) denotes the potential number of active nodes at time t, i.e., the total of active and inactive nodes that are in the room or the building. It should be noted that although the number of active nodes (i.e., nodes having at least one edge) (N_t) is always observable from the data, the potential number of nodes (N_mathrmp,t) is not. We do not usually know how many people were actually in the room at a given time because people could enter and exit the room at any time without being interacting with any other individual. We can observe the number of active nodes that appear in the record of contacts, but in many cases there is no record of nodes without any interaction. We assume that activity (a_i) is uniformly distributed on [0, 1], because (i) we do not have any prior information about the full distribution of the activity levels of all nodes including isolated ones, and (ii) introducing a more general distribution prohibits us from obtaining an analytical solution, which makes it difficult to implement parameter estimation.
The average numbers of active nodes N and edges M are analytically given as (see “Analytical expression for N and M” section in Methods for derivation):
$$beginaligned N&= N_mathrmpleft[ 1 frac2kappa N_mathrmpleft( 1left( 1frackappa 2right) ^N_mathrmpright) right] , endaligned$$
(2)
$$beginaligned M&= frac18 kappa N_mathrmp(N_mathrmp1), endaligned$$
(3)
where we drop time subscript t for brevity. From these expressions, it is clear that the two parameters (kappa) and (N_mathrmp) play different roles in the determination of N and M, but it is not clear how N and M correlate. To see the direct relationship between N and M, we eliminate one of the two parameters in Eq. (2), using Eq. (3). By doing this, we can effectively endogenise either (kappa) or (N_mathrmp). Depending on whether we endogenise (kappa) or (N_mathrmp), we obtain different functional forms that connect N and M.
Regime 1: (N_mathrmp)driven dynamics
First, let us consider the case of timevarying (N_mathrmp). This is a situation in which the dynamics of N and M are fully driven by changes in the population. We call this system as being in “Regime 1” or “state 1”:
Definition 1.
A system is in Regime 1 if (N_mathrmp) is timevarying and (kappa) is constant, in which case the dynamical relationship between N and M is given by:
$$beginaligned N_t&= N_mathrmp(M_t,kappa ) left[ 1 frac2kappa N_mathrmp(M_t,kappa )left( 1left( 1frackappa 2right) ^N_mathrmp(M_t,kappa )right) right] nonumber \&equiv h^1(M_t;kappa ), endaligned$$
(4)
where the timevarying (N_mathrmp) value is expressed as a function of (M_t) and (kappa): (N_mathrmp(M_t,kappa ) equiv frac1+ sqrt1+32M_t/kappa 2) (see, Eq. 3).
For the purpose of parameter estimation, we introduce an error term as (N_t=h^1(M_t;widehatkappa ) + varepsilon _1,t,) where (widehatkappa ) denotes the estimated value of (kappa), and (varepsilon _t^1) is a residual term following a normal distribution with mean zero and standard deviation (sigma _1). Estimated value of (N_mathrmp,t) when the system is in Regime 1 leads to:
$$beginaligned widehatN_mathrmp,t_S_t=1 = frac{1+ sqrt1+32M_t/widehatkappa }2, endaligned$$
(5)
where (S_t=1) denotes the fact that the system is in Regime 1 at time t. In Regime 1, network dynamics is totally driven by the timevarying nature of the population, what we call “(N_mathrmp)driven” dynamics. For a given (kappa), the slope of densification scaling is close to constant, while different (kappa) yield different slopes (Fig. 2, lower left).
Regime 2: (kappa)driven dynamics
Next, let us consider the case of timevarying (kappa). This corresponds to a situation in which the dynamics of the system is fully driven by changes in the overall activity of individuals. We call this system as being in “Regime 2” or “state 2”:
Definition 2.
A system is in Regime 2 if (kappa) is timevarying and (N_mathrmp) is constant, in which case the dynamical relationship between N and M is given by
$$beginaligned N&= N_mathrmpleft[ 1 frac2kappa (M,N_mathrmp)N_mathrmpleft( 1left( 1frackappa (M,N_mathrmp)2right) ^N_mathrmpright) right] nonumber \&equiv h^2(M;N_mathrmp), endaligned$$
(6)
where the timevarying value of (kappa) is expressed as a function of (M_t) and (N_mathrmp): (kappa (M_t,N_mathrmp) equiv frac8M_tN_mathrmp(N_mathrmp1)) (see, Eq. 3).
For estimating, we add an error term as (N_t=h^2(M_t;widehatN_mathrmp) + varepsilon _2,t,) where (widehatN_mathrmp) denotes the estimated value of (N_mathrmp), and (varepsilon _2,t) is a residual term following a normal distribution with mean zero and standard deviation (sigma _2). Estimated value of (kappa) at time t when the system is in Regime 2 leads to:
$$beginaligned widehatkappa _t_S_t=2 = frac8M_twidehatN_mathrmp(widehatN_mathrmp1). endaligned$$
(7)
In Regime 2, network dynamics is fully driven by the individuals’ timevarying activity levels, what we call “(kappa)driven” dynamics, and the slope of densification scaling in fact increases with N (Fig. 2, lower middle). This kind of accelerating growth of M naturally happens when edges are created in a fixedpopulation system, in which case the network tends to be denser as the number of inactive nodes vanishes.
Analysis of switching dynamics behind temporal densification and sparsification
A Markov regime switching model
In realworld networks, the mechanism of densification and sparsification may occasionally change depending on the context, such as working schedule, coffee breaks, lunch time, etc. To incorporate such a possibility, we propose a unified framework based on the Markov regime switching model in which the hidden state of a system can switch from Regime 1 to Regime 2 (respectively from Regime 2 to Regime 1) with probability (p_12) (resp. (p_21))^{16,17}. An important advantage of the regime switching model is that it allows us to calculate the probability of a system being in Regime (sin 1,2\) at time t for a given parameter set (varvectheta =N_mathrmp,kappa ,sigma _1,sigma _2,p_11,p_22\). This probability of the system being in Regime s can then be interpreted as the relevancy of each mechanism in explaining the densification dynamics at a given time (Fig. 2). We employ a Bayesian approach for the estimation of the parameters, using the Markov chain Monte Carlo (MCMC) to obtain posterior distributions (see, Methods “Bayesian estimation” for the estimation method). It should be noted that (varvectheta ) does not contain (N_mathrmp,t_S_t=1) and (kappa _t_S_t=2) since (varvectheta ) only contains constant parameters to be directly estimated by the Bayesian method.
In the following, we use the smoothed probability (mathrmPr(S_t=spsi _T;varvectheta )) which is calculated conditional on all the information available at time T, denoted by (psi _T) (see, Methods “Smoothed probability” for full derivation) ^{31}. Validation analyses using synthetic networks show that the proposed method correctly detects the switching of regimes and estimates the model parameters quite accurately (Table S1, Figs. S1 and S2 in Supporting Information (SI)). Given the probability of being in Regime (sin 1,2\), we can estimate the dynamical parameters (N_mathrmp,t) and (kappa _t) as:
$$beginaligned widehatN_mathrmp,t&= mathrmPr(S_t=1psi _T;widehatvarvectheta )cdot widehatN_mathrmp,t_S_t=1 ; +; mathrmPr(S_t=2psi _T;widehatvarvectheta )cdot widehatN_mathrmp, endaligned$$
(8)
$$beginaligned widehatkappa _t&= mathrmPr(S_t=1psi _T;widehatvarvectheta )cdot widehatkappa ; + ; mathrmPr(S_t=2psi _T;widehatvarvectheta )cdot widehatkappa _t_S_t=2, endaligned$$
(9)
where (widehatvarvectheta ) denotes the set of estimated parameters, which is summarised in Table 1 in Methods.
Classification of network dynamics
The Bayesian estimation of the parameters suggests that the empirical systems’ dynamics are indeed occasionally switching between (N_mathrmp)driven and (kappa)driven (Fig. 3, upper panels). For the conference data, a common feature is that the probability of being in Regime 1 is almost 1 prior to the first session and after the last keynote session of the day, and mostly zero in between (see Fig. 4 for the correspondence between the dynamics and the schedule of the conferences). For WS16, we see further fluctuations between the two regimes, one linked to the lunch break, the other to the poster session which closed the day. This suggests that the dynamics during the oral sessions, keynote talks and breaks are mainly driven by changes in the activity level of participants, while in the “opened” time slots, such as registration, closing and poster session, their dynamics are explained by timevarying population. The same patterns linked to the schedule are found on the other days of the conferences (see S3a–c).
For the Workplace data we see a roughly similar pattern (Fig. 3, top right). The dynamics in the early morning and evening are driven by a variation in (N_mathrmp), as well as around lunch time and coffee break, and changes in activity level are the main source of dynamics in between. This is, of course, not necessarily a general property of contact networks in physical space. We also see that the regime remains almost constant in most of the day (Fig. S3e), or there might be days in which the regime constantly changes throughout the day (Fig. S3f). In the case of the Hospital data, there is no clear tendency for the regimeswitching pattern (Fig. 3, third column and S3d), which seems natural for such an open environment with visitors and medical workers coming and going, and no general, fixed schedule for working hours.
We next attempt to classify the snapshot networks into two groups based on their probability of being in a particular regime. We identify a snapshot network at t as being in Regime 1 (resp. Regime 2) if more than 95 % of samples for the value of (mathrmPr(S_t=1psi _T;varvectheta )) generated by MCMC are greater than 0.5 (resp. lower than 0.5), i.e., in more than 95 % of parameter sampling the dynamics at t is considered to be attributed to Regime 1 (resp. Regime 2). Otherwise, the system is considered as being in an undetermined “gray area”. As seen in the lower panels of Fig. 3, the location of snapshot networks in the N–M space is strongly related to which regimes they belong to. As expected, the snapshots in Regime 1 exhibit a scaling whose slope is almost constant (i.e., (N_mathrmp)driven scaling), while the snapshots in Regime 2 exhibit accelerating growth patterns (i.e., (kappa)driven scaling). Classifying each time window according to the underlying dynamical mechanism is essentially equivalent to identifying patterns in the N–M space.
Temporal dynamics of population and activity level
We also examine the evolution of the dynamical parameters for both regimes (Fig. 4). For the two conferences (WS16 and IC2S217), the estimated population size (widehatN_mathrmp,t) increases at the beginning of the day and decreases at the end, consistent with the dynamics of participants entering and exiting the venue. The estimated activity parameter (widehatkappa _t) is high during these periods, and the level is consistent with those seen in highly active windows during social breaks. During the main program, the population is virtually constant and the size is consistent with the number of attendants ((sim) 120 for WS16, (sim) 200 for IC2S217). The variation of network size is thus mainly driven by the schedule, which constrains the participants’ networking activity. In the case of WS16, the fluctuation of (widehatN_mathrmp,t) during the lunch break and the poster session are worth noting since the variation of observed network size N seems to be driven by both mechanisms; we see slight reductions in the estimated population while the overall activity is still high in these time windows. This demonstrates the ability of the proposed method to extract mixedregime periods in which both of the two mechanisms are at work (see Fig. 2, right, for schematic). Similar patterns are also found in the other days (see S4).
In the Hospital data, the regimeswitching dynamics is much less periodic, with lots of transitions and mixed periods (Fig. 5a). This is however not surprising, because there is no fixed schedule regulating either the activity or the number of people present. For the Workplace data, we also do not expect a priori to see a clear segmentation of regimes because of the absence of a rigid schedule as in a hospital. However, the dynamics uncovered by our method indicates that the situation is much simpler than that for Hospital, as there seem to be less variation in population size, aside from the “opening” and “closing” effects and a reduction in population around the lunch time (Fig. 5b). The day that exhibited many regime switches presents however many episodes of small variations in population size (see S4), similar to the dynamics observed in a Hospital.
Nonmonotonic behaviour of network density
Since both types of scaling emerging from two different dynamics exhibit superlinearity, the average degree is always increasing in N. However, the density of networks, defined by (2M/(N(N1))), is not always increasing with N (Fig. 6). In fact, when the dynamics is (N_mathrmp)driven, the network density mostly decreases as the network size N increases (Fig. 6, blue circle). So, a rise in N causes the density to be smaller when the engine of dynamics is changes in population. In contrast, when changes in (kappa) play a dominant role, the network density may increase when the network size is sufficiently large (Fig. 6, palered cross). This is because when the number of active nodes N is close to its upper bound (N_mathrmp), at which the activity levels of remaining inactive nodes are fairly low, the overall activity (kappa) needs to be large enough for those lowactivity nodes to get at least one edge. This would necessarily increase the total number of edges in the network to a large extent, which leads to a “true” densification of networks.
These properties are also confirmed by the analytical equation for the average network density^{15}
$$beginaligned frac2MN(N1) = frackappa 4 left( frac11q_0(kappa ,N_mathrmp)right) ^2 left( 1+fracq_0(kappa ,N_mathrmp)N1right) , endaligned$$
(10)
where (q_0) denotes the fraction of isolated nodes in the system (see Eq. 16 in Methods “Analytical expression for N and M”). If the system is in Regime 1, in which (kappa) is constant, the density monotonically approaches (kappa /4) as (N_mathrmprightarrow infty) (i.e., (q_0rightarrow 0) and (Nrightarrow infty)). On the other hand, if the system is in Regime 2, in which (N_mathrmp) is constant, there is no a priori upper bound, and the density exhibits a nonmonotonic behaviour. In Regime 2, a change in (kappa) has two opposing effects on the network density. First, an increase in (kappa) directly increases density through a rise in the probability of edges being created. Second, a shift in (kappa) would also increase N, which reduces the density through the third term in Eq. 10. Since (q_0rightarrow 0) as N becomes sufficiently large, the latter effect is vanishing, and therefore the density begins to rise with N for a sufficiently large N.