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]. arXiv:1111.7243v1 [q-bio.QM] 30 Nov 2011 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 M oscillators coupled in an all-to-all network. In this model, the change in timeθ i of the phase of the i th oscillator, θ i , is governed byθ whereω i is the natural frequency of the i th oscillator and κ θ,θ ≥ 0 is the coupling strength of the oscillators [4]. Typically, theω i 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 as Here ψ is the average phase of the oscillators and the coherence r θ represents the spread of the oscillators from that average phase. Based upon eq. (2), r θ = 1 if each θ i = ψ and r θ = 0 if the values of θ i are distributed uniformly between [0, 2π) [45].

Glass networks of coupled switches
Coupled sets of N switches, which adopt one of a set of binary states, are modeled with Glass networks [1].
These models describe the evolution of the i th switch (x i ) as followṡ whereẋ i represents the change in time of the value of each x i , which are unobservable continuous variables that control the time of switching between observable, discrete states inx i . In this model, F i describes the change in state of the i th switch due to the coupling with the other N switches in the network [1].
In specified network structures and functions F i , 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 lets where w ij takes values between −1 and 1 representing the relative strength of the connection between switches i and j, κ x,x is the magnitude of coupling strengths, and τ i the threshold for switch activation.
Similar to the Kuramoto model, sets of the switches will synchronize for κ x,x above a thresholdκ x,x in appropriate network topologies.

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) andx i 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 where κ x,θ and κ θ,x are cross-component coupling strengths. In eq. (10), ω l is the time-varying frequency of the l th oscillator resulting from switch coupling, with initial values ω l (t = 0) =ω l , and In this system, zero values of the cross-coupling parameters κ x,θ and κ θ,x 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 τ i andω l 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 (κ x,θ ) interactions will cause an oscillator in the "up" part (θ l = 1) 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 (κ θ,x ) 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 parameter that 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 by Both 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 κ x,x , κ x,θ , κ θ,θ , κ θ,x , and σ ω . Specifically, we select κ x,x = 1 <κ x,x 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  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 τ , κ x,θ , 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 κ x,θ fixed while varying σ ω .   (Figure 2(a)). In these cases, the average decrease in oscillator natural frequencies caused by decreasing κ θ,x 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. Specifically, this frozen state can occur whenever κ x,x <κ x,x depending on the values ofx i , θ j , and ω j . However, the probability of selecting these initial states is decreased when the heterogeneity of the oscillators increases through incomplete synchronization (r θ (t) < 1) or increased σ ω (Figure 2(b)). In these cases, a single oscillator in the "up" phase (θ = 1) 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 κ θ,x (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" (κ x,x <κ x,x ), 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 switchoscillator timescales increases with decreasing κ θ,x and depends non-monotonically on σ ω . As we see in the plot of r ω (t) 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 κ θ,x 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 κ x,θ . 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 N = M = 10 to N = M = 500 nodes, described in the methods. For networks of all sizes, we observe that the dynamics On the other hand, when κ θ,x = 1, the system size has the greatest influence on the resulting dynamics for large values of σ ω (Supplemental Figure S7). 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 (θ = 1), causing the switches to turn "on" (x = 1) due to the structure of eq. (8) as was discussed previously. Furthermore, the rare oscillations observed in both switches and oscillators when κ θ,x = 1 occur only when the network is small. Intermediate values of κ θ,x = 0.1 show similar changes to those described for κ θ,x = 0.01 when σ ω = 1 and to those described for κ θ,x = 1 when σ ω = 10 (Supplemental Figure S6).
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. 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 2(c) and consistent with cell growth arising from re-wiring biochemical reactions in cancer cells [48].

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][50][51][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 κ θ,x 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 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][24][25][26][27][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][32][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 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 r θ (t) and ψ (t), r ω (t), and r x (t).

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 x for the switches in the G1/S and S/G2 modules are set at -0.5, τ to 1, and κ x,x to 2 to reproduce the approximate 10 minute duration of these switches in [47].