Emergence of network effects and predictability in the judicial system

As courts strive to simultaneously remain self-consistent and adapt to new legal challenges, a complex network of of citations between decided cases is established. Using network science methods to analyze the underlying patterns of citations between cases can help us understand the large-scale mechanisms which shape the judicial system. Here, we use the case-to-case citation structure of the Court of Justice of the European Union to examine this question. Using a link-prediction model, we show that over time the complex network of citations evolves in a way which improves our ability to predict new citations. Investigating the factors which enable prediction over time, we find that the content of the case documents plays a decreasing role, whereas both the predictive power and significance of the citation network structure itself show a consistent increase over time. Finally, our analysis enables us to validate existing citations and recommend potential citations for future cases within the court.


Introduction
Network representation of complex interactions among elements is an overarching framework heavily used in many fields of science [1][2][3] .For social systems, the dynamics of interactions between individuals (whether electronic, online or face-to-face) can be represented as time-varying networks, often called temporal networks, in which nodes come and go and edges are activated or deactivated as time goes on 4,5 .Many essential features of human behaviour encoded in the representation of temporal networks have been revealed over the past decade, such as burstiness [6][7][8] , circadian/diurnal rhythms 9 , temporal communities 10 , higher-order interactions 11, 12 , etc.While the studies of temporal networks shed light on the time-varying nature of interactions between nodes, dynamics in social systems emerge not only at the local level 13 , but also at the global level.In a wide variety of social contexts, network size (i.e., the number of active nodes) and the number of edges observed at a given point in time are very often not constant, and accordingly the average degree increases or decreases 14,15 .In fact, the numbers of aggregate nodes and edges have been shown to have a scaling relationship known as the densification power law or densification scaling 14 .In temporal networks (i.e., a sequence of snapshot networks), any variation in the number of active nodes N and the number of edges M can be a priori attributed to changes in (i) the population in the system (e.g., the number of students present in a school, the number of attendees in a conference, etc); (ii) the probability of two nodes being connected; or (iii) both.With a constant probability of edge creation, N and M will increase if more nodes enter the system, since each node will have a higher chance of finding partners.Likewise, for a given population, if the probability of two nodes being connected increases, M will surely increase, and N will rise as well as isolated nodes, if they exist, will be more likely to get connected.
These two mechanisms are fundamental factors that bring about the dynamics of N and M, yet separating their contributions based on the dynamical behaviour of N and M is a challenging problem.In a wide variety of social and economic systems, network dynamics are likely to be driven by a mixture of these two mechanisms, and moreover their relative importance may occasionally change as the network evolves 15 .In theory, each of these two mechanisms leads to a distinctive type of densification scaling; The first one, generated by the evolution of population, is a scaling behaviour similar to the typical densification scaling in which the number of edges M scales with the number of active nodes N with a constant exponent α, i.e., M ∝ N α 14 .The second one is an accelerating growth of M, which is caused by the evolution of the probability of edge creation 15 .In fact, for the human contact networks we study, neither of these two types of scaling is observed in their original form.Rather, we observe a "mixed" scaling behaviour which appears to be a composite of the two types and therefore cannot be explained by a single scaling law.
Here, we develop a Bayesian statistical method to identify the source of dynamics generating network densification and sparsification based on the sequence of N and M. To take into account possible changes in the source of dynamics, we derive two specifications (i.e., "regimes") for the solution of a simple generative model, namely a dynamic hidden variable model, each of which capturing one of the two fundamental mechanisms.By fitting the two specifications simultaneously to the observed mixed scaling relationship using a unified estimation framework, known as the Markov regime-switching model 16,17 , we are able to estimate the probability that the dynamical source of densification or sparsification at a given point in time is attributed to a particular mechanism.At the same time, the Bayesian inference also allows us to trace the paths of the time-varying parameters directly related to the dynamical source, i.e., the population in the system and the activity level of nodes.An important advantage of the regime-switching model is that it allows the "true" model specification to occasionally switch, possibly depending on the social context.
In this work we analyse networks of face-to-face human interactions collected by the SocioPatterns collaboration 18 .We focus on four datasets: contact networks in two scientific conferences, a hospital and a workplace.Such networks can indeed be affected by the two fundamental mechanisms at the same time, because (i) individuals can always enter and exit the system, and (ii) presence of a time schedule could facilitate or inhibit face-to-face interactions (e.g., attendees of a conference are more likely to have interactions during coffee breaks than during keynote talks).In particular, using data on academic conferences has an important advantage, as it allows us to compare the dynamical regimes detected by the proposed method with the "ground-truth" conference time schedules.We find indeed that during keynote talks, parallel sessions and coffee breaks, the temporal densification and sparsification in the contact networks formed by conference attendees are mainly related to shifts in the chance of contacts being made between attendees present at the venue.On the other hand, shifts in the population are the main driving force of densification and sparsification during registration and poster sessions.This result is consistent with our intuition that the number of attendees in the middle of the program would be mostly constant, while it may be more likely to change during registration, which is held in the morning, and poster sessions in which not all of the attendees participate.For contact networks in a hospital and a workplace, this kind of comparison with a prespecified time schedule is not possible because there is no such rigorous time constraints to follow.Nevertheless, in all the systems we examined, the proposed method reveals that the main driving force of network densification and sparsification is occasionally switching, suggesting that the formation of social ties in physical space generally involves multiple dynamical sources.

Empirical evidence on mixed densification scaling
We focus our analysis on temporal contact networks taken from the following four datasets: • WS-16: Contacts between participants of the Computational Social Science Winter Symposium 2016 at GESIS in Cologne on November 30, 2016 19 .
• IC2S2-17: 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, face-to-face 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 non-stationary 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 real-world 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.In upper panels, dynamical relationship between N and M is shown.Each dot represents a snapshot network created over a 10-minute time window.Gray dashed and dotted lines respectively denote N/2 (i.e., the lower bound for M) and N(N − 1)/2 (i.e., the upper bound for M).Lower panels show the behaviour of M over time.

Two dynamical regimes in the dynamic hidden-variable 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: P i j,t = κ t a i a j , i, j = 1, . . ., N p,t , t = 1, . . ., T. ( where a i is the "fitness" that represents the intrinsic activity level of node i 26-28 , and T denotes the last time window in the data. There are two time-varying parameters in the model.The first one is κ t > 0, which modulates the overall activity rhythm of nodes.A variation in κ would reflect the time-schedule 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 time-varying parameter N p,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 p,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 section 4.1 in Methods for derivation): where we drop time subscript t for brevity.From these expressions, it is clear that the two parameters κ and N p 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 κ or N p .Depending on whether we endogenise κ or N p , we obtain different functional forms that connect N and M.

Regime 1: N p -driven dynamics
First, let us consider the case of time-varying N p .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 p is time-varying and κ is constant, in which case the dynamical relationship between N and M is given by: where the time-varying N p value is expressed as a function of M t and κ: For the purpose of parameter estimation, we introduce an error term as , where κ denotes the estimated value of κ, and ε 1 t is a residual term following a normal distribution with mean zero and standard deviation σ 1 .Estimated value of N p,t when the system is in Regime 1 leads to: 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 time-varying nature of the population, what we call "N p -driven" dynamics.For a given κ, the slope of densification scaling is close to constant, while different κ yield different slopes.(Fig. 2, lower left).

Regime 2: κ-driven dynamics
Next, let us consider the case of time-varying κ.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 κ is time-varying and N p is constant, in which case the dynamical relationship between N and M is given by where the time-varying value of κ is expressed as a function of M t and N p : κ(M t , N p ) ≡ 8M t N p (N p −1) (see, Eq. 3).
For estimating, we add an error term as N t = h 2 (M t ; N p ) + ε 2,t , where N p denotes the estimated value of N p , and ε 2,t is a residual term following a normal distribution with mean zero and standard deviation σ 2 .Estimated value of κ at time t when the system is in Regime 2 leads to: In Regime 2, network dynamics is fully driven by the individuals' time-varying activity levels, what we call "κ-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 fixed-population system, in which case the network tends to be denser as the number of inactive nodes vanishes.

A Markov regime switching model
In real-world 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 s ∈ {1, 2} at time t for a given parameter set θ θ θ = {N p , κ, σ 1 , σ 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 4.2 for the estimation method).
In the following, we use the smoothed probability Pr(S t = s|ψ T ; θ θ θ ) which is calculated conditional on all the information available at time T , denoted by ψ T (see, Methods 4.3 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 s ∈ {1, 2}, we can estimate the dynamical parameters N p,t and κ t as: where θ θ θ 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 p -driven and κ-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 WS-16, 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 time-varying 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 (Figs. 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 regime-switching pattern (Figs. 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 Pr(S t = 1|ψ T ; θ θ θ ) 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 p -driven scaling), while the snapshots in Regime 2 exhibit accelerating growth patterns (i.e., κ-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 (WS-16 and IC2S2-17), the estimated population size N p,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 κ 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 (∼ 120 for WS-16, ∼ 200 for IC2S2-17).The variation of network size is thus mainly driven by the schedule, which constraints the participants' networking activity.In the case of WS-16, the fluctuation of N p,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 mixed-regime 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 regime-switching 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 S5), similar to the dynamics observed in a Hospital.

Non-monotonic 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(N − 1)), is not always increasing with N (Fig. 6).In fact, when the dynamics is N p -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 κ play a dominant role, the network density may increase when the network size is sufficiently large (Fig. 6, pale-red cross).This is because when the number of active nodes N is close to its upper bound N p , at which the activity levels of remaining inactive nodes are fairly low, the overall activity κ needs to be large enough for those low-activity 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 where q 0 denotes the fraction of isolated nodes in the system (see Eq. 16 in Methods 4.1).If the system is in Regime 1, in which κ is constant, the density monotonically approaches κ/4 as N p → ∞ (i.e., q 0 → 0 and N → ∞).On the other hand, if the system is in Regime 2, in which N p is constant, there is no a priori upper bound, and the density exhibits a non-monotonic behaviour.
In Regime 2, a change in κ has two opposing effects on the network density.First, an increase in κ directly increases density through a rise in the probability of edges being created.Second, a shift in κ would also increase N, which reduces the density through the third term in Eq. 10.Since q 0 → 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.

Discussion
Densification and sparsification of networks can occur for two reasons, namely a variation in the population N p and a variation in the overall activity level κ.A key finding of this work is that the relative importance of each of these two dynamical factors occasionally change, depending on the social context under study.By fitting the model to the observed scaling relations, we can detect the main factor that is relevant at a given point in time.Shifts in N p and/or κ affect the activity of all individuals equally, so these parameters could be considered as effective "temperatures" of the system.While in this work we studied face-to-face networks of individuals, by its versatility the proposed method could also be used for a wide variety of dynamical systems.There are some remaining issues for future research.First, the baseline model, a dynamic hidden variable model, relies on a "homogeneous mixing" hypothesis, which implies that nodes are connected to each other at random, given their activity levels.If we look at the structural properties of networks, such as triadic closure and community structure, such a hypothesis -especially for social contexts-would be unrealistic.However, the fact that the proposed method works remarkably well indicates that, as long as we look at network dynamics at a sufficiently coarse scale, keeping local properties aside, this homogeneous mixing assumption is a good approximation.In fact, introducing a non-random structure would easily make it impossible to obtain analytical expressions that would be needed for identification.Second, we assumed that the distribution of intrinsic node activities is uniform for simplicity.Ideally, one would need to set this distribution based on empirical evidence.However, measuring the empirical intrinsic activity levels of individuals is extremely difficult because one needs activity levels of totally inactive individuals as well.Furthermore, this parameter might very well have its own temporal evolution.If available, such a rich information would allow for a refinement of the method.Third, while the current method works well for temporal networks whose dynamical regime is occasionally switching, for fixed-regime systems in which the whole dynamics could be explained by either a N p -driven or a κ-driven regime, the proposed regime-switching model is unnecessary.In such cases, one would fit the empirical scaling to each of the two models separately, and then find out which model is better fitted 15 .
In many cases, examining the source of network dynamics from the level of each individual would be prohibitively difficult because each individual has his/her own circumstance, and privacy issues often prohibit researchers from obtaining enough information to reveal particular individuals' behaviour.In contrast, global quantities, such as the total numbers of nodes and edges, are much more widely accessible, and therefore utilising these quantities will be inevitable when high-resolution data are difficult to collect.A contribution of this work is that the proposed model allows us to detect the role of the two fundamental dynamical factors just by using information on the global network dynamics.Any dynamical processes occurring on networks, regardless of whether they are micro-or macro-phenomena, would be largely affected by the underlying dynamics of networks.This is in particular the case for spreading processes such as epidemics.A better understanding of the dynamics of densification and sparsification could thus benefit public health policies, which are of central importance for modern social systems.

Analytical expression for N and M
In this section we derive Eqs. ( 2) and ( 3).The numbers of active nodes N and edges M can be expressed as functions of parameters κ and N p (we drop time subscript t for brevity): where k(κ, N p ) denotes the average degree over all the existing nodes including isolated ones, and q 0 (κ, N p ) denotes the fraction of isolated nodes or equivalently the probability that a randomly chosen node being isolated.Let ρ(a) be the density of node activities, and define u(a, a ) as the probability that there is an edge between two nodes having activity levels a and a .The average degree k(κ, N p ) is given by the number of possible partners times the average of u(a, a ) (see, section S1 in SI for a full derivation): It should be noted that Eq. ( 12) is equivalent to the average degree in the standard fitness model 27 if N p − 1 is replaced with N, which is only asymptotically true in our model.The fraction of isolated nodes in the system is given by (see, section S1 in SI): Substituting ρ(a) = 1 (i.e., uniform distribution on [0, 1]) and u(a, a ) = κaa into Eq.( 12) leads to: Similarly, q 0 is given by: By defining a variable x ≡ 1 − κa 2 , we have: Combining these results with Eq. ( 11), we have: It should be noted that if |1 − κ/2| < 1 and N p is sufficiently large, then q 0 (κ, N p ) 0 and thereby N N p and M ∝ N 2 , as is shown in the study of the static fitness model [26][27][28] .
where {D D D t } denotes the sequence of observations D D D t = (N t , M t ), and f is given by: The log-likelihood function leads to: Bayesian inference is conducted based on the relationship p(θ θ θ |{D D D t }) ∝ L({D D D t }|θ θ θ )p(θ θ θ ), where p(θ θ θ |{D D D t }) and p(θ θ θ ) are posterior and prior densities, respectively.For each parameter we collect 20,000 samples (four chains, 5,000 samples after 5,000 burn-in for each chain) generated from the posterior using Markov chain Monte Carlo (MCMC).We implement MCMC using Pystan ver.2.19.0 32 , which runs the No-U-Turn sampler (NUTS) 33 .The mean parameter values are summarised in Table 1.Now we describe how information is updated in each period.The probability of being in state s conditional on information at time t is written as: where we drop argument θ θ θ in f for brevity.Given the initial guess for Pr(S 0 = r|ψ 0 ), we can recursively update the probability of being in state s.

Smoothed probability
The probability Pr(S t = s|ψ t ; θ θ θ ) obtained in Eq. ( 22) is based on information available at time t for a given parameter set θ θ θ .We can also obtain the probability based on all information, represented by information set ψ T .Let ξ ξ ξ t|T ≡ [Pr(S t = 1|ψ T ; θ θ θ ), Pr(S t = 2|ψ T ; θ θ θ )] be the vector of probabilities conditional on information at T .ξ ξ ξ t|T can be calculated by conducting backward iteration from T 16 : where and (÷) denote element-by-element multiplication and element-by-element division, respectively, and P P P = (p ss ) is the transition matrix.Note that all the terms in the RHS of the first equality are already known from the previous estimation procedure.After calculating ξ ξ ξ T −1|T , we use it to calculate the RHS of the second line.We repeat this until we obtain ξ ξ ξ t|T .

Validation
We check the accuracy of the inference method based on synthetic network data generated by the regime-switching hidden variable model.For given parameters N p , κ, p 11 , and p 22 , and time-varying variables {N p,t } and {κ t }, we generate sequences of {N t } and {M t } in a way prescribed in the model.When the network at t is in Regime 1 (Regime 2), the true N p,t (κ t ) is given by N p,t = 0.95N p,t−1 (κ t = 0.95κ t−1 ), and N p,t = N p (κ t = κ) otherwise.We assume that the initial probability of being in Regime 1 is set at 0.5, and N p,0 = N p and κ 0 = κ.For each parameter, we collect 20,000 samples by MCMC (5,000 samples from four chains after 5,000 burn-in iterations).
The estimated parameters under different sets of ground-truth θ θ θ are summarised in Table S1 in SI.The estimated smoothed probabilities well match the true states of the generated networks (Fig. S1 in SI).We also group the generated networks based on the probability of being in Regime 1; For each time period, if more than 95% of the sampled values for Pr(S t = 1|ψ t ; θ θ θ ) are higher (lower) than 0.5, then we classify the corresponding snapshot as being in Regime 1 (Regime 2).If it is not classified as Regime 1 or 2, the network is considered to be in a "gray area".As shown in the middle and the right columns of Fig. S1, the classification of generated networks based on estimated parameters is consistent with the ground truth, while there are some networks that are in gray zones especially when the observed pairs of (N t , M t ) are overlapped between the two regimes.A comparison between the estimated and the true paths of N p,t and κ t is also shown in Fig. S2.

Figure 1 .
Figure 1.Dynamical behaviour of number of active nodes N and number of active edges M. In upper panels, dynamical relationship between N and M is shown.Each dot represents a snapshot network created over a 10-minute time window.Gray dashed and dotted lines respectively denote N/2 (i.e., the lower bound for M) and N(N − 1)/2 (i.e., the upper bound for M).Lower panels show the behaviour of M over time.

Figure 2 .
Figure 2. Schematic of the identification method.Empirical densification is fitted to the regime switching model in which the model switches from Regime 1 to Regime 2 (resp.from Regime 2 to Regime 1) with probability p 12 (resp.p 21 ).Then, the estimated parameters are used to infer the probability of the system being in Regime 1 at a given time t.For panels at lower left and lower middle, different colours denote different N p , and different symbols denote different κ (see Eqs. 2 and 3).If the scaling is N p -driven (resp.κ-driven), the time variation of N and M is fully caused by shifts in N p (resp.κ).

Figure 4 .
Figure 4. Estimation of N p,t and κ t for (a) WS-16 and (b) IC2S2-17.N p,t and κ t are shown in the upper and the lower panels, respectively, and the 95 % credible interval is indicated by shading.Upper panels show the number of active nodes (dashed blue line) at each time, thus the difference between the two lines represents the number of isolated nodes.Lower panels also show the number of edges at each time (dashed red line).Vertical dotted lines indicate the time windows of the scheduled sessions, with the labels in the middle.

Figure S2 .
Figure S2.Validation of estimated parameters.True parameter value is denoted by red circle, and the mean of 20,000 parameter sampling and their 95% credible intervals are denoted by black solid and gray shading, respectively.
3, top right).The dynamics in the early morning and evening are driven by a variation in N p , as well as around lunch time and coffee break, and changes in activity level are the main source Identification of dynamical regime.Upper panels show the smoothed probability of being in Regime 1 (i.e., N p -driven dynamics) at each time window.95 % credible interval is indicated by shading.Lower panels show N-M plots with classified regimes being denoted by different colours and symbols.We identify a snapshot network as being in Regime 1 (resp. TimeNFigure3.Regime 2) if the estimated probability of being in Regime 1 (resp.Regime 2) is greater than 0.5 in more than 95 % of MCMC sampling.Otherwise, a network is considered as being in an undetermined "gray area".

Table 1 .
Estimated parameters.For each parameter, mean and 95% credible interval obtained by MCMC are shown at the upper and lower rows, respectively.N max denotes max t {N t }.This section describes how we can infer the model parameters and the dynamical regime at a given time interval t.Let Pr(S t = s|ψ t−1 ; θ θ θ ) be the probability that a network is in state s (i.e., in Regime s) conditional on information available at the end of time interval t − 1, denoted by ψ t−1 , for a given set of parameters θ θ θ = {N p , κ, σ 1 , σ 2 , p 11 , p 22 }.The likelihood function is then given by:

Table S1 .
Estimated parameters based on synthetic networks.For each parameter, mean and 95% credible interval obtained by MCMC are shown at the upper and lower rows, respectively.Parameters without hats, annotated in the first row and the first column, denote true values used in the generation of synthetic networks.