Quantifying the dynamics of coupled networks of switches and oscillators

Matthew R. Francis1, Elana J. Fertig2,∗,

1 Physics Department, Randolph-Macon College, Ashland, VA, USA
2 Department of Oncology and Division of Oncology Biostatistics and
Bioinformatics, Sidney Kimmel Comprehensive Cancer Center, School of Medicine,
Johns Hopkins University, Baltimore, MD, USA
$\ast$ E-mail: Corresponding<EMAIL_ADDRESS>
## Abstract

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 [14, 15, 16, 12,
17] and cardiac arrhythmias [18, 19, 20, 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, 24, 25, 26, 27, 28], dynamic frequencies [15, 29, 30], or
phase delays [16, 31, 32, 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
$M$ oscillators coupled in an all-to-all network. In this model, the change in
time $\dot{\theta}_{i}$ of the phase of the $i^{\mbox{th}}$ oscillator,
$\theta_{i}$, is governed by

$\dot{\theta_{i}}=\hat{\omega}_{i}+\frac{\kappa_{\theta,\theta}}{M}\sum_{j=1}^{M}\sin\left(\theta_{j}-\theta_{i}\right),$
(1)

where $\hat{\omega}_{i}$ is the natural frequency of the $i^{\mbox{th}}$
oscillator and $\kappa_{\theta,\theta}\geq 0$ is the coupling strength of the
oscillators [4]. Typically, the $\hat{\omega}_{i}$ values are drawn from a
normal distribution centered at 0 with variance $\sigma_{\omega}$.

In the Kuramoto model, the phases of the oscillators will synchronize if
$\kappa_{\theta,\theta}$ is above a threshold coupling strength
$\hat{\kappa}_{\theta,\theta}$. Such synchronization is quantified with the
mean field of the oscillators as

$r_{\theta}e^{i\psi}=\frac{1}{M}\sum_{j=1}^{M}e^{i\theta_{j}}.$ (2)

Here $\psi$ is the average phase of the oscillators and the coherence
$r_{\theta}$ represents the spread of the oscillators from that average phase.
Based upon eq. (2), $r_{\theta}=1$ if each $\theta_{i}=\psi$ and
$r_{\theta}=0$ if the values of $\theta_{i}$ are distributed uniformly between
$[0,2\pi)$ [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^{\mbox{th}}$ switch ($\tilde{x}_{i}$) as follows

$\dot{x}_{i}=-x_{i}+F_{i}\left(\tilde{x}_{1},\tilde{x}_{2},\ldots,\tilde{x}_{n}\right)\mbox{,
and}$ (3) $\tilde{x}_{i}=0\mbox{ if }x_{i}<0\mbox{; }1\mbox{ otherwise},$ (4)

where $\dot{x}_{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 in $\tilde{x}_{i}$. In this
model, $F_{i}$ describes the change in state of the $i^{\mbox{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

$F_{i}=\kappa_{x,x}\sum_{j=1}^{N}w_{ij}\tilde{x}_{j}-\tau_{i},$ (5)

where $w_{ij}$ takes values between $-1$ and $1$ representing the relative
strength of the connection between switches $i$ and $j$, $\kappa_{x,x}$ is the
magnitude of coupling strengths, and $\tau_{i}$ the threshold for switch
activation. Similar to the Kuramoto model, sets of the switches will
synchronize for $\kappa_{x,x}$ above a threshold $\hat{\kappa}_{x,x}$ 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:

$\displaystyle\dot{x}_{i}$ $\displaystyle=$ $\displaystyle-
x_{i}+G_{i}\left(\tilde{x}_{1},\tilde{x}_{2},\ldots,\tilde{x}_{N},\theta_{1},\theta_{2},\ldots,\theta_{M}\right)$
(6) $\displaystyle\dot{\theta}_{l}$ $\displaystyle=$
$\displaystyle\omega_{l}\left(\tilde{x}_{1},\tilde{x}_{2},\ldots,\tilde{x}_{N}\right)$
$\displaystyle+$ $\displaystyle
H_{l}\left(\tilde{x}_{1},\tilde{x}_{2},\ldots,\tilde{x}_{N},\theta_{1},\theta_{2},\ldots,\theta_{M}\right).$

Here, eq. (6) is analogous to the Glass network in eq. (3) and $\tilde{x}_{i}$
is defined according to eq. (4).

In this study, we explore a case of the switch-oscillator model in eqs. (6)
and (Network model of coupled oscillators and switches) which contains an all-
to-all network that couples the Kuramoto model, eq. (1), and Hopfield network,
eqs. (3) - (5), as follows

$\displaystyle\dot{x}_{i}$ $\displaystyle=$ $\displaystyle-
x_{i}+\frac{\kappa_{x,x}}{N}\sum_{j\neq.~{}i}^{N}\tilde{x}_{j}+\frac{\kappa_{x,\theta}}{M}\sum_{k=1}^{M}\tilde{\theta}_{k}-\tau_{i},$
(8) $\displaystyle\dot{\theta}_{l}$ $\displaystyle=$
$\displaystyle\omega_{l}+\frac{\kappa_{\theta,\theta}}{M}\sum_{k=1}^{M}\sin\left(\theta_{k}-\theta_{l}\right),$
(9) $\displaystyle\dot{\omega}_{l}$ $\displaystyle=$
$\displaystyle\frac{\kappa_{\theta,x}}{N}\sum_{j=1}^{N}(\tilde{x}_{j}\hat{\omega}_{l}-\omega_{l}),$
(10)

where $\kappa_{x,\theta}$ and $\kappa_{\theta,x}$ are cross-component coupling
strengths. In eq. (10), $\omega_{l}$ is the time-varying frequency of the
$l^{\mbox{th}}$ oscillator resulting from switch coupling, with initial values
$\omega_{l}(t=0)=\hat{\omega}_{l}$, and

$\tilde{\theta}_{l}=\left\\{\begin{array}[]{ll}1&\mbox{ if
}0\leq\theta_{l}<\pi\\\ 0&\mbox{ otherwise}\\\ \end{array}\right..$ (11)

In this system, zero values of the cross-coupling parameters
$\kappa_{x,\theta}$ and $\kappa_{\theta,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 $\tau_{i}$ and $\hat{\omega}_{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 ($\kappa_{x,\theta}$) interactions will cause an
oscillator in the “up” part ($\tilde{\theta}_{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 ($\kappa_{\theta,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

$r_{\omega}(t)=\frac{1}{M}\sum_{\mu}\frac{\omega_{\mu}(t)}{\hat{\omega}_{\mu}}$
(12)

that tracks how closely each individual oscillator’s frequency $\omega_{\mu}$
corresponds to the natural frequency $\hat{\omega}_{\mu}$. 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

$r_{x}(t)=\frac{1}{N}\sum_{j}\tilde{x}_{j}(t).$ (13)

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 $\kappa_{x,x}$,
$\kappa_{x,\theta}$, $\kappa_{\theta,\theta}$, $\kappa_{\theta,x}$, and
$\sigma_{\omega}$. Specifically, we select $\kappa_{x,x}=1<\hat{\kappa}_{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 $\kappa_{\theta,\theta}$ both above and below the Kuramoto threshold
$\hat{\kappa}_{\theta,\theta}$. 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. Supplemental
videos 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 $\tau$, $\kappa_{x,\theta}$, and $\sigma_{\omega}$
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
$\tau$ and $\kappa_{x,\theta}$ fixed while varying $\sigma_{\omega}$. Figure 2
shows the probability of observing the states in Figures 1-1 in 100
simulations of all-to-all systems containing 100 switches and oscillators as a
function of $\kappa_{\theta,x}$ and $\sigma_{\omega}$. Because of their common
control of system timing, we would obtain comparable distributions when
varying either $\tau$ or $\kappa_{x,\theta}$ instead of $\sigma_{\omega}$.

#### The coupled system preserves synchronization in both oscillators and
switches.

Figure 1 shows a state of the model in which the switches are all in the “on”
state and oscillators are synchronized ($r_{x}$ near 1, $r_{\theta}$ near 1,
and $\psi$ oscillating between $[0,2\pi)$ 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
$\kappa_{x,\theta}$ can induce sustained switch activity for parameter values
of $\kappa_{x,x}$ in which switches would decay in the uncoupled system.
Furthermore, this switch synchronization will occur for all values of
$\kappa_{x,x}$ in which synchronization occurs in the uncoupled Hopfield model
(i.e., all $\kappa_{x,x}$ larger than a threshold value $\hat{\kappa}_{x,x}$)
because the oscillators only contribute positively to the derivative in eq.
(8) in our model. On the other hand, no value of $\kappa_{\theta,x}$ will
cause oscillator phases to synchronize if $\kappa_{\theta,\theta}$ is below
the critical coupling parameter for the pure Kuramoto model
($\hat{\kappa}_{\theta,\theta}$). However, there are parameter regimes in
which this synchronization occurs stochastically, depending on the initial
values selected for $x_{i}$, $\theta_{j}$, and $\hat{\omega}_{j}$ (Figure 2).
In these cases, the average decrease in oscillator natural frequencies caused
by decreasing $\kappa_{\theta,x}$ or $\sigma_{\omega}$ 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 depicts a model state in which switches are all “off” ($r_{x}$ near
zero) and oscillators “freeze”: each
$\theta_{j}\left(t\right)=\psi\left(t\right)=\Psi$ for some constant values
$\Psi$ for all $t$ beyond the preliminary freezing time $t_{f}$. 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
$\omega_{j}$ in eq. (10). Specifically, this frozen state can occur whenever
$\kappa_{x,x}<\hat{\kappa}_{x,x}$ depending on the values of $\tilde{x}_{i}$,
$\theta_{j}$, and $\hat{\omega}_{j}$. However, the probability of selecting
these initial states is decreased when the heterogeneity of the oscillators
increases through incomplete synchronization ($r_{\theta}(t)<1$) or increased
$\sigma_{\omega}$ (Figure 2). In these cases, a single oscillator in the “up”
phase ($\tilde{\theta}=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 $\kappa_{\theta,x}$ (Figure 2).

#### 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). This oscillatory synchronization
occurs when the pure Hopfield model would turn switches “off”
($\kappa_{x,x}<\hat{\kappa}_{x,x}$), the pure Kuramoto model would induce
oscillator synchronization
($\kappa_{\theta,\theta}>\hat{\kappa}_{\theta,\theta}$), 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 shows that this balance in switch-oscillator
timescales increases with decreasing $\kappa_{\theta,x}$ and depends non-
monotonically on $\sigma_{\omega}$. As we see in the plot of
$r_{\omega}\left(t\right)$ in Figure 1, 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., $\sigma_{\omega}$ is too large), the system will
tend to be locked in the “on” state (Figure 2); if too fast (i.e.,
$\sigma_{\omega}$ too small) the system will tend to be locked in the “off”
state (Figure 2).

#### Synchronization of network-wide oscillations may be transitory.

Oscillatory behavior in the switches is also observed for unsynchronized
oscillators ($\kappa_{\theta,\theta}<\hat{\kappa}_{\theta,\theta}$) as
depicted in Figure 1. In this case, the value of $\kappa_{\theta,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
$\kappa_{x,\theta}$. 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, 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 $\kappa_{\theta,\theta}>\hat{\kappa}_{\theta,\theta}$ 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 of the system was limited
to the four qualitative behaviors observed for networks of size $N=M=100$
depicted in Figure 1. However, the system size does have a notable effect on
the frequency with which each of these behaviors occurs. Supplemental Figures
S5-S7 plot the observed frequencies for each of the network sizes as a
function of the $\kappa_{\theta,x}$ and $\sigma_{\omega}$ values considered in
Figure 2.

When $\kappa_{\theta,x}=0.01$, the observed frequencies of the system states
depend most strongly on network size in simulations using the smallest value
of $\sigma_{\omega}=1$ is also small (Supplemental Figure S5). 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.

On the other hand, when $\kappa_{\theta,x}=1$, the system size has the
greatest influence on the resulting dynamics for large values of
$\sigma_{\omega}$ (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
($\tilde{\theta}=1$), causing the switches to turn “on” ($\tilde{x}=1$) due to
the structure of eq. (8) as was discussed previously. Furthermore, the rare
oscillations observed in both switches and oscillators when
$\kappa_{\theta,x}=1$ occur only when the network is small. Intermediate
values of $\kappa_{\theta,x}=0.1$ show similar changes to those described for
$\kappa_{\theta,x}=0.01$ when $\sigma_{\omega}=1$ and to those described for
$\kappa_{\theta,x}=1$ when $\sigma_{\omega}=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. 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 [47,
44, 40]. 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 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 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 $\kappa_{\theta,x}$
and $\sigma_{\omega}$ 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, 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 [44, 34, 35]. 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
($N=M$) of sizes 10, 50, 100, 200, and 500. All simulations are run one
hundred times from random initial conditions for the state of switches
($x_{i},i=1,\ldots,N$) and oscillators ($\theta_{j},j=1,\ldots,M$), drawn from
a Gaussian distribution and a uniform distribution on $[0,2\pi)$,
respectively. Similarly, oscillator natural frequencies are drawn randomly
from a Gaussian distribution of mean zero and standard deviation parameter
$\sigma_{\omega}$. 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 $r_{\theta}\left(t\right)$ and
$\psi\left(t\right)$, $r_{\omega}\left(t\right)$, and $r_{x}\left(t\right)$.

### 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, $\tau$ to 1, and $\kappa_{x,x}$ 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
$\frac{2\pi}{60}\mbox{min}^{-1}$ to likewise reflect the timescale reported in
[47], while the remainder of the coupling parameters are left untuned, set to
$\kappa_{\theta,x}=\kappa_{x,\theta}=\kappa_{\theta,\theta}=2$ 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.

## Acknowledgments

This work was sponsored by NCI (CA141053). We thank LV Danilova, B Fertig, AV
Favorov, BR Hunt, MJ Stern, MF Ochs, E Ott, E Webster, and LM Weiner for
advice. Code available upon request.

## References

  * 1. Glass L, Kauffman S (1973) The logical analysis of continuous, non-linear biochemical control networks.  Journal of Theoretical Biology 39: 103-129.
  * 2. Hopfield JJ (1982) Neural networks and physical systems with emergent collective computational abilities.  Proc Natl Acad Sci U S A 79: 2554-8.
  * 3. Shmulevich I, Dougherty ER, Zhang W (2002) From boolean to probabilistic boolean networks as models of genetic regulatory networks.  Proceedings of the IEEE 90: 1778-1792.
  * 4. Kuramoto Y (1975) International Symposium on Mathematical Problems in Theoretical Physics, volume 39 of _Lecture Notes in Physics_.  Springer, 420 pp.
  * 5. Strogatz SH (2003) Sync: the emerging science of spontaneous order.  New York: Hyperion, 1st ed edition.
  * 6. Garcia-Ojalvo J, Elowitz MB, Strogatz SH (2004) Modeling a synthetic multicellular clock: repressilators coupled by quorum sensing.  Proc Natl Acad Sci U S A 101: 10955-60.
  * 7. Danino T, Mondragón-Palomino O, Tsimring L, Hasty J (2010) A synchronized quorum of genetic clocks.  Nature 463: 326-30.
  * 8. Tu BP, Kudlicki A, Rowicka M, McKnight SL (2005) Logic of the yeast metabolic cycle: temporal compartmentalization of cellular processes.  Science 310: 1152-8.
  * 9. Klevecz RR, Li CM, Marcus I, Frankel PH (2008) Collective behavior in gene regulation: the cell is an oscillator, the cell cycle a developmental process.  FEBS J 275: 2372-84.
  * 10. Huang S, Eichler G, Bar-Yam Y, Ingber D (2005) Cell fates as high-dimensional attractor states of a complex gene regulatory network.  Phys Rev Lett 94: 128701.
  * 11. Peskin C (1975) Mathematical aspects of heart physiology.  Courant institute lecture notes, Courant Institute of Mathematical Sciences.
  * 12. Breakspear M, Heitmann S, Daffertshofer A (2010) Generative models of cortical oscillations: neurobiological implications of the kuramoto model.  Front Hum Neurosci 4: 190.
  * 13. Strogatz SH (1987) Human sleep and circadian rhythms: a simple model based on two coupled oscillators.  J Math Biol 25: 327-47.
  * 14. Fox JJ, Jayaprakash C, Wang D, Campbell SR (2001) Synchronization in relaxation oscillator networks with conduction delays.  Neural Comput 13: 1003-21.
  * 15. Cumin D, Unsworth C (2007) Generalizing the Kuramoto model for the study of neuronal synchronization in the brain.  Physica D 226: 181-196.
  * 16. Wang Q, Duan Z, Perc M, Chen G (2008) Synchronization transitions on small-world neuronal networks: Effects of information transmission delay and rewiring probability.  EPL 83: 50008.
  * 17. Pérez T, Garcia GC, Eguíluz VM, Vicente R, Pipa G, et al. (2011) Effect of the topology and delayed interactions in neuronal networks synchronization.  PLoS One 6: e19900.
  * 18. Kunysz A, Shrier A, Glass L (1997) The challenge of complex dynamics for theoretical models of cardiac activity.  Journal of Theoretical Medicine 1: 79–90.
  * 19. Chang HS, Staras K, Gilbey MP (2000) Multiple oscillators provide metastability in rhythm generation.  J Neurosci 20: 5135-43.
  * 20. Sato D, Xie LH, Sovari AA, Tran DX, Morita N, et al. (2009) Synchronization of chaotic early afterdepolarizations in the genesis of cardiac arrhythmias.  Proc Natl Acad Sci U S A 106: 2983-8.
  * 21. Zahanich I, Sirenko SG, Maltseva LA, Tarasova YS, Spurgeon HA, et al. (2011) Rhythmic beating of stem cell-derived cardiac cells requires dynamic coupling of electrophysiology and ca cycling.  J Mol Cell Cardiol 50: 66-76.
  * 22. Tiana-Alsina J, Garcia-Lopez JH, Torrent M, Garcia-Ojalvo J (2011) Dual-lag synchronization between coupled chaotic lasers due to path-delay interference.  Chaos 21: 043102.
  * 23. Mondal A, Sinha S, Kurths J (2008) Rapidly switched random links enhance spatiotemporal regularity.  Phys Rev E Stat Nonlin Soft Matter Phys 78: 066209.
  * 24. Sorrentino F, Ott E (2008) Adaptive synchronization of dynamics on evolving complex networks.  Phys Rev Lett 100: 114101.
  * 25. Jiang M, Ma P (2009) Coherence resonance induced by rewiring in complex networks.  Chaos 19: 013115.
  * 26. Volman V, Perc M (2010) Fast random rewiring and strong connectivity impair fast random rewiring and strong connectivity impair subthreshold signal detection in excitable networks.  New Journal of Physics 12: 043013.
  * 27. Adhikari BM, Prasad A, Dhamala M (2011) Time-delay-induced phase-transition to synchrony in coupled bursting neurons.  Chaos 21: 023116.
  * 28. Gómez-Gardeñes J, Gómez S, Arenas A, Moreno Y (2011) Explosive synchronization transitions in scale-free networks.  Phys Rev Lett 106: 128701.
  * 29. Martens E, Barreto E, Strogatz S, Ott E, So P, et al. (2009) Exact results for the kuramoto model with a bimodal frequency distribution.  Phys Rev E 79: 026204.
  * 30. Taylor D, Ott E, Restrepo JG (2010) Spontaneous synchronization of coupled oscillator systems with frequency adaptation.  Phys Rev E Stat Nonlin Soft Matter Phys 81: 046214.
  * 31. Ermentrout B, Ko TW (2009) Delays and weakly coupled neuronal oscillators.  Philos Transact A Math Phys Eng Sci 367: 1097-115.
  * 32. Wang Q, Perc M, Duan Z, Chen G (2009) Synchronization transitions on scale-free neuronal networks due to finite information transmission delays.  Phys Rev E Stat Nonlin Soft Matter Phys 80: 026206.
  * 33. Wang Q, Chen G, Perc M (2011) Synchronous bursts on scale-free neuronal networks with attractive and repulsive coupling.  PLoS One 6: e15851.
  * 34. Alon U (2007) Network motifs: theory and experimental approaches.  Nat Rev Genet 8: 450–461.
  * 35. Novak B, Tyson JJ (2008) Design principles of biochemical oscillators.  Nat Rev Mol Cell Biol 9: 981–991.
  * 36. Milo R, Shen-Orr S, Itzkovitz S, Kashtan N, Chklovskii D, et al. (2002) Network motifs: simple building blocks of complex networks.  Science 298: 824-7.
  * 37. Shen-Orr SS, Milo R, Mangan S, Alon U (2002) Network motifs in the transcriptional regulation network of escherichia coli.  Nat Genet 31: 64-8.
  * 38. Rao CV, Arkin AP (2001) Control motifs for intracellular regulatory networks.  Annu Rev Biomed Eng 3: 391-419.
  * 39. Prill RJ, Iglesias PA, Levchenko A (2005) Dynamic properties of network motifs contribute to biological network organization.  PLoS Biol 3: e343.
  * 40. Csikász-Nagy A, Novák B, Tyson JJ (2008) Reverse engineering models of cell cycle regulation.  Adv Exp Med Biol 641: 88-97.
  * 41. Purcell O, Savery NJ, Grierson CS, di Bernardo M (2010) A comparative analysis of synthetic genetic oscillators.  J R Soc Interface 7: 1503-24.
  * 42. Taylor D, Restrepo JG (2011) Network connectivity during mergers and growth: optimizing the addition of a module.  Phys Rev E Stat Nonlin Soft Matter Phys 83: 066112.
  * 43. Kim JR, Yoon Y, Cho KH (2008) Coupled feedback loops form dynamic motifs of cellular networks.  Biophys J 94: 359-65.
  * 44. Tyson JJ, Chen KC, Novak B (2003) Sniffers, buzzers, toggles and blinkers: dynamics of regulatory and signaling pathways in the cell.  Curr Opin Cell Biol 15: 221–231.
  * 45. Strogatz SH (2000) From Kuramoto to Crawford: exploring the onset of synchronization in populations of coupled oscillators.  Physica D 143: 1–20.
  * 46. Edwards R (2000) Analysis of continuous-time switching networks.  Physica D 146: 165–199.
  * 47. Tyson JJ, Chen K, Novak B (2001) Network dynamics and cell physiology.  Nat Rev Mol Cell Biol 2: 908.
  * 48. Hanahan D, Weinberg R (2000) The hallmarks of cancer.  Cell 100: 57–70.
  * 49. Restrepo JG, Ott E, Hunt B (2005) Onset of synchronization in large networks of coupled oscillators.  Phys Rev E 71: 036151.
  * 50. Restrepo JG, Ott E, Hunt B (2005) Synchronization in large directed networks of coupled phase oscillators.  Chaos 16: 015107.
  * 51. Restrepo JG, Ott E, Hunt B (2006) Emergence of coherence in complex networks of heterogeneous dynamical systems.  Phys Rev Lett 96: 254103.
  * 52. Restrepo JG, Ott E, Hunt B (2006) Emergence of synchronization in complex networks of interacting dynamical systems.  Physica D 224: 114-122.
  * 53. Barabási AL, Oltvai ZN (2004) Network biology: understanding the cell’s functional organization.  Nat Rev Genet 5: 101-13.
  * 54. Pikovsky A, Zaikin A, de la Casa MA (2002) System size resonance in coupled noisy systems and in the ising model.  Phys Rev Lett 88: 050601.
  * 55. Dorogovtsev S (2010) Lectures on Complex Networks.  Oxford University Press.
  * 56. Kaern M, Elston TC, Blake WJ, Collins JJ (2005) Stochasticity in gene expression: from theories to phenotypes.  Nat Rev Genet 6: 451-464.
  * 57. Stein RB, Gossen ER, Jones KE (2005) Neuronal variability: noise or part of the signal?  Nat Rev Neurosci 6: 389-97.
  * 58. Faisal AA, Selen LPJ, Wolpert DM (2008) Noise in the nervous system.  Nat Rev Neurosci 9: 292-303.
  * 59. Gillespie DT (2007) Stochastic simulation of chemical kinetics.  Annu Rev Phys Chem 58: 35-55.
  * 60. Fertig EJ, Danilova LV, Favorov AV, Ochs MF (2011) Hybrid modeling of cell signaling and transcriptional reprogramming and its application in c. elegans development.  Frontiers in Bioinformatics and Computational Biology 2: 77.

## Figure Legends

###### List of Figures

  1. 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 ($r_{\theta}$ black, solid; $r_{x}$ green, dashed; and $r_{\omega}$ blue, dash-dotted) and the bottom-panel shows the evolution of the mean phase $\psi$ (red, solid). (a) Oscillators synchronize and switches stay “on” ($\kappa_{x,x}=11$, $\kappa_{x,\theta}=1.5$, $\kappa_{\theta,x}=1$, $\kappa_{\theta,\theta}=40$, and $\sigma_{\omega}=10$), (b) oscillators freeze (as evidenced by unchanging $\psi$) and switches stay “off” ($\kappa_{x,x}=1$, $\kappa_{x,\theta}=1.5$, $\kappa_{\theta,x}=1$, $\kappa_{\theta,\theta}=40$, and $\sigma_{\omega}=10$), (c) oscillators synchronize and switches oscillate ($\kappa_{x,x}=1$, $\kappa_{x,\theta}=160$, $\kappa_{\theta,x}=0.2$, $\kappa_{\theta,\theta}=42$, and $\sigma_{\omega}=3$), and (d) transitory oscillations in oscillators and switches ($\kappa_{x,x}=0.1$, $\kappa_{x,\theta}=1.4$, $\kappa_{\theta,x}=2$, $\kappa_{\theta,\theta}=1.8$, and $\sigma_{\omega}=10$).
  2. 2 Percentage of simulations in which the qualitative dynamics in Figure 1 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 $\sigma_{\omega}$ for $\kappa_{\theta,x}=0.01$ (solid), $\kappa_{\theta,x}=0.1$ (dotted) and $\kappa_{\theta,x}=1$ (dashed). $\kappa_{x,x}=1<\hat{\kappa}_{x,x}$, $\kappa_{x,\theta}=1.5$, and $\kappa_{\theta,\theta}=40>\hat{\kappa}_{\theta,\theta}$.
  3. 3 Dependence on network size for qualitative states for $\kappa_{\theta,x}=0.01$ and $\sigma_{\omega}=1$. 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).
  4. 4 Dependence on network size for qualitative states for $\kappa_{\theta,x}=1$ and $\sigma_{\omega}=10$. Percentage of simulations with qualitative dynamics plotted as described in Figure 3.
  5. 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)
  6. 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.

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 ($r_{\theta}$ black, solid; $r_{x}$ green,
dashed; and $r_{\omega}$ blue, dash-dotted) and the bottom-panel shows the
evolution of the mean phase $\psi$ (red, solid). (a) Oscillators synchronize
and switches stay “on” ($\kappa_{x,x}=11$, $\kappa_{x,\theta}=1.5$,
$\kappa_{\theta,x}=1$, $\kappa_{\theta,\theta}=40$, and $\sigma_{\omega}=10$),
(b) oscillators freeze (as evidenced by unchanging $\psi$) and switches stay
“off” ($\kappa_{x,x}=1$, $\kappa_{x,\theta}=1.5$, $\kappa_{\theta,x}=1$,
$\kappa_{\theta,\theta}=40$, and $\sigma_{\omega}=10$), (c) oscillators
synchronize and switches oscillate ($\kappa_{x,x}=1$, $\kappa_{x,\theta}=160$,
$\kappa_{\theta,x}=0.2$, $\kappa_{\theta,\theta}=42$, and
$\sigma_{\omega}=3$), and (d) transitory oscillations in oscillators and
switches ($\kappa_{x,x}=0.1$, $\kappa_{x,\theta}=1.4$, $\kappa_{\theta,x}=2$,
$\kappa_{\theta,\theta}=1.8$, and $\sigma_{\omega}=10$).

Figure 2: Percentage of simulations in which the qualitative dynamics in
Figure 1 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 $\sigma_{\omega}$ for $\kappa_{\theta,x}=0.01$ (solid),
$\kappa_{\theta,x}=0.1$ (dotted) and $\kappa_{\theta,x}=1$ (dashed).
$\kappa_{x,x}=1<\hat{\kappa}_{x,x}$, $\kappa_{x,\theta}=1.5$, and
$\kappa_{\theta,\theta}=40>\hat{\kappa}_{\theta,\theta}$. Figure 3: Dependence
on network size for qualitative states for $\kappa_{\theta,x}=0.01$ and
$\sigma_{\omega}=1$. 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). Figure 4: Dependence on network size
for qualitative states for $\kappa_{\theta,x}=1$ and $\sigma_{\omega}=10$.
Percentage of simulations with qualitative dynamics plotted as described in
Figure 3. 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.

### Supplemental Figure Captions

Figure S5. Dependence of dynamics on network size for
$\kappa_{\theta,x}=0.01$. 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).

Figure S6. Dependence of dynamics on network size for $\kappa_{\theta,x}=0.1$.
As for Supplemental Figure S5.

Figure S7. Dependence of dynamics on network size for $\kappa_{\theta,x}=1$.
As for Supplemental Figure S5.


\entry{nisanRonen}{article}{}
\name{author}{2}{}{%
\strng{namehash}{fa6784364edfe871c5fe05cd4186574b}
\strng{fullhash}{fa6784364edfe871c5fe05cd4186574b}
\strng{bibnamehash}{fa6784364edfe871c5fe05cd4186574b}
\strng{authorbibnamehash}{fa6784364edfe871c5fe05cd4186574b}
\strng{authornamehash}{fa6784364edfe871c5fe05cd4186574b}
\strng{authorfullhash}{fa6784364edfe871c5fe05cd4186574b}
\field{labelalpha}{NR01}
\field{sortinit}{N}
\field{sortinithash}{22369a73d5f88983a108b63f07f37084}
\field{labelnamesource}{author}
\field{labeltitlesource}{title}
\field{issn}{0899-8256}
\field{journaltitle}{Games and Economic Behavior}
\field{number}{1}
\field{title}{Algorithmic Mechanism Design}
\field{volume}{35}
\field{year}{2001}
\field{pages}{166\bibrangedash 196}
\range{pages}{31}
\verb{doi}
\verb https://doi.org/10.1006/game.1999.0790
\endverb
\verb{urlraw}
\verb https://www.sciencedirect.com/science/article/pii/S089982569990790X
\endverb
\verb{url}
\verb https://www.sciencedirect.com/science/article/pii/S089982569990790X
\endverb
\endentry
\entry{rachmilevitchFirstPrice}{article}{}
\name{author}{1}{}{%
\strng{namehash}{0e0bffbe81ff3d0fa543da1695499fa0}
\strng{fullhash}{0e0bffbe81ff3d0fa543da1695499fa0}
\strng{bibnamehash}{0e0bffbe81ff3d0fa543da1695499fa0}
\strng{authorbibnamehash}{0e0bffbe81ff3d0fa543da1695499fa0}
\strng{authornamehash}{0e0bffbe81ff3d0fa543da1695499fa0}
\strng{authorfullhash}{0e0bffbe81ff3d0fa543da1695499fa0}
\field{extraname}{1}
\field{labelalpha}{Rac14}
\field{sortinit}{R}
\field{sortinithash}{5e1c39a9d46ffb6bebd8f801023a9486}
\field{labelnamesource}{author}
\field{labeltitlesource}{title}
\field{issn}{0899-8256}
\field{journaltitle}{Games and Economic Behavior}
\field{title}{First-best collusion without communication}
\field{volume}{83}
\field{year}{2014}
\field{pages}{224\bibrangedash 230}
\range{pages}{7}
\verb{doi}
\verb https://doi.org/10.1016/j.geb.2013.11.007
\endverb
\verb{urlraw}
\verb https://www.sciencedirect.com/science/article/pii/S0899825613001607
\endverb
\verb{url}
\verb https://www.sciencedirect.com/science/article/pii/S0899825613001607
\endverb
\keyw{Auctions,Collusion,Repeated games}
\endentry
\entry{rachmilevitchSecondPrice}{article}{}
\name{author}{1}{}{%
\strng{namehash}{0e0bffbe81ff3d0fa543da1695499fa0}
\strng{fullhash}{0e0bffbe81ff3d0fa543da1695499fa0}
\strng{bibnamehash}{0e0bffbe81ff3d0fa543da1695499fa0}
\strng{authorbibnamehash}{0e0bffbe81ff3d0fa543da1695499fa0}
\strng{authornamehash}{0e0bffbe81ff3d0fa543da1695499fa0}
\strng{authorfullhash}{0e0bffbe81ff3d0fa543da1695499fa0}
\field{extraname}{2}
\field{labelalpha}{Rac15}
\field{sortinit}{R}
\field{sortinithash}{5e1c39a9d46ffb6bebd8f801023a9486}
\field{labelnamesource}{author}
\field{labeltitlesource}{title}
\field{issn}{0899-8256}
\field{journaltitle}{Games and Economic Behavior}
\field{title}{Bribing in second-price auctions}
\field{volume}{92}
\field{year}{2015}
\field{pages}{191\bibrangedash 205}
\range{pages}{15}
\verb{doi}
\verb https://doi.org/10.1016/j.geb.2015.06.008
\endverb
\verb{urlraw}
\verb https://www.sciencedirect.com/science/article/pii/S0899825615000962
\endverb
\verb{url}
\verb https://www.sciencedirect.com/science/article/pii/S0899825615000962
\endverb
\keyw{Second-price auctions,Bribing,Collusion,Signaling}
\endentry
\entry{roughgarden2020transaction}{article}{}
\name{author}{1}{}{%
\strng{namehash}{efce6f6d496cdcc1bb07fb47dfecb3aa}
\strng{fullhash}{efce6f6d496cdcc1bb07fb47dfecb3aa}
\strng{bibnamehash}{efce6f6d496cdcc1bb07fb47dfecb3aa}
\strng{authorbibnamehash}{efce6f6d496cdcc1bb07fb47dfecb3aa}
\strng{authornamehash}{efce6f6d496cdcc1bb07fb47dfecb3aa}
\strng{authorfullhash}{efce6f6d496cdcc1bb07fb47dfecb3aa}
\field{extraname}{1}
\field{labelalpha}{Rou20}
\field{sortinit}{R}
\field{sortinithash}{5e1c39a9d46ffb6bebd8f801023a9486}
\field{labelnamesource}{author}
\field{labeltitlesource}{title}
\field{eprinttype}{arXiv}
\field{journaltitle}{CoRR}
\field{title}{Transaction Fee Mechanism Design for the Ethereum Blockchain: An Economic Analysis of {EIP-1559}}
\field{volume}{abs/2012.00854}
\field{year}{2020}
\verb{eprint}
\verb 2012.00854
\endverb
\verb{file}
\verb :roughgarden2020transaction - Transaction Fee Mechanism Design for the Ethereum Blockchain_ an Economic Analysis of EIP 1559.pdf:PDF
\endverb
\endentry
\entry{roughgarden2021transaction}{article}{}
\name{author}{1}{}{%
\list{publisher}{1}{%
{ACM New York, NY, USA}%
\strng{namehash}{efce6f6d496cdcc1bb07fb47dfecb3aa}
\strng{fullhash}{efce6f6d496cdcc1bb07fb47dfecb3aa}
\strng{bibnamehash}{efce6f6d496cdcc1bb07fb47dfecb3aa}
\strng{authorbibnamehash}{efce6f6d496cdcc1bb07fb47dfecb3aa}
\strng{authornamehash}{efce6f6d496cdcc1bb07fb47dfecb3aa}
\strng{authorfullhash}{efce6f6d496cdcc1bb07fb47dfecb3aa}
\field{extraname}{2}
\field{labelalpha}{Rou21}
\field{sortinit}{R}
\field{sortinithash}{5e1c39a9d46ffb6bebd8f801023a9486}
\field{labelnamesource}{author}
\field{labeltitlesource}{title}
\field{journaltitle}{ACM SIGecom Exchanges}
\field{number}{1}
\field{title}{Transaction fee mechanism design}
\field{volume}{19}
\field{year}{2021}
\field{pages}{52\bibrangedash 55}
\range{pages}{4}
\endentry
\entry{reijsbergen2021transaction}{article}{}
\name{author}{6}{}{%
\strng{namehash}{e8f86d4323a86d02453d559ed4b8dc3d}
\strng{fullhash}{0e4c980cdcc825b76054efe0a6e52b6d}
\strng{bibnamehash}{0e4c980cdcc825b76054efe0a6e52b6d}
\strng{authorbibnamehash}{0e4c980cdcc825b76054efe0a6e52b6d}
\strng{authornamehash}{e8f86d4323a86d02453d559ed4b8dc3d}
\strng{authorfullhash}{0e4c980cdcc825b76054efe0a6e52b6d}
\field{labelalpha}{RSMLSP21}
\field{sortinit}{R}
\field{sortinithash}{5e1c39a9d46ffb6bebd8f801023a9486}
\field{labelnamesource}{author}
\field{labeltitlesource}{title}
\field{journaltitle}{2021 IEEE International Conference on Blockchain (Blockchain)}
\field{title}{Transaction Fees on a Honeymoon: Ethereum's EIP-1559 One Month Later}
\field{year}{2021}
\field{pages}{196\bibrangedash 204}
\range{pages}{9}
\verb{file}
\verb :reijsbergen2021transaction - Transaction Fees on a Honeymoon_ Ethereum's EIP 1559 One Month Later.pdf:PDF
\endverb
\endentry
\entry{shi2023what}{inproceedings}{}
\name{author}{3}{}{%
\name{editor}{1}{}{%
given={Yael\bibnamedelima Tauman},
giveni={Y\bibinitperiod\bibinitdelim T\bibinitperiod}}}%
\list{location}{1}{%
{Saarbrücken/Wadern, Germany}%
\list{publisher}{1}{%
{Schloss Dagstuhl - Leibniz-Zentrum für Informatik}%
\strng{namehash}{11e4cbb1a78c7a729dd3c519d132c75b}
\strng{fullhash}{11e4cbb1a78c7a729dd3c519d132c75b}
\strng{bibnamehash}{11e4cbb1a78c7a729dd3c519d132c75b}
\strng{authorbibnamehash}{11e4cbb1a78c7a729dd3c519d132c75b}
\strng{authornamehash}{11e4cbb1a78c7a729dd3c519d132c75b}
\strng{authorfullhash}{11e4cbb1a78c7a729dd3c519d132c75b}
\strng{editorbibnamehash}{482df7d0138a4e484da58b3748958a2b}
\strng{editornamehash}{482df7d0138a4e484da58b3748958a2b}
\strng{editorfullhash}{482df7d0138a4e484da58b3748958a2b}
\field{labelalpha}{SCW23}
\field{sortinit}{S}
\field{sortinithash}{b164b07b29984b41daf1e85279fbc5ab}
\field{labelnamesource}{author}
\field{labeltitlesource}{title}
\field{booktitle}{14th Innovations in Theoretical Computer Science Conference, {ITCS} 2023, January 10-13, 2023, MIT, Cambridge, Massachusetts, {USA}}
\field{series}{LIPIcs}
\field{title}{What Can Cryptography Do for Decentralized Mechanism Design?}
\field{volume}{251}
\field{year}{2023}
\field{pages}{97:1\bibrangedash 97:22}
\range{pages}{-1}
\verb{doi}
\verb 10.4230/LIPIcs.ITCS.2023.97
\endverb
\verb{urlraw}
\verb https://doi.org/10.4230/LIPIcs.ITCS.2023.97
\endverb
\verb{url}
\verb https://doi.org/10.4230/LIPIcs.ITCS.2023.97
\endverb
\endentry
\entry{Stigler1964}{article}{}
\name{author}{1}{}{%
given={George\bibnamedelima J.},
giveni={G\bibinitperiod\bibinitdelim J\bibinitperiod}}}%
\list{publisher}{1}{%
{University of Chicago Press}%
\strng{namehash}{3e27ef99bf72db975408034556b655be}
\strng{fullhash}{3e27ef99bf72db975408034556b655be}
\strng{bibnamehash}{3e27ef99bf72db975408034556b655be}
\strng{authorbibnamehash}{3e27ef99bf72db975408034556b655be}
\strng{authornamehash}{3e27ef99bf72db975408034556b655be}
\strng{authorfullhash}{3e27ef99bf72db975408034556b655be}
\field{labelalpha}{Sti64}
\field{sortinit}{S}
\field{sortinithash}{b164b07b29984b41daf1e85279fbc5ab}
\field{labelnamesource}{author}
\field{labeltitlesource}{title}
\field{issn}{00223808, 1537534X}
\field{journaltitle}{Journal of Political Economy}
\field{number}{1}
\field{title}{A Theory of Oligopoly}
\field{urlday}{6}
\field{urlmonth}{1}
\field{urlyear}{2024}
\field{volume}{72}
\field{year}{1964}
\field{urldateera}{ce}
\field{pages}{44\bibrangedash 61}
\range{pages}{18}
\verb{urlraw}
\verb http://www.jstor.org/stable/1828791
\endverb
\verb{url}
\verb http://www.jstor.org/stable/1828791
\endverb
\endentry
\entry{valiant2007collusion}{article}{}
\name{author}{2}{}{%
\strng{namehash}{8d554140a381b7417c410b457bb79808}
\strng{fullhash}{8d554140a381b7417c410b457bb79808}
\strng{bibnamehash}{8d554140a381b7417c410b457bb79808}
\strng{authorbibnamehash}{8d554140a381b7417c410b457bb79808}
\strng{authornamehash}{8d554140a381b7417c410b457bb79808}
\strng{authorfullhash}{8d554140a381b7417c410b457bb79808}
\field{labelalpha}{VM07}
\field{sortinit}{V}
\field{sortinithash}{afb52128e5b4dc4b843768c0113d673b}
\field{labelnamesource}{author}
\field{labeltitlesource}{title}
\field{title}{Collusion-Resilient Revenue In Combinatorial Auctions}
\field{year}{2007}
\endentry
\entry{maximizingMinerRevenue}{misc}{}
\name{author}{3}{}{%
\strng{namehash}{7d28ac6bcdb1fd7696ffc98a6bb44b94}
\strng{fullhash}{7d28ac6bcdb1fd7696ffc98a6bb44b94}
\strng{bibnamehash}{7d28ac6bcdb1fd7696ffc98a6bb44b94}
\strng{authorbibnamehash}{7d28ac6bcdb1fd7696ffc98a6bb44b94}
\strng{authornamehash}{7d28ac6bcdb1fd7696ffc98a6bb44b94}
\strng{authorfullhash}{7d28ac6bcdb1fd7696ffc98a6bb44b94}
\field{labelalpha}{WSC23}
\field{sortinit}{W}
\field{sortinithash}{4315d78024d0cea9b57a0c6f0e35ed0d}
\field{labelnamesource}{author}
\field{labeltitlesource}{title}
\field{howpublished}{Cryptology ePrint Archive, Paper 2023/283}
\field{note}{\url{https://eprint.iacr.org/2023/283}}
\field{title}{Maximizing Miner Revenue in Transaction Fee Mechanism Design}
\field{year}{2023}
\verb{urlraw}
\verb https://eprint.iacr.org/2023/283
\endverb
\verb{url}
\verb https://eprint.iacr.org/2023/283
\endverb
\endentry
\entry{yao2018incentive}{misc}{}
\name{author}{1}{}{%
given={Andrew\bibnamedelima Chi-Chih},
giveni={A\bibinitperiod\bibinitdelim C\bibinithyphendelim C\bibinitperiod}}}%
\list{publisher}{1}{%
\strng{namehash}{bdf7dbd523c5db366c4ffbc0ea3dc192}
\strng{fullhash}{bdf7dbd523c5db366c4ffbc0ea3dc192}
\strng{bibnamehash}{bdf7dbd523c5db366c4ffbc0ea3dc192}
\strng{authorbibnamehash}{bdf7dbd523c5db366c4ffbc0ea3dc192}
\strng{authornamehash}{bdf7dbd523c5db366c4ffbc0ea3dc192}
\strng{authorfullhash}{bdf7dbd523c5db366c4ffbc0ea3dc192}
\field{labelalpha}{Yao18}
\field{sortinit}{Y}
\field{sortinithash}{fd67ad5a9ef0f7456bdd9aab10fe1495}
\field{labelnamesource}{author}
\field{labeltitlesource}{title}
\field{title}{An Incentive Analysis of some Bitcoin Fee Designs}
\field{year}{2018}
\verb{doi}
\verb 10.48550/ARXIV.1811.02351
\endverb
\keyw{Computer Science and Game Theory (cs.GT),FOS: Computer and information sciences,FOS: Computer and information sciences}
\endentry
\entry{yaish2023uncle}{inproceedings}{}
\name{author}{3}{}{%
\list{location}{1}{%
{Copenhagen, Denmark}%
\list{publisher}{1}{%
{Association for Computing Machinery}%
\strng{namehash}{b64eb14998ca6456a60a2fe55fcf1cd3}
\strng{fullhash}{b64eb14998ca6456a60a2fe55fcf1cd3}
\strng{bibnamehash}{b64eb14998ca6456a60a2fe55fcf1cd3}
\strng{authorbibnamehash}{b64eb14998ca6456a60a2fe55fcf1cd3}
\strng{authornamehash}{b64eb14998ca6456a60a2fe55fcf1cd3}
\strng{authorfullhash}{b64eb14998ca6456a60a2fe55fcf1cd3}
\field{labelalpha}{YSZ23}
\field{sortinit}{Y}
\field{sortinithash}{fd67ad5a9ef0f7456bdd9aab10fe1495}
\field{labelnamesource}{author}
\field{labeltitlesource}{title}
\field{abstract}{We present and analyze an attack on Ethereum 1's consensus mechanism, which allows miners to obtain higher mining rewards compared to their honest peers. This attack is novel in that it relies on manipulating block timestamps and the difficulty-adjustment algorithm (DAA) to give the miner an advantage whenever block races ensue. We call our attack Uncle Maker, as it induces a higher rate of uncle blocks. We describe several variants of the attack. Among these, one that is risk-free for miners. Our attack differs from past attacks such as Selfish Mining, that have been shown to be profitable but were never observed in practice: We analyze data from Ethereum's blockchain and show that some of Ethereum's miners have been actively running a variant of this attack for several years without being detected, making this the first evidence of miner manipulation of a major consensus mechanism. We present our evidence, as well as estimates of the profits gained by attackers, at the expense of honest miners. Since several blockchains are still running Ethereum 1's protocol, we suggest concrete fixes and implement them as a patch for geth.}
\field{booktitle}{Proceedings of the 2023 ACM SIGSAC Conference on Computerand Communications Security (CCS '23)}
\field{isbn}{9798400700507}
\field{series}{CCS '23}
\field{title}{Uncle Maker: (Time)Stamping Out The Competition in Ethereum}
\field{year}{2023}
\verb{doi}
\verb 10.1145/3576915.3616674
\endverb
\verb{file}
\verb :yaish2023uncle - Uncle Maker_ (Time)Stamping Out the Competition in Ethereum.pdf:PDF
\endverb
\keyw{cryptocurrency,blockchain,proof of work,consensus,security}
\endentry
\entry{zhou2023sok}{inproceedings}{}
\name{author}{10}{}{%
\list{location}{1}{%
{Los Alamitos, CA, USA}%
\list{publisher}{1}{%
{IEEE Computer Society}%
\strng{namehash}{beb3b9b64078b19d09dc02044a5e0752}
\strng{fullhash}{adacf7473b0e72c369adc72b50a2b81c}
\strng{bibnamehash}{adacf7473b0e72c369adc72b50a2b81c}
\strng{authorbibnamehash}{adacf7473b0e72c369adc72b50a2b81c}
\strng{authornamehash}{beb3b9b64078b19d09dc02044a5e0752}
\strng{authorfullhash}{adacf7473b0e72c369adc72b50a2b81c}
\field{labelalpha}{ZXECWWQWSG23}
\field{sortinit}{Z}
\field{sortinithash}{96892c0b0a36bb8557c40c49813d48b3}
\field{labelnamesource}{author}
\field{labeltitlesource}{title}
\field{abstract}{Within just four years, the blockchain-based Decentralized Finance (DeFi) ecosystem has accumulated a peak total value locked (TVL) of more than 253 billion USD. This surge in DeFi’s popularity has, unfortunately, been accompanied by many impactful incidents. According to our data, users, liquidity providers, speculators, and protocol operators suffered a total loss of at least 3.24 billion USD from Apr 30, 2018 to Apr 30, 2022. Given the blockchain’s transparency and increasing incident frequency, two questions arise: How can we systematically measure, evaluate, and compare DeFi incidents? How can we learn from past attacks to strengthen DeFi security? In this paper, we introduce a common reference frame to systematically evaluate and compare DeFi incidents, including both attacks and accidents. We investigate 77 academic papers, 30 audit reports, and 181 real-world incidents. Our data reveals several gaps between academia and the practitioners’ community. For example, few academic papers address "price oracle attacks" and "permissonless interactions", while our data suggests that they are the two most frequent incident types (15\% and 10.5\% correspondingly). We also investigate potential defenses, and find that: (i) 103 (56\%) of the attacks are not executed atomically, granting a rescue time frame for defenders; (ii) bytecode similarity analysis can at least detect 31 vulnerable/23 adversarial contracts; and (iii) 33 (15.3\%) of the adversaries leak potentially identifiable information by interacting with centralized exchanges.}
\field{booktitle}{2023 IEEE Symposium on Security and Privacy (SP)}
\field{month}{5}
\field{title}{SoK: Decentralized Finance (DeFi) Attacks}
\field{year}{2023}
\field{pages}{2444\bibrangedash 2461}
\range{pages}{18}
\verb{doi}
\verb 10.1109/SP46215.2023.10179435
\endverb
\verb{file}
\verb :zhou2023sok - SoK_ Decentralized Finance (DeFi) Attacks.pdf:PDF
\endverb
\verb{urlraw}
\verb https://doi.ieeecomputersociety.org/10.1109/SP46215.2023.10179435
\endverb
\verb{url}
\verb https://doi.ieeecomputersociety.org/10.1109/SP46215.2023.10179435
\endverb
\keyw{privacy;smart contracts;ecosystems;finance;decentralized applications;frequency measurement;security}
\endentry
\enddatalist
\endrefsection
\endinput

\appendix

\section{Missing Proofs}
\label[appendix]{app:missing_proofs}

\SimpleJointBound*
\begin{proof}
\[
\begin{split}
& u_{joint}(\mathbf{b}) = \sum_{i=1}^n \left( v_i \alloc(b_i, \mathbf{b}_{-i}) - \burn(b_i, \mathbf{b}_{-i})\right) - \sum_{i=n+1}^{n'} \burn(b_i, \mathbf{b}_{-i}) %\stackrel{\text{\cref{eq:BB_positive}}}{\leq} \\
\stackrel{\text{\cref{eq:BB}}}{\leq} \\
& \sum_{i=1}^n v_i \alloc(b_i, \mathbf{b}_{-i}) \leq \max \mathbf{v} \sum_{i=1}^n \alloc(b_i, \mathbf{b}_{-i}) \stackrel{\text{\cref{eq:feasibility}}}{\leq} \max \mathbf{v}.
\end{split}
\]
\end{proof}

\CoalitionJointBound*

\begin{proof}
The expression for the joint utility under the off-chain agreement $\Omega_{\mathbf{v}}$ takes a similar form to that of the utility of the miner and the $c$ colluding bidders, but also includes the utility of all other bidders.
By the individual rationality condition of \cref{eq:IR_util}, it therefore must be at least as large.

\[
\begin{split}
u_{joint}(\Omega_{\mathbf{v}} ; \mathbf{v})
\stackrel{\text{\cref{eq:joint_util}}}{=}
u_{miner}(\Omega_{\mathbf{v}} ; \mathbf{v}) + \sum_{i=1}^n u_i(\Omega_{\mathbf{v}_i}, \Omega_{\mathbf{v}_{-i}} ; v_i)
\\&
u_{miner}(\Omega_{\mathbf{v}} ; \mathbf{v}) + \sum_{i \in C} u_i(\Omega_{\mathbf{v}_i}, \Omega_{\mathbf{v}_{-i}} ; v_i) + \sum_{i \not \in C} u_i(\Omega_{\mathbf{v}_i}, \Omega_{\mathbf{v}_{-i}} ; v_i)
\\&
\stackrel{\text{\cref{eq:IR_util}}}{\geq}
u_{miner}(\Omega_{\mathbf{v}} ; \mathbf{v}) + \sum_{i \in C} u_i(\Omega_{\mathbf{v}_i}, \Omega_{\mathbf{v}_{-i}} ; v_i).
\end{split}
\]
\end{proof}

\OCAJoint*
\begin{proof}
(OCA-proofness $\implies$ \cref{eq:general_oca})

If \gls{OCA}-proofness holds for any coalition $C$, then it also hold for the grand coalition of all $n$ users.
For any $\mathbf{v}$ and corresponding $\Omega_{\mathbf{v}}$ choose $C = [n]$, then the condition of \cref{eq:oca_cond} takes the form:

$$u_{joint}(\mathbf{v} ; \mathbf{v}) \geq u_{miner}(\Omega_{\mathbf{v}} ; \mathbf{v}) + \sum_{i \in C} u_i(\Omega_{\mathbf{v}_i}, \Omega_{\mathbf{v}_{-i}} ; v_i) = u_{miner}(\Omega_{\mathbf{v}} ; \mathbf{v}) + \sum_{i=1}^n u_i(\Omega_{\mathbf{v}_i}, \Omega_{\mathbf{v}_{-i} } ; v_i)\stackrel{\text{\cref{eq:joint_util}}}{=} u_{joint}(\Omega_{\mathbf{v}} ; \mathbf{v}).$$

(\cref{eq:general_oca} $\implies$ OCA-proofness)

This follows as an immediate corollary of Claim~\ref{clm:coalition_joint_bound}: For any $\mathbf{v}$, corresponding $\Omega_{\mathbf{v}}$ and coalition $C$, it then holds that:

$$u_{joint}(\mathbf{v} ; \mathbf{v}) \geq u_{joint}(\Omega_{\mathbf{v}} ; \mathbf{v}) \stackrel{\text{Claim~\ref{clm:coalition_joint_bound}}}{\geq} u_{miner}(\Omega_{\mathbf{v}} ; \mathbf{v}) + \sum_{i \in C} u_i(\Omega_{\mathbf{v}_i}, \Omega_{\mathbf{v}_{-i}} ; v_i),$$

but this is exactly the condition of \cref{eq:oca_cond}.

\SingleBidderRevenue*
\begin{proof}
By Myerson's lemma (given in~\cref{fct:MyersonLemma}), the allocation rule is monotone and the payment is uniquely determined by it. %In the general case, the winning coalition (which is comprised of the winner and the miner) joint utility is given by $u_{winner}(\mathbf{b}) + u_{miner}(\mathbf{b}) = E[\sum_j v_j 1[a(\mathbf{b}) = j]] - \bar{p}(\mathbf{b}) + \bar{p}(\mathbf{b}) - \bar{\beta}(\mathbf{b}) =  E[\sum_j v_j 1[a(\mathbf{b}) = j]] - \bar{\beta}(\mathbf{b})$ \yotam{Introduce the bar notation for expected payments and burn}.

In the single bidder case, for a bidder with value $v_1$ and bid $b_1$ (i.e., $\mathbf{v} = (v_1)$, with an off-chain agreement of the form $\Omega_{\mathbf{v}} = (b_1)$, no fake bids or omitting by the miner, and the $1$-bidder coalition $C = \{1\}$), the OCA-proofness condition of \cref{eq:oca_cond} takes the form:

$$v_1 \cdot \alloc(v_1) - \burn(v_1) = u_{joint}(v_1 ;  v_1) \geq u_{miner}(b_1 ; v_1) + u_1(b_1 ; v_1) = v_1 \cdot \alloc(b_1) - \burn(b_1).$$

Notice that this is exactly the incentive-compatibility condition for the single bidder, if we replace payment with burn.
Thus, we can follow Myerson's arguments to determine that the burn rule is uniquely determined by the allocation due to the OCA-proofness condition, in exactly the same way as the payment rule is determined due to DSIC. %Thus, since they must be equal in expectation for every bid, and the burn rule is bounded by the payment rule by burn balance (\cref{eq:BB}), we deduce that the payment and burn rule are equal.

%From the critical bid property, $r$ is also the payment of the winning bidder.
%If $r = 0$ then the payment is always $0$, and nothing is burnt.
%If $r > 0$ and there is some $b^* \geq r$ so that $\burn(b^*) < r$, then given a truthful bidder with value $v \in (\burn(b^*), r)$, the item is not allocated.
%But, this bidder and the miner can form a beneficial 1-\gls{OCA} where the bidder bids $b^*$.
%Under this \gls{OCA}, the bidder receives a utility of $v - \pay(b^*)$, and the miner earns $\pay(b^*) - \burn(b^*)$.
%In total, the \gls{OCA} yields a joint utility of $v - \pay(b^*) + \pay(b^*) - \burn(b^*) = v - \burn(b^*) > 0$.
\end{proof}

\end{proof}

\DSICOCA*
\begin{proof}
($\implies$ All DSIC+$1$-OCA-proof mechanisms must be of the form specified)

Given the characterization of \cref{lem:general_oca_char}, the only allocation and burn rules possible are as described. However, since we assume the DSIC property, the payment is uniquely determined by the allocation, and since the allocations all fall within the class of second-price auctions with reserve $r$ allocations, then so do the corresponding payments.

($\implies$ All $2$-nd price auctions with reserve $r$ and constant burn $r$ are DSIC+$1$-OCA-proof)

We know that the $2$-nd price auction with reserve $r$ is \gls{DSIC} \cite{myerson1981optimal}, and we furthermore know that the $1$-OCA-proofness property is satisfied due to the characterization of \cref{lem:general_oca_char}.

\end{proof}

\MMICOCA*
\begin{proof}
($\implies$ All MMIC+$1$-OCA-proof mechanisms must be of the form specified)

Given the characterization of \cref{lem:general_oca_char}, the only allocation and burn rules possible are as described. Regarding payment, consider bidder $i$ with valuation $v_i$, and two sets of bids of other bidders, $\mathbf{v_{-i}}, \mathbf{b_{-i}}$ so that $\tilde{\alloc}(v_i, \mathbf{v}_{-i}) = \tilde{\alloc}(v_i, \mathbf{b}_{-i}) = i$. Then, it must hold that $\pay(v_i, \mathbf{v}_{-i}) - \burn(v_i, \mathbf{v}_{-i}) = \pay(v_i, \mathbf{b}_{-i}) - \burn(v_i, \mathbf{b}_{-i}) $, or otherwise the miner would omit all bids $\mathbf{v}_{-i}$ and add fake bids $\mathbf{b}_{-i}$ (or vice versa). Since the burn is a constant $r$, this implies that $\pay(v_i, \mathbf{v}_{-i}) = \pay(v_i, \mathbf{b}_{-i})$. I.e., the payment is uniquely determined by the winner's bid $v_i$ (and must be higher than the burn $r$ by \cref{eq:BB}). We can thus write $\tilde{\pay}(\mathbf{b}) = f(\max \mathbf{b})$ for some function $f(v) \geq r$.

We now show that $f(v)$ must be monotone whenever $v \geq r$. If not, then there are such $v' > v \geq r$ so that $f(v') < f(v)$. Consider two bidders who bid $v', v$. Without miner intervention, the miner would receive payment $f(v') - r$. If the miner omits the bid $v'$, it receives payment $f(v) - r > f(v') - r$, and so this contradicts MMIC.

($\implies$ All auctions of the form specified are MMIC+$1$-OCA-proof)

We have $1$-OCA-proofness by the characterization of \cref{lem:general_oca_char}. Regarding MMIC, the miner cannot gain by adding fake bids since the burn is constant and the payment only depends on the winning bid: If the added fake bids do not change the winning bid, then the miner revenue does not change, and if they do, since the item is allocated to the highest bidder, this means that one of the fake bids wins, and the miner revenue is then non-positive. The miner cannot gain by omitting bids since the burn is constant and the payment is monotone in the highest bid, and the highest bid can only decrease by omitting bids.

\end{proof}

\homogeneousPayments*

\begin{proof}
For a DSIC mechanism we can write, by Fact~\ref{fct:MyersonLemma}:
\[
\begin{split}
p(\alpha b_i, \alpha \mathbf{b}_{-i})
a(\alpha b_i, \alpha \mathbf{b}_{-i}) \cdot \alpha b_i - \int_0^{\alpha b_i} a(t, \alpha \mathbf{b}_{-i}) dt
\\&
\alpha \cdot a(b_i, \mathbf{b}_{-i}) \cdot b_i - \alpha \int_0^{b_i} a(\alpha \cdot z, \alpha \mathbf{b}_{-i}) dz
\\&
\alpha \left( a(b_i, \mathbf{b}_{-i}) \cdot b_i - \int_0^{b_i} a(z, \mathbf{b}_{-i}) dz \right)
\\&
\alpha p(b_i, \mathbf{b}_{-i}).
\end{split}
\]
\end{proof}

\CTPA*
\begin{proof}
Let $\alpha$ be the constant allocation probability of the CTPA mechanism. Since in the single bidder case the allocation probability is $\alpha$ regardless of the bid, the burn must be $0$ as for any bid $b_1$:
$$\burn(b_1) \stackrel{\text{Claim~\ref{lem:SingleBidderRevenue}}}{=} \pay(b_1) \stackrel{\text{Fact.~\ref{fct:MyersonLemma}}}{=} \alloc(b_1) \cdot b_1 - \int_0^{b_1} \alloc(t) dt = \alpha \cdot b_1 - \int_0^{b_1} \alpha dt = 0.$$
Since the item's allocation probability is always $\alpha$, we can extend the argument of Claim~\ref{clm:simple_joint_bound} to show:
\begin{equation}
\label{eq:alpha_joint_bound}
\begin{split}& u_{joint}(\mathbf{b} ; \mathbf{v}) = \sum_{i=1}^n \left( v_i \alloc(b_i, \mathbf{b}_{-i}) - \burn(b_i, \mathbf{b}_{-i})\right) - \sum_{i=n+1}^{n'} \burn(b_i, \mathbf{b}_{-i}) %\stackrel{\text{\cref{eq:BB_positive}}}{\leq} \\
\stackrel{\text{\cref{eq:BB}}}{\leq} \\
& \sum_{i=1}^n v_i \alloc(b_i, \mathbf{b}_{-i}) \leq \max \mathbf{v} \sum_{i=1}^n \alloc(b_i, \mathbf{b}_{-i}) \leq \alpha \cdot \max \mathbf{v} .
\end{split}
\end{equation}
Thus, for any vector of valuations $\mathbf{v}$,
$$u_{joint}(\mathbf{v} ; \mathbf{v}) \stackrel{\text{\cref{eq:alpha_joint_bound}}}{\leq} \alpha \cdot \max \mathbf{v} = u_{joint}(\max \mathbf{v} ; \mathbf{v}). $$
If the inequality is strict, this violates the OCA-proofness condition of \cref{lem:oca_joint}.
Thus, we must always have $$u_{joint}(\mathbf{v} ; \mathbf{v}) = \alpha \max \mathbf{v},$$
i.e., the item must be allocated with the entire constant probability $\alpha$ to the highest bidder. However, this is a form of a second-price auction (but with limited allocation probability), and because of DSIC this means the payment must increase with the second highest bid, which contradicts MMIC.
\end{proof}

\section{Non-anonymous Deterministic Mechanisms}
\label{app:nonAnon}

We quickly revisit our auction model to allow for non-anonymous mechanisms.
We now have a fixed set of $N$ possible named bidders, out of which some subset $I⊆N$ actually participates in the auction. The allocation, payment, and burn rule $,,$ now operate on the domain $R^I$ and may incorporate the identities of the bidders into the outcome. The DSIC, MMIC, OCA-proofness, and SCP notions remain the same, with the important emphasis that we still allow for the miner to issue fake bids. This is justified by the fact that an agent may hold multiple bidder identities, and so in the analysis perspective it is possible that the miner is in control of virtually any of the bidder identities.

Importantly, \cref{lem:SingleBidderRevenue} ($0$ revenue in the single bidder case), Fact~\ref{fct:MyersonLemma} (Myerson's Lemma), \cref{thm:zero_revenue} ($0$ revenue in the general case) transfer to the non-anonymous case and we subsequently use them.

We now characterize the space of DSIC+MMIC+$1$-OCA-proof non-anonymous mechanisms.
We do so by proving progressive refinements of the space of possible mechanisms. %Thus, we allow the auction functions $(\alloc, \pay, \burn)$ to possibly incorporate the unique ID of the bidder into their operation. For this reason, we consider them not only as functions from $R_+^n$ (where any $n$ bidders with the same bid values are treated in the same way), but rather as functions from $R^I$, where $I \subseteq N$ is some specific set of bidders out of all possible bidders $N$.

\begin{lemma}
\label{lem:non_anonymous}
All \gls{DSIC}+\gls{MMIC}+$1$-\gls{OCA}-proof mechanisms satisfy the following properties:

\begin{enumerate}[leftmargin=*]
\item Each bidder $i$ has a personalized burn rate $r_i$ so that whenever bidder $i$ wins, exactly $r_i$ must be paid and burned.

\item There is at most one unique bidder $i^*$ with $r_{i^*} \neq \infty$.

\item For any given set of bidders $I$, the item is allocated only if $i^* \in I$, and the allocation rule is of the form $\tilde{a}(\mathbf{b}) = \begin{cases} i^* & b_{i^*} \geq r_{i^*}, \\
\emptyset & Otherwise.\end{cases}$
\end{enumerate}
\end{lemma}

\begin{proof}
We give a proof sketch, as the proof resembles the one for the anonymous case.

We prove each property separately.
\begin{enumerate}[leftmargin=*]
\item By \cref{lem:SingleBidderRevenue} and Fact~\ref{fct:MyersonLemma}, each bidder $i \in N$ has some $r_i$ so that in case $i$ is the single bidder, $\alloc(v_i) = \mathds{1}_{v_i \geq r_i}, \tilde{\pay}(v_i) = \tilde{\burn}(v_i) = r_i$. Thus, it must hold that whenever $i\in I$ and for any $v_i$ so that $\tilde{\alloc}(v_i, \mathbf{v}_{-i}) = i$, it also holds that $\tilde{\burn}(v_i, \mathbf{v}_{-i}) = r_i$. Since there is $0$ revenue, it also holds that $\tilde{\pay}(v_i, \mathbf{v}_{-i}) = r_i$.

\item Assume towards contradiction that there are two different bidders $i,j$ with $r_i \neq \infty, r_j \neq \infty$. Consider $I = \{i, j\}$. For $\mathbf{v}$ so that $v_i = r_i + 2, v_j = r_j + 1$, it must hold that $\tilde{\alloc}(\mathbf{v}) = i$, otherwise there is a $1$-OCA where the miner omits $v_j$. But since we already concluded that bidder $i$ pays $r_i$ in this case, by the critical bid property, for $\mathbf{v}'$ where $v'_i = r_i, v'_j = r_j + 1$, it also holds that $\tilde{\alloc}(\mathbf{v}) = i$. But this means there is a $1$-OCA where the miner omits $v'_i$.

\item No other bidder besides $i^*$ can be allocated because all have infinite burn rates if allocated. Thus, if $i^* \not \in I$, then no bidder is ever allocated. If $i^* \in I$, and $v_{i^*} > r_{i^*}$, assume towards contradiction that $\tilde{\alloc}(\mathbf{v}) \neq i^*$, then it must be $\tilde{\alloc}(\mathbf{v}) = \emptyset$. But then there is a $1$-OCA where the miner omits all other bids and the joint utility is $v_{i^*} - r_{i^*} > 0$. If  $v_{i^*} < r_{i^*}$, assume towards contradiction that $\tilde{\alloc}(\mathbf{v}) = i^*$. Since the burn must be at most $v_{i^*}$, there is a $1$-OCA from the single $i^*$ bidder case where the miner adds the fake bids $\mathbf{v}_{-i}$ to reduce the burn.
\end{enumerate}

\end{proof}

It is then straightforward to show that the class of mechanisms we identified cannot be further refined, and in fact satisfies a stronger property:

\begin{lemma}
The class of mechanisms characterized in \cref{lem:non_anonymous} is \gls{DSIC}+MMIC+SCP.
\end{lemma}
\begin{proof}
Consider a mechanism contained in the class.
It is \gls{DSIC}, as per the definition of the class in \cref{lem:non_anonymous}, the mechanism satisfies the monotone allocation and critical bid payment conditions of Myerson's Lemma~\cref{fct:MyersonLemma}.
The MMIC property is satisfied as the payment and burn rules are equal, implying that the miner always receives $0$ revenue, thereby preventing the miner from being able to unilaterally increase its revenue without potentially colluding with users.

Finally, note that no beneficial collusion exists: the individual burn rate is fixed for each user, and cannot change as dependent on the bids of the other transactions included in the block, or the identities of the transactions' creators.

SCP is satisfied since if the coalition $C$ does not contain $i^*$, the aggregate utility of the coalition is $0$ with truthful reports, while non-positive under any manipulation. If the coalition does contain $i^*$, the SCP condition becomes the DSIC condition for agent $i^*$, which is satisfied as the mechanism is DSIC.

\end{proof}

\end{document}


# The Interactive Modeling of a Binary Star System

Canberk Soytekin
Çakabey Schools
Izmir, Turkey
<EMAIL_ADDRESS>
Ayşe Pelin Dedeler
Çakabey Schools
Izmir, Turkey
<EMAIL_ADDRESS>

###### Abstract

An interactive binary star system simulation was developed to be showcased on
an educational platform. The main purpose of the project is to provide insight
into the orbital mechanics of such star systems with the help of a three-body
simulation. The initial simulation script was written in the Python
programming language with the help of the VPython addition. The custom-made
models were created on Blender and exported. For the final implementation of
the simulation on the Godot game engine, the Python code was converted into
GDScript and the Blender models were re-textured.

_K_ eywords Three-Body Problem $\cdot$ Binary Stars $\cdot$ Interactive
Learning

##  1 Introduction

The three-body problem was posed by Newton in Principia [1] , where he
considered the motion of the Earth and the Moon around the Sun. Since then,
the three-body problem has been one of the most well-known problems in
dynamical astronomy [2].Extensive theoretical and numerical research has been
dedicated to the study of this attractive problem; however, for simplicity’s
sake, the majority of research has focused on the circular constrained problem
[3]. The development of high-speed computers has increased interest in
simulating the general three-body problem. Applications of the three-body
problem in astronomy are essential for determining how three stars, a star
with a planet that has a moon, or any other set of three celestial objects can
maintain a stable orbit. Three-body problem based foundational models of the
gravitational influence of the Sun on Earth and the Moon were investigated
while constructing significant space missions like the James Webb Space
Telescope [4]. The Hubble Space Telescope, the Chandra X-Ray Observatory, the
SPITZER Space Telescope, and the Kepler Space Telescope are some of the most
noteworthy discovery missions that NASA has conducted; each of these missions
have used a distinct approach to the three-body problem[5].

Figure 1: The TOI-1338 Binary Star System [6]

Scientists are increasingly favoring active learning tools [7]. Simulations
may not only serve to motivate students, but also to enhance their intuitive
knowledge of abstract physics problems.[8], To achieve our goal of turning the
three-body problem into an interactable topic for students, researchers, and
enthusiasts, we have created a simulated system consisting of two stars and
one observatory satellite which are in elliptical orbits determined by a
differentially simulated application of classical mechanics principles. The
initial parameters of the observatory satellite: position (x, y, z) and
velocity (x, y, z) vectors can be given by the user to initiate the
simulation.

##  2 The Scientific Groundwork

The scientific bases of the project were obtained from Orbital Mechanics by
Vladimir Chobotov [9], Elements of Newtonian Mechanics by Jens M. Knudsen and
Poul G. Hjorth [10], and the MIT 8.01SC Classical Mechanics Course.

$m_{2}\frac{v_{0}^{2}}{R_{12}}=G\frac{m_{1}m_{2}}{R_{12}^{2}},\quad\centering
v_{0}^{2}=G\frac{m_{1}}{R_{12}},\quad
v_{0}=\sqrt{\frac{Gm_{1}}{R_{12}}};\quad\@add@centering$ (1)

$\overrightarrow{p_{02}}=\left[\begin{array}[]{c}0\\\
\mathrm{~{}m}_{2}\cdot\mathrm{v}_{0}\end{array}\right],\quad\sum\vec{p}=0,\quad\overrightarrow{p_{01}}-\overrightarrow{p_{02}}$
(2)

Conservation of momentum and Kepler’s First Law for celestial bodies were
utilized in determining the initial velocities of the binary stars as well as
satisfying the orbital prerequisites.

$\overrightarrow{r_{H0}}=R\textsubscript{12}\left[\begin{array}[]{l}r_{x}\\\
r_{y}\\\
r_{z}\end{array}\right],\quad\overrightarrow{v_{H0}}=\mathrm{v}_{f}\left[\begin{array}[]{l}v_{x}\\\
v_{y}\\\ v_{z}\end{array}\right]$ (3)

The scalar properties of the satellite are mass (m0) and initial speed (vf).
The vectoral properties are its position vector and the components of its
velocity.

Figure 2: The diagram of the initial conditions.

##  3 Python Simulation

In the mathematical modeling of the scientific bases of the project, a
differential calculation method was utilized. In the Python simulation, the
positions and the momenta of the stars and the satellite were updated
approximately every 200 microseconds (dt = 200 µs).

Throughout each dt time interval, the forces on the bodies and their momenta
were assumed to be constant.

The momenta were updated differentially according to the following equations:

$\overrightarrow{p_{new}}=\overrightarrow{p_{old}}+\overrightarrow{\Delta
p},\quad\overrightarrow{\Delta p}=\int_{t_{0}}^{t_{0}+\Delta t}\vec{F}dt$ (4)

$\lim_{\Delta t\rightarrow 0}\quad\overrightarrow{\Delta p}=\vec{F}\cdot\Delta
t,\quad\vec{p}_{new}^{\prime}=\overrightarrow{p_{old}}+\vec{F}\cdot\Delta t$
(5)

The positions were updated differentially according to the following
equations:

$\vec{r}_{new}^{\prime}=\overrightarrow{r_{old}}+\overrightarrow{\Delta
r},\quad\overrightarrow{\Delta r}=\int_{t_{0}}^{t_{0}+\Delta
t}\vec{v}dt,\quad\vec{p}=m\vec{v},\quad\vec{F}=\frac{d\vec{p}}{dt}$ (6)

$\lim_{\Delta t\rightarrow 0}\quad\overrightarrow{\Delta
r}=\frac{\vec{p}}{m}\cdot\Delta
t,\quad\vec{r}_{new}^{\prime}=\overrightarrow{r_{old}}+\frac{\vec{p}}{m}\cdot\Delta
t$ (7)

Figure 3: The Python Simulation

Figure 4: Time-dependent diagram of the  $\overrightarrow{F_{net}}(t)$ vector
of the satellite.

The vector analysis and visual Modeling of the bodies were done with the help
of the VPython module.

Two varieties of code were written for the Python simulation: a visual code
and an analysis code. In the visual code, the script runs at a rate of 500 (a
maximum 500 operations per second, a property of VPython) and is focused on
the visual depiction of the simulation; in the analysis code, an intensive
vectoral analysis module utilizing numpy and matplotlib was integrated into
the while loop in addition to the visualization module. Even through the
analysis code runs at 10 times the rate (a maximum of 5000 operations per
second) of the visual code, the increased workload causes it to run slower
compared to its less intensive counterpart. Both simulation scripts are
available on GitHub.

https://github.com/KameSanKaro/The-Interactive-Modeling-of-a-Binary-Star-
System---Supplementary-Code.git

Figure 5: The satellite model.

##  4 Godot Implementation

In the Godot implementation of the simulation, the simulation code (visual
code) was converted from Python to GDScript for increased Godot compatibility.
In addition to the translation of the primary simulation code, various signal
scripts were implemented for navigability measures, such as menu screens.

The transfer of variables (initial parameters and the score) between scenes
were conducted with the help of the Autoload property of Godot. The simulation
objects were implemented as rigid bodies, and their kinematic and inertial
properties were made compatible to those of Godot. In the Godot implementation
process, the code was changed from momentum-based to velocity-based because of
Godot not being able to handle the large masses of the objects in the
calculations:

$\overrightarrow{v_{new}}=\overrightarrow{v_{old}}+\overrightarrow{\Delta
v},\quad\overrightarrow{\Delta v}=\int_{t_{0}}^{t_{0}+\Delta t}\vec{a}dt$ (8)

$\lim_{\Delta t\rightarrow 0}\quad\overrightarrow{\Delta v}=\vec{a}\cdot\Delta
t,\quad\vec{v}_{new}^{\prime}=\overrightarrow{v_{old}}+\vec{a}\cdot\Delta t$
(9)

$\vec{r}_{new}^{\prime}=\overrightarrow{r_{old}}+\overrightarrow{\Delta
r},\quad\overrightarrow{\Delta r}=\int_{t_{0}}^{t_{0}+\Delta t}\vec{v}dt,$
(10)

$\lim_{\Delta t\rightarrow 0}\quad\overrightarrow{\Delta
r}={\vec{v}}\cdot\Delta
t,\quad\vec{r}_{new}^{\prime}=\overrightarrow{r_{old}}+{\vec{v}}\cdot\Delta t$
(11)

The simulation was changed from 2D in Python to 3D in Godot.

The stars and the satellite were modeled in Blender and exported in the .obj
format. All textures on the exported models were in the Principled BDSF
format.

The models were re-textured, and the glow of the stars were re-added with a
WorldEnvironment node, under the main simulation scene (System.tscn). The
corona effect of the stars was added with the help of the PulseGlow shader for
Godot (a shader add-on by @fbcosentino).

Figure 6: The Stars modeled in Blender.

Figure 7: Close-ups of the Godot implemented star models.

###  4.1 Scoring Function

In order to ascertain the relative success of each simulation trial, a scoring
function was added onto the final Godot simulation. The scoring function adds
the increment value (one minus the ratio of the distance to the closest star
over 1000) in each delta time frame, where delta is a parameter of the process
procedure of the Godot game engine. If the satellite “has made observations of
both stars” (has entered within a radius of 200 of each star), the overall
score is multiplied by two. The scoring function terminates at the end of each
simulation run.

$\textbf{Score}=[(r_{1}flag^{\wedge}r_{2}flag).2]\cdot\sum_{i=delta}^{SimulationTime}1-\frac{r_{close}}{1000}$
(12)

##  5 Conclusion

In conclusion, a realistic and interactive binary star simulation based on
classical mechanics was developed. The simulation aims to provide an intuition
of orbital mechanics in the form of a game; to that end, easily navigable menu
screens were developed.

Figure 8: The "Start" Screen

Figure 9: The "Parameters" Screen and an ongoing simulation.

Figure 10: The “Crash” game-over screen and “Escaped orbit!” game-over
screens.

In the simulation, the aim is to get Helios-3, an observatory satellite, in an
orbit as close as possible to each of the stars to make the necessary
scientific observations. The player provides the initial parameters for the
simulation (the velocity and position vectors of the satellite). In each
trial, the “usefulness of each run to the scientists” is evaluated with the
help of a scoring function. After each run, the player is able to see their
score, the replay of their run and why their run has ended.

## References

  * [1] Isaac Newton.  Philosophiæ Naturalis Principia Mathematica.  Londini Apud G., J. Innys Second, 1726.
  * [2] Steven N. Shore.  Celestial mechanics.  Encyclopedia of Physical Science and Technology (Third Edition), pages 527–540, 2003.
  * [3] Allan Kardec de Almeida Junior and Antonio Fernando Bertachini de Almeida Prado.  Comparisons between the circular restricted three-body and bi-circular four body problems for transfers between the two smaller primaries.  Scientific Reports, 12(1):4148, Mar 2022.
  * [4] Clampin M. et al. Gardner J.P., Mather J.C.  The james webb space telescope.  Space Sci Rev 123, pages 485––606, 2006.
  * [5] : Z E Musielak and B Quarles.  The three-body problem.  Rep. Prog. Phys. 77 065901, 2014.
  * [6] Yvette Smith.  Discovering circumbinary star systems, 2021.  [Online; accessed October 2, 2022].
  * [7] Wen-Chung Shih, Shian-Shyong Tseng, and Chao-Tung Yang.  Deployment of interactive games in learning management systems on cloud environments for diagnostic assessments.  In Maiga Chang, Wu-Yuin Hwang, Ming-Puu Chen, and Wolfgang Müller, editors, Edutainment Technologies. Educational Games and Virtual Reality/Augmented Reality Applications, pages 492–496, Berlin, Heidelberg, 2011. Springer Berlin Heidelberg.
  * [8] Kuo-En Chang, Yu-Lung Chen, He-Yan Lin, and Yao-Ting Sung.  Effects of learning support in simulation-based physics learning.  Computers & Education, 51(4):1486–1498, 2008.
  * [9] Vladimir A. Chobotov.  Orbital Mechanics, Third Edition.  American Institute of Aeronautics and Astronautics, Inc., 2002.
  * [10] Knudsen J M and Hjorth P G.  Elements of Newtonian Mechanics: Including Nonlinear Dynamics.  (Heidelberg: Springer), 2000.