Quantifying the dynamics of coupled networks of switches and oscillators.

Complex network dynamics have been analyzed with models of systems of coupled switches or systems of coupled oscillators. However, many complex systems are composed of components with diverse dynamics whose interactions drive the system's evolution. We, therefore, introduce a new modeling framework that describes the dynamics of networks composed of both oscillators and switches. Both oscillator synchronization and switch stability are preserved in these heterogeneous, coupled networks. Furthermore, this model recapitulates the qualitative dynamics for the yeast cell cycle consistent with the hypothesized dynamics resulting from decomposition of the regulatory network into dynamic motifs. Introducing feedback into the cell-cycle network induces qualitative dynamics analogous to limitless replicative potential that is a hallmark of cancer. As a result, the proposed model of switch and oscillator coupling provides the ability to incorporate mechanisms that underlie the synchronized stimulus response ubiquitous in biochemical systems.

Introduction

The dynamics in systems ranging from intercellular gene regulation to organogenesis are driven by complex interactions (represented as edges) in subcomponents (represented as nodes) in networks. If the structure of these networks is known, network-wide models of coupled systems have been applied to predict their qualitative dynamics. For example, models of coupled switches based upon Glass networks [1] have been applied to model systems such as neuronal synapses [2] and gene regulatory networks [3]. Similarly, models of coupled oscillators along networks based upon the Kuramoto model [4] have been used to model synchronization of oscillators in diverse systems reviewed in [5]. In biochemical systems, in vivo oscillator synchronization has been observed in synthetic oscillatory fluorescent bacteria [6], [7], yeast gene regulatory networks [8], [9], and human cell fate decisions [10]. Such spontaneous synchronization has also been attributed to the development of the mammalian cardiac pacemaker cells (reviewed in [11]) and cortical systems (reviewed in [12]) including notably the circadian pacemaker (e.g., [13]). More recently, these network models have been found to be insufficient to model more complex dynamics in neuronal information transfer [12], [14]–[17] and cardiac arrhythmias [18]–[21]. These limitations extend to physical systems, such as the coupled lasers studied in [22]. Therefore, numerous studies have modified these network models to account for evolving networks [15], [23]–[28], dynamic frequencies [15], [29], [30], or phase delays [16], [31]–[33]. However, these mathematical modifications typically do not encode the mechanism underlying the limitations in the Kuramoto and Glass network models.

We hypothesize that the observed limitations in the standard Kuramoto and Glass models arise from their exclusion of coupling components with qualitatively different dynamics. Several studies have inferred that biochemical systems contain “network motifs” with both oscillatory and switch-like dynamics [34], [35]. The dynamics of these motifs are inferred from the topology of subgraphs in the networks of these systems. Their structures are statistically overrepresented in biochemical networks [36], [37] such as intracellular regulatory networks [38], implicating evolutionary preservation (and thus utility) of these network motifs [39]. The dynamics of these motifs have been used to model yeast cell cycle regulation [40] and have been further confirmed in synthetic, designed biochemical circuits (reviewed in [41]). Because these heterogeneous network motifs are all identified as components within a single biochemical network, their interactions must drive the global dynamics of the network [42]. Previously, [43] have shown that coupling small sets of heterogeneous network motifs ensures the robustness of motif dynamics and [42] have shown that coupling networks changes their dynamics in isolation. However, the network-level dynamics that result from coupling oscillatory and switch-like components have not been studied comprehensively.

In this paper, we develop a theoretical framework to quantify the network-wide dynamics resulting from coupling switches and oscillators. This model is based upon introducing cross-coupling between the Kuramoto and Glass models, due to their wide success in modeling the dynamics in networks of oscillators and networks of switches, respectively. Simulations with the proposed model across state-space in an all-to-all network yields four operational states: (1) switches remain “on” and oscillators synchronize, (2) switches are “off” and oscillators freeze, (3) switches fluctuate in sync with oscillators, and (4) switches fluctuate transitionally until oscillators freeze. Further application of our model to the network motifs identified in yeast in [44] recapitulates the qualitative dynamics of the system observed in that study. However, a simple rewiring of this cell-cycle network that introduces feedback causes a cancer-like sustained re-activation of the cell cycle machinery without regard for external signal growth signals. These dynamics suggest that modeling cross-motif coupling may predict critical processes in the dynamics of biochemical networks with minimal parameterization.

The Kuramoto model of coupled oscillators

Quantitative studies of coupled oscillators often apply the Kuramoto model of oscillators coupled in an all-to-all network. In this model, the change in time of the phase of the oscillator, , is governed bywhere is the natural frequency of the oscillator and is the coupling strength of the oscillators [4]. Typically, the values are drawn from a normal distribution centered at 0 with variance .

In the Kuramoto model, the phases of the oscillators will synchronize if is above a threshold coupling strength . Such synchronization is quantified with the mean field of the oscillators asHere is the average phase of the oscillators and the coherence represents the spread of the oscillators from that average phase. Based upon eq. (2), if each and if the values of are distributed uniformly between [45].

Glass networks of coupled switches

Coupled sets of switches, which adopt one of a set of binary states, are modeled with Glass networks [1]. These models describe the evolution of the switch () as follows where represents the change in time of the value of each , which are unobservable continuous variables that control the time of switching between observable, discrete states in . In this model, describes the change in state of the switch due to the coupling with the other switches in the network [1]. In specified network structures and functions , such Glass networks can exhibit complex dynamics, including periodic and aperiodic orbits (e.g., [46]).

One type of Glass network, called a Hopfield network [2], has dynamics applicable to the smooth-decay of signal in biochemical switches [2]. The Hopfield model letswhere takes values between and representing the relative strength of the connection between switches and , is the magnitude of coupling strengths, and the threshold for switch activation. Similar to the Kuramoto model, sets of the switches will synchronize for above a threshold in appropriate network topologies.

Results

Network model of coupled oscillators and switches

By combining the established models for switches and oscillators, we model the dynamics of the heterogeneous system of coupled switches and oscillators in systems including biochemical networks with the following set of equations: Here, eq. (6) is analogous to the Glass network in eq. (3) and is defined according to eq. (4).

In this study, we explore a case of the switch-oscillator model in eqs. (6) and (7) which contains an all-to-all network that couples the Kuramoto model, eq. (1), and Hopfield network, eqs. (3)–(5), as follows where and are cross-component coupling strengths. In eq. (10), is the time-varying frequency of the oscillator resulting from switch coupling, with initial values , andIn this system, zero values of the cross-coupling parameters and cause the model to reduce to the standard uncoupled Kuramoto and Hopfield models. Similar decoupling of the models occurs if the switch and oscillator systems are at vastly different timescales, determined by the and parameters, respectively. The transformation in eq. (11) facilitates comparable switch-like dynamics in the oscillators when they interact with switches in eq. (8). Nonzero switch-oscillator () interactions will cause an oscillator in the “up” part () of its cycle to feed energy into the switch in question, nudging it towards the “on” state if off or delaying its decay if already on. Similarly, an “on” switch with a nonzero oscillator-switch interaction () will feed energy into the oscillators causing them to cycle at their natural frequency if coupled to that switch.

By thus incorporating coupling between switches and oscillators within the framework established by the standard Kuramoto and Hopfield models, the dynamics of our model in eqs. (8)–(10) can be analyzed within the framework of these well established models. Similar to analysis of the Kuramoto model and Glass network, we summarize the dynamics of our system using order parameters. For oscillators, we utilize the order parameter defined in eq. (2). We introduce a new order parameterthat tracks how closely each individual oscillator's frequency corresponds to the natural frequency . Analogously, we measure the fraction of switches that are in the “on” position at a given time using a switch-switch order parameter defined byBoth of these functions will have a maximum of 1 when all switches are on, and minimum 0 if all switches are off.

Simulation results in all-to-all networks

We first explore the qualitative dynamics of the heterogeneous system through numerical simulations in all-to-all networks. We limited these simulations to all-to-all networks, because of the ability of this network topology to describe the qualitative dynamics from the Kuramoto model. These simulations explore the majority of parameter space defined by , , , , and . Specifically, we select to ensure that switches are able to turn off without appropriate stimulation from the oscillators. We consider the effects of switches on oscillators for values of both above and below the Kuramoto threshold . Figure 1 plots the time-dependent order parameters observed in the four qualitative states observed in simulations of the coupled model eqs. (8)–(10) that are reflective of the qualitative dynamics observed in simulations with these parameter values. Movies S1– S4 further summarize the results of these simulations. We note that these four states were the only qualitative states observed for our coupled model in all-to-all networks simulated according to the description in the Methods section. Because , , and all control the relative timing of switches and oscillators, their values were selected in these simulations to optimize visualization in the supplemental videos. When exploring the effect of timing on the system dynamics, we hold and fixed while varying . Figure 2 shows the probability of observing the states in Figures 1(a)–1(c) in 100 simulations of all-to-all systems containing 100 switches and oscillators as a function of and . Because of their common control of system timing, we would obtain comparable distributions when varying either or instead of .

Figure 1

Summary of the qualitative dynamics of the heterogeneous network model of eqs. (8)–(10).

In all figures, top-panel shows temporal evolution of the mean field statistics ( black, solid; green, dashed; and blue, dash-dotted) and the bottom-panel shows the evolution of the mean phase (red, solid). (a) Oscillators synchronize and switches stay “on” (, , , , and ), (b) oscillators freeze (as evidenced by unchanging ) and switches stay “off” (, , , , and ), (c) oscillators synchronize and switches oscillate (, , , , and ), and (d) transitory oscillations in oscillators and switches (, , , , and ).

Figure 2

Percentage of simulations in which the qualitative dynamics in occur.

In (a) oscillators synchronize and switches are “on”, in (b) oscillators freeze and switches are “off”, and in (c) switches vary with oscillators vs for (solid), (dotted) and (dashed). , , and .

The coupled system preserves synchronization in both oscillators and switches

Figure 1(a) shows a state of the model in which the switches are all in the “on” state and oscillators are synchronized ( near 1, near 1, and oscillating between periodically). While such synchronization is observed in the uncoupled Hopfield and Kuramoto models, the oscillator-switch cross coupling extends the region of parameter space over which this synchronization occurs. Specifically, modest values of can induce sustained switch activity for parameter values of in which switches would decay in the uncoupled system. Furthermore, this switch synchronization will occur for all values of in which synchronization occurs in the uncoupled Hopfield model (i.e., all larger than a threshold value ) because the oscillators only contribute positively to the derivative in eq. (8) in our model. On the other hand, no value of will cause oscillator phases to synchronize if is below the critical coupling parameter for the pure Kuramoto model (). However, there are parameter regimes in which this synchronization occurs stochastically, depending on the initial values selected for , , and (Figure 2(a)). In these cases, the average decrease in oscillator natural frequencies caused by decreasing or will increase the effective period of oscillators, thereby increasing the probability of switches being locked in the “on” state and oscillator synchronization in the heterogeneous system.

Coupling switches to unsynchronized oscillators can freeze network-wide dynamics

Figure 1(b) depicts a model state in which switches are all “off” ( near zero) and oscillators “freeze”: each for some constant values for all beyond the preliminary freezing time . While the decaying switches are observed in an uncoupled Hopfield model, the freezing oscillators cannot be simulated in the uncoupled Kuramoto model. Such oscillator freezing will occur whenever the oscillators decay to the “off” state by virtue of the coupling of the oscillators to switches through the in eq. (10). Specifically, this frozen state can occur whenever depending on the values of , , and . However, the probability of selecting these initial states is decreased when the heterogeneity of the oscillators increases through incomplete synchronization () or increased (Figure 2(b)). In these cases, a single oscillator in the “up” phase () can contribute positively to the switch states, forcing the system out of this frozen state. The probability of obtaining this frozen model state further depends on the relative timing of switch decay and oscillator freezing. Specifically, the probability of obtaining the frozen state decreases with the average oscillator frequency, determined predominantly by the parameter (Figure 2(b)).

Coupling switches to synchronized oscillators can induce synchronized oscillations in switches

An additional consequence of coupling switches and oscillators in a state in which switches vacillate between all “on” and all “off” along with the synchronized oscillator frequency (Figure 1(c)). This oscillatory synchronization occurs when the pure Hopfield model would turn switches “off” (), the pure Kuramoto model would induce oscillator synchronization (), and the timing between the oscillators and switches are balanced such that the average period of the coupled oscillators is slightly less than the average decay time of the system of switches. Figure 2(c) shows that this balance in switch-oscillator timescales increases with decreasing and depends non-monotonically on . As we see in the plot of in Figure 1(c), the average oscillator natural frequencies will decrease towards the end of the “down” phase in response to switches turning off, and then increase to their full natural values in the “up” phase as switches turn back “on”. Therefore, if synchronized oscillator period is too slow (i.e., is too large), the system will tend to be locked in the “on” state (Figure 2(a)); if too fast (i.e., too small) the system will tend to be locked in the “off” state (Figure 2(b)).

Synchronization of network-wide oscillations may be transitory

Oscillatory behavior in the switches is also observed for unsynchronized oscillators () as depicted in Figure 1(d). In this case, the value of must be large enough to enable switches to freeze the oscillators' phases. However, because the oscillators are uncoupled, a small subset of oscillators in the “up” phase can drive the switches to turn on for large-enough values of . These switch oscillations are transitory, ending when at last the switch coupling dominates the system and induces all of the oscillators to freeze. For unsynchronized oscillators in the parameter range of Figure 1(d), the transitional oscillations in the switch state occurs regularly in 21 of 100 simulations. In 8 of these 21 simulations, the switch state turns “on” after decaying at least twice. More rarely, transitory changes in switch state may be induced by a similar mechanism in simulations for which and switches ultimately settle on the all “on” or all “off” states.

System size affects the distribution of qualitative dynamics

We also explored the dynamics of the coupled system for networks of sizes ranging from to nodes, described in the methods. For networks of all sizes, we observe that the dynamics of the system was limited to the four qualitative behaviors observed for networks of size depicted in Figure 1. However, the system size does have a notable effect on the frequency with which each of these behaviors occurs. The figures in Figures S1, S2, and S3 plot the observed frequencies for each of the network sizes as a function of the and values considered in Figure 2.

When , the observed frequencies of the system states depend most strongly on network size in simulations using the smallest value of is also small (Figure S1). In this case, the probability of observing the system with synchronized oscillatory dynamics in both switches and oscillators decays as the network grows. Both the state in which the switches are on and oscillators are synchronized and the state in which the switches are off and oscillators are frozen have with compensatory increases in probability (Figure 3). The relative probability of obtaining the frozen state increases, with notable decay in the probability of obtaining the state in which switches are “on” and oscillators are synchronized in large networks.

Figure 3

Dependence on network size for qualitative states for

and . Percentage of simulations in which the qualitative dynamics have switches off and oscillators frozen (blue, solid), switches on and oscillators synchronized (green, dashed), and oscillatory switches and synchronized oscillators (red, dotted).

On the other hand, when , the system size has the greatest influence on the resulting dynamics for large values of (Figure S3). In this case, the system changes from containing mostly switches in the on state and synchronized oscillators to switches that are entirely in the “frozen” state for large network sizes (Figure 4). We hypothesize that the system is forced into the frozen state in larger networks because of increased oscillator synchronization in large networks. Therefore, small networks would have a higher probability of having few oscillators that are unsynchronized and in the “up” phase (), causing the switches to turn “on” () due to the structure of eq. (8) as was discussed previously. Furthermore, the rare oscillations observed in both switches and oscillators when occur only when the network is small. Intermediate values of show similar changes to those described for when and to those described for when (Figure S2).

Figure 4

Dependence on network size for qualitative states for

and . Percentage of simulations with qualitative dynamics plotted as described in Figure 3.

The heterogeneous network models qualitative dynamics of the yeast cell cycle derived from network motifs

Previous work by [47] make the cell cycle processes controlling mitotic division of fission yeast Schizosaccharomyces pombe cells provides an optimal system in which to apply our model. The biochemical reactions responsible for driving the cell cycle are well understood and the resulting dynamics in each of the stages of the cell cycle have been characterized extensively in [40], [44], [47]. The cell cycle machinery in mitosis is divided into four, sequential stages: phase 1 is a gap or rest phase (G1); phase 2 is a DNA synthesis stage (S); phase 3 is an additional gap stage (G2); and phase 4 is the mitotic division stage (M). Previously, [47] observed that the dynamics of the yeast cell cycle can be divided into three sequentially interacting modules, triggered by a signal based upon cell size: (1) G1/S transitions with a toggle-switch, (2) S/G2 transitions with a toggle-switch, and (3) G2/M transitions with an oscillator. Although the specific timing differs from [47], we observe similar qualitative dynamics to that observed in [47] when applying our heterogeneous model to evolve the state of these cell cycle stages (Figure 5) as described in the Methods section. We note that the response in this system is consistent with the transitory oscillations observed in Figure 1(d) in the case of all-to-all coupling. We also modeled this cell-cycle system in a rewired-network, in which the G2/M transitions feedback into G1/S (Figure 6). In this case, we observe sustained reactivation of the cell cycle regardless of the external signal. These dynamics are analogous to the synchronized dynamics in Figure 1(c) and consistent with cell growth arising from re-wiring biochemical reactions in cancer cells [48].

Figure 5

Simulated dynamics of the yeast cell cycle.

Evolution of the states of the cell cycle modules (G1/S top, green dashed; S/G2 top, red dotted; G2/M bottom, black) in response to an external stimulus to initiate the cell cycle (top, blue solid).

Figure 6

Simulated dynamics of an aberrant cell cycle network.

As for Figure 5 with a network topology linking the G2/M module to the G1/S module.

Discussion

Our model of coupled switches and oscillators in all-to-all networks demonstrates that networks with components having heterogeneous dynamics can exhibit synchronization similar to that observed in homogeneous systems. As is the case in homogeneous models (e.g., [49]–[52]), we expect analogous synchronization to hold in small-world, biochemical network topologies (e.g., [53]). However, these alternative topologies would likely change the probability of observing each of the qualitative model behaviors similar to the observed dependence of probabilities in network size. In this alternative network topologies, the qualitative states of the network model may have greater variability in small network sizes in accordance with the findings of [54]. Finally, in these topologies the heterogeneous model could yield additional, complex qualitative dynamic states, resulting from the complex dynamics that they cause in models of coupled switches alone [46].

While uncoupled network motifs may adopt switch-like or oscillatory dynamics, coupling between these components can induce switch-like behavior in oscillators and oscillatory behavior in switches. These qualitative changes in component dynamics occur stochastically, depending on the distribution of frequencies and switch states. They are more likely to occur in simulations with an imbalance in relative timescales, in which the dynamics of the faster network motif will dominate the system. Similarly, when and are both small, the coordinated oscillations in the switches and oscillators that occur in frequently small networks are largely eliminated in larger networks. We hypothesize that this larger network effectively increases the range of natural frequencies and phases, making the simulation less likely to have the constrained distribution required to obtain such synchronized oscillations. We can expect that biological systems have evolved components according to these distributions to ensure the robustness of the dynamics in the system. For example, multiple proteins can often serve similar functions in cell signaling pathways, which would increase the system size and decrease the probability of transient behaviors in our model. This robustness will be further ensured through the sheer size of most biochemical systems. For example, in humans yeast two-hybrid maps and metabolic network maps both contain on the order of thousands of interactions between thousands of species [53].

Furthermore, we have also observed that the heterogeneous network model will freeze the oscillator dynamics in the presence of inactive switches and then subsequently activate in synchrony in the presence of active switches. As a result, our model provides a natural mechanism for the coordination of complex machinery such as the initiation of cell-cycle dynamics. For example, when we apply our model to the yeast cell cycle motifs in [47], we recapitulate the qualitative dynamics of delayed initiation of stages of the cell cycle observed in simulations with differential equations of the regulatory dynamics in [47]. Additional tuning of the model parameters or incorporation of additional cell cycle checkpoints would facilitate a precise match of the timing of [47]. Because parameters are defined for modules and their interaction, our model requires far fewer rate parameters than any differential equation model of sets of biochemical reactions of the yeast cell cycle. Generally, the oscillator in the final G2/M step of the cell cycle is active only when the series of switches in the previous steps of the cell cycle are activated, consistent with the transient dynamics observed in our network model. However, rewiring the network to introduce feedback from the G2/M stage to the G1/S stage of the cell cycle will cause the modeled cell cycle machinery to engage continually without regard to the external growth signals, consistent with the malignant rewring in cancer cells [48]. Similar to the oscillatory behavior induced in switches in simulations in all-to-all networks, this small modification to the topology of cell cycle interactions altered the resulting dynamics of the network motifs for the G1/S and S/G2 motifs. We, therefore, hypothesize that motif dynamics predicted by the structure of subgraphs may not accurately describe their in vivo dynamics if considered in isolation, consistent with the hypothesis in [55] and findings of [42].

We observed that the switches in the cell cycle block activation of the yeast cell cycle when no external signal is present. Similarly, when part of the larger but sparse networks that compose biochemical systems [53], inactive switches would effectively destroy links between nodes on the network. As a result, the proposed heterogeneous model provides a potential mechanism for Kuramoto-based models with evolving network topologies such as [15], [23]–[28]. Similarly, we observed that the intermediate switches delay the oscillations in the final G2/M motif in the simulated yeast cell cycle. As a result, we hypothesize that coupling switches to oscillators through their frequencies in this model also provides a natural mechanism for extensions of the Kuramoto model with dynamic frequencies [15], [29], [30] or phase delays [16], [31]–[33].

The heterogeneous network model described in this paper facilitates characterization of the dynamics of complex, biochemical systems by abstracting the dynamics of their composite motifs such as the yeast cell cycle based upon [47]. We note that the proposed heterogeneous network model is deterministic once the initial values of all the switches and oscillator frequencies have been specified. However, many intracellular reactions (e.g., [56]) and neuronal systems (reviewed in [57], [58]) evolve stochastically. In these cases, the Hopfield networks used to model the switches could be replaced with probabilistic Boolean networks [3] and the oscillators evolved with stochastic solvers such as the stochastic simulation algorithm (reviewed in [59]), integrated with the methodology developed in [60]. Similar modifications could also extend the heterogeneous model to incorporate coupling with components of additional dynamics pertinent to biochemical systems, such as those of the network motifs enumerated in [34], [35], [44]. These studies would also ideally consider the dynamics of the heterogeneous network model in additional small-world and random network topologies, as well as the topologies defined by neuronal systems and gene regulatory networks.

Materials and Methods

Numerical simulations in the all-to-all network

In this study, we analyze the range of possible dynamics of the coupled, heterogeneous networks by applying this model to all-to-all networks. Analyses were performed for networks with equal number of switches and oscillators () of sizes 10, 50, 100, 200, and 500. All simulations are run one hundred times from random initial conditions for the state of switches () and oscillators (), drawn from a Gaussian distribution and a uniform distribution on , respectively. Similarly, oscillator natural frequencies are drawn randomly from a Gaussian distribution of mean zero and standard deviation parameter . Simulations of 100 seconds (in the arbitrary units of the model), with a time step of 0.01 seconds were found sufficient to reflect the range of possible model behaviors and verify consistency across initial conditions. The behavior of each simulation is summarized based on the time-dependent order parameters and , , and .

Numerical simulations of the yeast cell cycle

Based upon [47], we model the yeast cell cycle as an initiating external signal (namely the cell size), coupled to a toggle switch representing the transition between G1/S, a toggle switch representing the transition between S/G2, and an oscillator representing the transition from G2/M. While the external signal is incorporated into the model with coupling to the other switches in eq. (6), its state is not updated by the model. The duration of this external signal is set at 10 simulated minutes, based upon [47]. Similarly, the initial values of the hidden variable for the switches in the G1/S and S/G2 modules are set at −0.5, to 1, and to 2 to reproduce the approximate 10 minute duration of these switches in [47]. The natural frequency is for the G2/M module set to to likewise reflect the timescale reported in [47], while the remainder of the coupling parameters are left untuned, set to because we sought only to reproduce the qualitative dynamics of the [47] model. The rewiring in the system with enduring cell cycle activation is achieved by adding an edge from the module for G2/M to the switch in the G1/S module.

Movie of the dynamics in . Left panel displays the evolution of the total number of switches (out of 100) in the on state (left) and off state (right). Right panel displays the evolution of the phase of each of the one hundred oscillators.

(AVI)

Click here for additional data file.

Movie of the dynamics in . As for Movie S1.

(AVI)

Click here for additional data file.

Movie of the dynamics in . As for Movie S1.

(AVI)

Click here for additional data file.

Movie of the dynamics in . As for Movie S1.

(AVI)

Click here for additional data file.

Dependence of dynamics on network size for . Number of simulations (of 100) for which the switches are off and oscillators are frozen (left panel), the switches are on and the oscillators are synchronized (center panel), and both the oscillators and switches have synchronized oscillations (right).

(PDF)

Click here for additional data file.

Dependence of dynamics on network size for . As for Figure S1.

(PDF)

Click here for additional data file.

Dependence of dynamics on network size for . As for Figure S1.

(PDF)

Click here for additional data file.