iSciencell

OPEN ACCESS

ArticleA design principle for posttranslational chaotic
oscillators
Hiroto Q.
Yamaguchi, Koji L.
Ode, Hiroki R.
Uedauedah-tky@umin.ac.jp

Highlights
Two substrates withreversible two-site
phosphorylation can
exhibit chaos behavior
The chaos does not
require autocatalysis or
allosteric regulation of
enzymes
The chaos is a result of thecoupling of two substrates
via enzyme availability

Yamaguchi et al., iScience 24,101946
January 22, 2021 ª 2020 The
Authors.
https://doi.org/10.1016/
j.isci.2020.101946

iSciencellOPEN ACCESS

ArticleA design principle for posttranslational
chaotic oscillators
Hiroto Q. Yamaguchi,1,2,3 Koji L. Ode,1,4 and Hiroki R. Ueda1,4,5,*
Summary

Chaos behavior has been observed in various cellular and molecular processes.
Here, we modeled reversible phosphorylation dynamics to elucidate a designprinciple for autonomous chaos generation that may arise from generic enzymatic reactions. A comprehensive parameter search demonstrated that the reaction system composed of a set of kinases and phosphatases and two substrates
with two modification sites exhibits chaos behavior. All reactions are describedaccording to the Michaelis-Menten reaction scheme without exotic functions being applied to enzymes and substrates. Clustering analysis of parameter sets thatcan generate chaos behavior revealed the existence of motif structures. Thesechaos motifs allow the two-substrate species to interact via enzyme availability
and constrain the two substrates’ dynamic changes in phosphorylation status
so that they occur at different timescales. This chaos motif structure is found inseveral enzymatic reactions, suggesting that chaos behavior may underlie cellular
autonomy in a variety of biochemical systems.
Introduction
The spontaneous activation of chemical reactions is the basis of cellular autonomy. Biological oscillatorsplay an important role in the control of the autonomous and temporal activities that occur in cells. Examples
include autonomous rhythmicity observed in the circadian cellular clock process (Nagoshi et al., 2004), cell
cycle progression (Murray and Kirschner, 1989), and rhythmic neuronal firing (Izhikevich, 2007). Such activities can span timescales from milliseconds to hours. Biological oscillators are also thought to play animportant role in autonomous pattern formation from body segmentation in mammals (cm scale) (Palmeirim et al., 1997) and the positioning of cell division machinery inside unicellular prokaryotes (mm scale)(Raskin and de Boer, 1999). These oscillators are often modeled as feedback loops composed of genetic,
biochemical, and cellular events that exhibit limit cycle oscillation over a wide period range (Goldbeter,
2002). Biological systems are composed of multiple oscillators, which in many cases interact with eachother. For example, it was recently demonstrated that the segmentation clock period is modulated byintercellular coupling of genetic oscillators (Yoshioka-Kobayashi et al., 2020). Such coupling also occurs
within cells among oscillators that drive various physiological functions, e.g., circadian clocks influencecell cycle progression (Matsuo et al., 2003; Yang et al., 2010) and, in turn, cell cycle phase influences thecircadian clock (Bieler et al., 2014). Such coupling of multiple oscillators may help produce robust oscillation (Droin et al., 2019; Yan and Goldbeter, 2019b).
In certain conditions, coupled oscillators may exhibit more complex dynamics including multirhythmicity
and chaos behavior (Gerard and Goldbeter, 2012; Yan and Goldbeter, 2019a). Indeed, chaos behavior
may play a fundamental role in the functioning of biological systems; it has been suggested that autonomous and complex systems (such as living systems) converge on the edge of chaos (Kauffman, 1992; Langton, 1990). Autonomous chaos behavior has been observed and modeled extensively in various biologicallayers; early examples include the peroxidase enzymatic reaction (Olsen and Degn, 1977), neuronal firing(Aihara et al., 1985; Skarda and Freeman, 1987), and cAMP signaling in Dictyostelium (Durston, 1974; Martiel and Goldbeter, 1985). Chaos behavior continues to be recognized for its importance in modern
biology. For example, a recent modeling study of the sleep-wake cycle suggested that the transition between chaos behavior and regular oscillation may underlie the state alteration from sleep-associated to
wake-associated neuronal firing patterns (Rasmussen et al., 2017). In such studies, mathematical models
are used to understand and predict how complex and autonomous behavior occurs when cellular and molecular activities are modulated. However, it is often difficult to understand how chaos behavior arises from

1Department of Systems

Pharmacology, Graduate
School of Medicine, theUniversity of Tokyo, Hongo
7-3-1, Bunkyo-ku, Tokyo,
113-0033, Japan
2Center for Information and

Neural Networks, National
Institute of Information andCommunications
Technology, 1-4 Yamadaoka,
Suita, Osaka 565-0871, Japan
3Graduate School of Frontier

Biosciences, Osaka
University, 1-3 Yamadaoka,
Suita, Osaka 565-0871, Japan
4Laboratory for Synthetic

Biology, RIKEN Center for
Biosystems Dynamics
Research, 1-3 Yamadaoka,
Suita, Osaka 565-0871, Japan
5Lead Contact*Correspondence:uedah-tky@umin.ac.jp
https://doi.org/10.1016/j.isci.
2020.101946

iScience 24, 101946, January 22, 2021 ª 2020 The Authors.
This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).

1ll

OPEN ACCESS

a biochemical reaction system because the biological networks that give rise to chaos behavior can involve
multiple gene expression networks, cell-cell interactions, and the complex regulation of molecules.
In modern research, a synthetic biology approach is increasingly being employed to gain mechanistic insights into the autonomous generation of complex behaviors from biological components. In thisapproach, the aim is to design complex responses from relatively simple model systems. Protein networks
shaped by enzymatic reactions and protein-protein interactions have recently been studied using a synthetic biology approach, leading to the de novo design of protein logical circuits (Chen et al., 2020; Fink
et al., 2019; Gao et al., 2018). Complex responses have also been observed in nuclear-free cytosolic extractsystems, where protein networks, rather than genetic circuits, play a role to elicit autonomous responses(Pomerening et al., 2003, 2005). Therefore, the design of protein networks in cells may improve understanding of how the complex responses exhibited by individual cells originated. Mathematical models can facilitate the design of synthetic protein-based circuits for autonomous activity including chaos. These models
must be detailed enough to map to measurable biochemical parameters but also contain a small enoughnumber of components to facilitate experimental design.
In the present study, we use modeling of biological autonomous chaos to interrogate the general designprinciples of biological chaos that can be realized in a set of ubiquitous molecular systems. Reversible
phosphorylation is a well-understood mechanism by which protein activity is regulated. A typical modelof reversible phosphorylation assumes a set of antagonizing enzymes (a kinase and a phosphatase) anda substrate (Conradi and Shiu, 2018; Kholodenko, 2006). The enzymes catalyze the phosphorylation anddephosphorylation of the substrate following the Michaelis-Menten reaction scheme. Despite its simplicity
and generality, this reaction scheme can generate various complex behaviors; in particular, reversible
phosphorylation with two or more modification sites in one-substrate molecule can exhibit biostability(Markevich et al., 2004), multistability (Thomson and Gunawardena, 2009), and autonomous oscillation(Conradi et al., 2019; Jolley et al., 2012; Suwanmajo and Krishnan, 2015). In addition, the spatiotemporal
dynamics of this system include traveling waves (Markevich et al., 2006) and the Turing pattern (Sugai
et al., 2017). The generic applicability of the reversible phosphorylation model to the modeling of various
types of posttranslational modification enables researchers to investigate a typical molecular mechanismthat produces complex behaviors from a combination of simple enzymatic reactions.
Specifically, we demonstrate here that chaos behavior can arise from a reversible phosphorylation model
entirely described by the Michaelis-Menten reaction scheme. The system does not assume typical allosteric
regulation of enzymes in which the enzymes’ reaction rates are actively modulated by the substrate, explicit
autocatalytic activities, or any other exotic control of the enzymes. The complexity of our model system
arises from its composition, which includes two-substrate species each of which has two modification sites.
We found that typical combinations of reaction parameters create chaos motifs that allow each of the two
substrates to behave as an oscillator and to couple with each other through the availability of sharedenzymes. This motif structure can be found in several well-known kinases. Therefore, our model will help
researchers to identify and perhaps design possible enzymatic reaction components that can induce
autonomous complex outputs in cells.

Results
A coupled dual-phosphorylation system shows chaos behavior
The simplest type of reversible phosphorylation scheme is a system with a single substrate that has one
phosphorylation site (Figure S1A). However, this scheme can be modeled as a one-dimensional system
by applying the Michaelis-Menten approximation, and thus, oscillation and chaos are impossible (Jolley
et al., 2012; Strogatz, 1994). Therefore, we next consider a similar reversible phosphorylation system but
single substrate contains two phosphorylation sites. In this system, we previously found that reversible
two-site phosphorylation in a single substrate system can have a limit cycle oscillation (Jolley et al.,2012). We assessed whether the same system could exhibit chaos behavior. The system consists of a single
substrate with two phosphorylation sites, which are phosphorylated and dephosphorylated by a kinase anda phosphatase, respectively (Figure S1A). All the reactions follow the Michaelis-Menten reaction scheme.
Parameters were randomly sampled from an exponential distribution bounded between 1 and 1000 min 1
for reaction rate constants (k) and between 0.01 and 1000 mM for binding constants (Km). Note that a low Km
value indicates strong binding between an enzyme and the substrate. The parameter range corresponds toreasonable values typically used in other simulation studies (Chickarmane et al., 2007; Jolley et al., 2012; Liu

2

iScience 24, 101946, January 22, 2021

iScience
ArticleiSciencell

ArticleOPEN ACCESS

et al., 2011; Markevich et al., 2004; Qiao et al., 2007; Shankaran et al., 2009; Sugai et al., 2017) and within the
distribution of kcat and Km values observed in various enzymes (Bar-Even et al., 2011). Through the numerical integration of the system with 10,000,000,000 randomly sampled parameter sets, we identified
12,636,210 parameter sets that exhibit limit cycle oscillation but did not find any parameters that exhibit
chaos behavior, suggesting that the single-substrate system is not complex enough to include chaos
behavior.
We therefore considered a dual-substrate system consisting of two types of substrate, named Sa and Sb(Figure 1A). Each substrate has two phosphorylation sites that are modified reversibly by a shared kinase(E) and a shared phosphatase (F). Sa has four phosphorylation states, named Sa00, Sa01, Sa10, and Sa11; Sb
also has four phosphorylation states, named Sb00, Sb01, Sb10, and Sb11. The states Sa00 and Sb00 are unphosphorylated substrate states, whereas Sa01, Sa10, Sb01, and Sb10 are substrate states with single phosphorylation at either one of the modification sites. The remaining states, Sa11 and Sb11, are substrate states withphosphorylation at both of the two modification sites.
This dual-substrate system showed chaos behavior with some parameter sets (hereafter called ‘‘chaosparameter sets’’). With chaos parameter sets, the system never converged to a stable state or limit cycle
oscillation (Figure 1B). The chaos behavior was further validated by showing that small changes in the initialstate could lead to a variety of different states after the system had evolved to a certain time (Figure 1C).
The power spectrum of chaos orbitals showed a broad spectrum rather than a few peaks with strong power
intensity (Figure 1D). These features distinguish the chaos parameter sets from parameters resulting in limit
cycle oscillations (hereafter called ‘‘oscillation parameter sets’’) (Figures S1B–S1D).

Conserved parameter motifs in the chaos parameter sets
To ascertain whether a specific parameter region was related to chaos behavior in the parameter space,chaos parameter sets were collected by a random parameter search within the same distribution of kcat
and Km values described above. This random search revealed 4,818 chaos and 286,255 oscillation parameter sets from 13,269,134,692 randomly generated parameter sets (0.0000363% for chaos and 0.00216%
for oscillation parameter sets) (Figure 2A).
To identify a typical motif structure of the chaos parameter sets, hierarchical clustering analysis was performed on the collected chaos parameter sets. Seven clusters were defined by setting a threshold to theclustering dendrogram (Figure 2B). To elucidate the characteristic combination in each cluster, we identified the fast reaction rate based on the criterion that the probability of the corresponding k parameter exceeds a certain threshold (0.2; for details, see Transparent methods) at the highest parameter bin(1000 min 1). Similarly, we identified the stable substrate-enzyme complex based on the criterion thatthe probability of corresponding Km exceeds the threshold at the lowest parameter bin (0.01 mM). Figure 2C(upper section) shows a schematic view of the chaos parameter sets classified as cluster 1. This scheme
highlights the positions of large k values (fast reaction speeds) and small Km values (stable enzyme-substrate complex formation). The parameter distributions of clusters 2–7 are shown in Figure S2. The histograms in each panel indicate the parameter distributions (Figures 2C and S2).
We found that all the clusters included a shared motif structure in at least one part of the phosphorylation/dephosphorylation pathway (as highlighted by yellow coloring trimmed with green or orange in Figures 2C
and S2 [discussed later]). The shared motif structure shown in Figure 2D is characterized by a fast reactionthat produces stable binding of the substrate (Spost) with either the kinase or phosphatase that is not thesame enzyme catalyzed in the fast reaction (i.e., if the fast reaction is a phosphorylation reaction, the product is stably bound with a phosphatase not a kinase.) For example, in cluster 1 (Figure 2C), a kinase reaction
with a high ka4 value produces Sa11, which involves stable binding with the phosphatase as indicated by the
low Kma7 value. In the same way, the other motif is composed of a phosphatase reaction with a high ka5
value that produces Sa00, which involves stable binding with the kinase as indicated by a low Kma2 value.

Motif arrangements necessary for chaos behavior
The motifs can be categorized into two types (named motif A and B, respectively) by considering the symmetry of a state transition diagram (Figure 3A). In motif A, a fast reaction is assigned such that a fully phosphorylated (S00) or fully dephosphorylated product (S11) is produced. On the other hand, motif B is
composed of a fast reaction that produces a single phosphorylated substrate (S01 or S10). S00 and S11

iScience 24, 101946, January 22, 2021

3

iSciencell

ArticleOPEN ACCESS

A

B

C

D

Figure 1. A coupled dual-phosphorylation system shows chaos behavior(A) A reaction network comprising a two-substrate (i.e., Sa and Sb) system. Each substrate undergoes a reversible
phosphorylation cycle at two residues. The phosphorylation (green arrows) and dephosphorylation (blue arrows)
reactions are catalyzed by a shared kinase and phosphatase, respectively. All the reactions follow a Michaelis-Menten
reaction scheme characterized by parameters for reaction speed k and enzyme-substrate affinity Km.
(B) An example time course of chaotic phosphorylation dynamics.
(C) The example chaotic trajectory from t = 0 min to t = 5,000 min is projected to Sa00-Sa01-Sa10 space (green line). Red dots
indicate the phosphorylation states at t = 5,000 min obtained from 200 runs of the simulation with slightly different initialconditions.
(D) Power spectrum of the example chaotic trajectory analyzed by fast Fourier transform (FFT).
See also Figure S1.

can be converted by converting the kinase and phosphatase; thus, they have a symmetrical relationship. S01
and S10 are also symmetrical because they can be converted by converting the position of two phosphorylation residues. However, S00 cannot be converted into S01 or S10 by a symmetrical conversion because S00

4

iScience 24, 101946, January 22, 2021

iSciencell

ArticleA

OPEN ACCESS

C

B

D

Figure 2. Shared motifs in chaotic parameter sets(A) Workflow for chaotic motif identification.
(B) Clustering of chaotic parameter sets with Ward’s algorithm. The chaotic parameter sets were standardized byconsidering phosphor-site symmetries.
(C) Schematic representation and parameter histograms for the largest cluster. Reaction constants (k) for which the
distribution exceeds the probability 0.2 (red dotted lines) at the highest value are denoted as fast reactions. Michaelis
constants (Km) for which the distribution exceeds the probability 0.2 at the lowest value are denoted as high-affinity
enzyme-substrate bindings. The fast reactions are shown as thick arrows, while the high-affinity bindings are indicated byattaching the kinase (E) or phosphatase (F) to the edge of the corresponding substrate species.
(D) A shared design motif found in all clusters. A fast reaction produces a substrate species that has high affinity with either
E or F.
See also Figure S2.

binds to and is a substrate for a kinase only. In contrast, S01 and S10 can bind to and be a substrate for eithera kinase or a phosphatase. Similarly, S11 cannot be symmetrically converted to S01 or S10. These asymmetries make motif A and motif B independent. Either type of motif can be found in each of the chaos clusters.
In Figures 1C and S2, motif A is highlighted in yellow trimmed with green and motif B is highlighted inyellow trimmed with orange.
To investigate whether the presence of a motif structure is important to chaos behavior in the two-substrate
system, we performed a biased parameter search. Each motif could be placed into four symmetrical positions on a reaction diagram of each Sa and Sb substrate (Figure 3A). We tested a scenario in which each
substrate could have up to two motif As and two motif Bs, respectively. According to the symmetry within
the reaction diagram, 46 patterns of different motif biases cover all possible arrangements of motifs (Figure 3B). Among the 46 motif arrangements, 35 were more likely to show chaos behavior than a system
without any motif bias. These 35 motif arrangements included those that correspond to the clustering analysis results (e.g., cluster 1–7). For example, a structure that has motif A in the positions Sa01/Sa00 and
Sa10/Sa11 corresponds to cluster 1, 2, and 4 and shows a higher probability of having chaos behavior
than in the parameter search without any motif bias (Figure 3B: arrangement #4).
The biased parameter search suggested that motifs with specific arrangements contribute to chaos
behavior (Figure 3B). First, the probability of chaos behavior increases when motif A or motif B are placedin rotationally symmetrical positions (e.g., conversion between S004S01 or between S114S10); arrangements #4, #6, and #7 show higher chaos probabilities as they contain two motif As in rotationallyiScience 24, 101946, January 22, 2021

5

ll

OPEN ACCESS

iScience
ArticleA

B

Figure 3. Functional significance of the number and arrangement of chaos motifs(A) The chaos motif can be subdivided into two types. Motif A comprises either an unphosphorylated substrate (S00) or fully phosphorylated substrate (S11) asa high-affinity substrate that tightly binds to kinase E or phosphatase F, respectively. Motif B comprises a single phosphorylated substrate (S01 or S10) as ahigh-affinity substrate species. In both motifs, the fast reactions that produce high-affinity species are catalyzed by enzymes that are not tightly bound to the
produced high-affinity species. In Figures 2 and S2, motif A is shown in yellow with a green border, while motif B is shown in yellow with an orange border.
(B) The motif structures were imposed in the parameter search by applying biased k and Km values. The positions of imposed motifs are indicated in the
diagrams below each bar chart. A rhombus divided into four subspaces represents the reaction structure shown in (A). For example, the upper left spaceindicates the conversion between S00 and S01. When a green circle is shown in any of the four subspaces, motif A is imposed at the reversible phosphorylation

6

iScience 24, 101946, January 22, 2021

iSciencell

ArticleOPEN ACCESS

Figure 3. Continued
scheme corresponding to the subspace. Similarly, an orange circle denotes motif B imposed at the indicated subspace. Bar charts indicate theprobability of chaotic parameter sets found in the presence of each imposed motif. In the bar charts, motif arrangements correspond to the clusters
found in the random parameter search (Figures 2 and S2), which are indicated as cluster numbers. The red arrow indicates the motif arrangement withthe highest probability of chaos behavior.
See also Figure S3.

symmetrical positions, while arrangements #11, #13, and #14 also exhibit higher chaos probabilities by containing two motif Bs in rotationally symmetrical positions. Second, comparisons between arrangementswith only motif A or motif B show that arrangements with motif A tend to have chaos behavior morethan those with motif B. Third, a certain arrangement of both motif A and motif B substantially increases
the probability of chaos behavior. This is evident by comparison of the probabilities of arrangements #6,#13, and #25; all of these arrangements have motif A or B in the positions Sa004Sa01, Sa114Sa10, andSb004Sb01; however, arrangement #25, which contains two motif As in substrate Sa and one motif B in substrate Sb, has a higher probability of chaos behavior than arrangements #6 and #13. It should also be notedthat even if the number of included motifs increases, the inappropriate motif arrangement decreases chaos
probability to a probability lower than that observed with an unbiased parameter search (e.g., arrangement#38). Among the 46 arrangements, the highest chaos probability was observed for an arrangement of Sawith two motif As in rotationally symmetrical positions and Sb with two motif Bs in rotationally symmetrical
positions (arrangement #46: red arrow in Figure 3B).
With two motif As or two motif Bs in rotationally symmetrical positions, the substrate phosphorylationdephosphorylation reaction tends to flow in a cyclic order. For example, if motif A is placed in the symmetrical positions Sa01/Sa00 and Sa10/Sa11 (e.g., arrangement #4), the fast dephosphorylation reaction leads
to a counter-clockwise phosphorylation-dephosphorylation cycle (Sa01/Sa00/Sa10/Sa11/Sa01 .) (seeFigure 4A). This cyclical bias with stable enzyme-substrate binding is a typical motif arrangement wherelimit cycle oscillation occurs in a single-substrate system (Jolley et al., 2012). Accordingly, two-substrate
systems with two motif As or motif Bs in rotationally symmetrical positions have a higher probability of limit
cycle oscillation, as well as chaos behavior (Figure S3A). Such correlation can also be observed for othermotif arrangements; for example, a motif arrangement with a high probability of chaos behavior tends
to have a high probability of oscillation (Figures S3A and S3B). This relationship is to be expected, giventhat chaos behavior is often caused by the coupling of oscillators. However, a high probability of oscillation
does not fully account for a high probability of chaos behavior because arrangement #46, which has thehighest probability of chaos behavior, does not have the highest probability of oscillation (Figure S3B,
red arrow). Thus, arrangement #46 retains reaction features that are important for chaos behavior.

The chaos motif arrangements create two oscillator structures
To investigate the mechanism underlying the generation of chaos behavior, we considered stereotypicalchaos parameter sets (n = 97) based on motif arrangement #46 (Figures 3B and 4A). For comparison, we
also considered stereotypical oscillation parameter sets (n = 97) based on motif arrangement #46. Both stereotypical parameter sets had symmetrical structures because the setting parameters in the symmetrical
positions (e.g., ka1 and ka8) had identical values.
Bifurcation analyses were performed with the stereotypical parameter sets. Bifurcation maps for stereotypicalchaos sets support the importance of motif structure (Figures 4B and 4C). The motif A structure in substrate
Sa is implemented by high ka4/5 and low Kma2/7 values. Indeed, in the bifurcation maps, chaos behavior isrestricted in the range of the highest ka4/5 and lowest Kma2/7 values. These constrains appear to be important
for avoiding the system’s convergence to a fixed steady state because high ka4/5 and low Kma2/7 values arealso critical for oscillation in the stereotypical oscillation parameter sets (Figure S4). For substrate Sb, motif B
is implemented by high kb1/8 and low Kmb4/6 values. Similar to substrate Sa, chaos behavior is restricted within
the range of high kb1/8 and low Kmb4/6 values, although this constrain appears to be less important than for motif
A in substrate Sa. The importance of a motif B in substrate Sb is characteristic of chaos behavior because theseparameters have a negligible impact in the bifurcation map of the stereotypical oscillation parameter sets (Figure S4). In other words, the appropriate arrangement and parameter space of both motif A in substrate Sa andmotif B in substrate Sb are important for chaos behavior.
Another constrain can be found in ka2. Chaos behavior and oscillation are restricted in the narrow range ofka2/7 values (Figures 4 and S4). Indeed, ka2/7 is critical for period determination in an oscillator composed of

iScience 24, 101946, January 22, 2021

7

ll

OPEN ACCESS

A

B

C

Figure 4. Critical parameters for the oscillation-chaos transition(A) Schematic diagram of the stereotypical chaos/oscillation parameter sets.
(B) Bifurcation analysis of stereotypical chaos parameter sets (n = 97), all of which preserve a motif arrangement shown inthe schematic representation. A heatmap shows the results of bifurcation analysis. White circles indicate the centroid of
the 97 stereotypical parameter sets. The color of the heatmap represents, for a set of stereotyped parameters, theprobability of chaos behavior (red), oscillation (green), or convergence to some steady state (blue) when each indicated8

iScience 24, 101946, January 22, 2021

iScience
ArticleiSciencell

ArticleOPEN ACCESS

Figure 4. Continued
parameter is fixed to an indicated value. For example, if all 97 parameter sets had chaos behavior, the color would be
red, but if half the parameter sets had chaos behavior and the other half had oscillation, the color code would bebrown (50% red +50% green).
(C) The analysis method is the same as that in (B), except that the indicated parameter and its symmetrical pair (e.g., ka1
and ka8, Kmb4 and Kmb5) are simultaneously fixed to an indicated value.
See also Figure S4.

a similar motif A arrangement because this rate constant determines the time required to convert all the
enzyme-trapping statuses (Sa00 and Sa11) to nontrapping statuses (Jolley et al., 2012). In contrast, theperiod-determining parameter for a motif B oscillator is kb3/6. By comparing the parameter ranges of chaos
behavior for ka2/7 and kb3/6 with those of an oscillation, it was found that the faster cycle (higher ka2/7) for the
Sa oscillator relative to that of the Sb oscillator appears to be important for chaos behavior. This expectedperiod difference appears not to be important for oscillation because changes to kb3/6 do not affect thesteady-state behavior of the stereotypical oscillation parameter sets (Figure S4).

Coupling of fast and slow oscillators is raised by a combined motif A and motif B arrangement
To further investigate the underlying mechanism of chaos behavior, we chose one example among the stereotypical chaos parameter sets (Table S1). The system with this example parameter set showed chaos
behavior characterized by a mixture of fast and slow dynamics in the trajectory of phosphorylation statuses(Figure 5A). In Figure 5A, the approximate period of slower dynamics is shown by a blue band, whereas thatof faster dynamics is shown by an orange band (see the Sb plot in Figure 5A). These short and long periods
are sensitive to changes in ka2/7 and kb3/6 values, which are reaction parameters suggested to be critical for
periods of Sa- and Sb-based oscillation (Jolley et al., 2012). A two-fold increase in ka2/7 further shortens theshort period (23 ka2/7: orange band in Figure 5A) while the long period is further lengthened (23 ka2/7: blue
band in Figure 5A). On the other hand, a two-fold decrease in kb3/6 lengthens the long period (0.53 kb3/6:
blue band in Figure 5A) without markedly affecting the short period (0.53 kb3/6: orange band in Figure 5A).
Because the short period is sensitive to ka2/7, we suggest that the dynamics of substrate Sa are primarily
responsible for the fast oscillatory component of chaos behavior. The slow oscillatory component may
be largely determined by the coupling of Sa and Sb because both ka2/7 and kb3/6 affect the longer period.
The different roles of Sa and Sb were further validated by removing one of the Sa or Sb substrates from thechaos system. If the amount of Sb is set to zero with all other parameters maintained, the remaining Sa substrate shows a limit cycle oscillation with a period comparable to the fast oscillatory component of thechaos system (Figure 5B). On the other hand, if the amount of Sa is set to zero, the remaining Sb substrateshows an oscillation with a longer period than that of a Sa-only oscillation (Figure 5C). The different intrinsic
periods of Sa-only and Sb-only oscillations are conserved among the stereotypical parameter sets;
among the 97 stereotypical parameter sets, 85 showed oscillations in the absence of Sa or Sb. In most cases,
a Sa-only oscillation had a higher frequency (i.e., shorter period) than a Sb-only oscillation (Figure 5D). The
unimodal distribution of frequency difference suggests that the coupled Sa and Sb oscillators tend tobehave chaotically when the frequency ratio is not too low or high. These results support the hypothesisthat a Sa oscillator with an intrinsically shorter period contributes to the faster dynamics of chaos behavior
and a Sb oscillator with an intrinsically longer period contributes to the slower dynamics, while the coupling
between Sa and Sb dynamics is required for the slower dynamics of chaos behavior.
In agreement, a mixture of fast and slow dynamics was observed under increased Sa-enzyme affinity (i.e.,decreased Kma2/7) conditions (Figure S5A). Fast dynamics are absent if the affinity of the Sa-enzyme is significantly decreased (Figure S5B), suggesting that Sa-enzyme affinity is important for producing fast dynamics.
In addition, increasing Sb-enzyme affinity disrupts fast dynamics (Figure S5C), while decreasing Sb-enzyme
affinity disrupts slower dynamics (Figure S5D); this suggests that setting Sb-enzyme affinity at an appropriate strength is important for the coexistence of fast and slow dynamics.
The phosphorylation dynamics of Sa and Sb affect each other through the availability of a shared kinase andphosphatase. The availability of free enzymes is restricted by the presence of substrate species with high affinities to these enzymes (i.e., enzyme sequestration). Therefore, the results presented above suggest that enzyme
sequestration by substrate Sa may be important for fast dynamics, whereas slower dynamics are likely related toenzyme sequestration by substrate Sb. Given this possibility, we aimed to directly visualize the behavior of a freekinase, free phosphatase, and the sequestered enzymes using a stochastic simulation (Figure 5E). This stochastic

iScience 24, 101946, January 22, 2021

9

iSciencell

ArticleOPEN ACCESS

A

B

C

D

E

10

iScience 24, 101946, January 22, 2021

iSciencell

ArticleOPEN ACCESS

Figure 5. Dissection of timescales encoded by the chaos motifs(A) Time course of chaotic dynamics for a representative chaos parameter set, the representative parameter set with a
two-fold increase in the parameter values for symmetrical ka2 and ka7, or the representative parameter set with a 0.5-folddecrease in the parameter values for symmetrical kb3 and kb6. The dynamics are mainly composed of oscillatory dynamics
of high and low frequency. A typical period of high-frequency oscillation is shown as an orange stripe, whereas that of alow-frequency oscillation is shown as a blue stripe.
(B) Time course of oscillatory dynamics for the representative parameter set when the amount of Sb is fixed to zero. In theabsence of Sb, the parameter set gives oscillation in Sa with high frequency.
(C) Time course of oscillatory dynamics for the representative parameter set when the amount of Sa is fixed to zero. In theabsence of Sa, the parameter set gives oscillation in Sb with low frequency.
(D) Comparison of frequency differences encoded by oscillatory dynamics intrinsically driven by Sa and Sb. For the 97
stereotypical chaotic parameter sets, the amount of Sa or Sb was fixed to zero. In total, 85 of the stereotypical parameter
sets showed oscillation for the remaining Sb or Sa dynamics. The ratio of Sa-frequency and Sb-frequency indicates that the
Sa substrate is responsible for the high-frequency oscillation in the stereotypical parameter sets.
(E) The upper two plots indicate time evolution of Sa and Sb at each phosphorylation state based on the stochastic
simulation. The lower heatmaps illustrate the fraction of time that each sampled enzyme or substrate molecule(numbered as #1 to #20) spend in an indicated form (e.g., substrate-unbound free E, ES complex from Sb10/Sb11) during a
0.1 s epoch. Sa molecules trap E and F reciprocally at high frequency, while Sb molecules slowly trap and release the E andF, thereby driving low-frequency rhythmicity.
See also Figure S5, Tables S1 and S2.

simulation did not presuppose the quasi-steady-state assumption of the Michaels-Menten reaction scheme.
Instead, the dynamics of the enzyme-substrate complex were simulated explicitly. In the stochastic simulation,
fast and slow dynamics were still observable in Sa and Sb phosphorylation states. Consistent with this, the availability of the free kinase (E) and free phosphatase (F) fluctuated in fast and slow dynamics. This fluctuation in theslow dynamics of E and F was represented as E-depleted and F-depleted phases, which were well matched withthe presence of ESb10/b11 and FSb01/b00 complexes, respectively, both of which are Sb complexes with highaffinities to E and F, respectively. On the other hand, ESa00/a10 and FSa11/a01 complexes, which are Sa complexes with high affinities to E and F, respectively, showed clear oscillations that were well matched with fast dynamics. In summary, the stochastic simulation demonstrated that enzyme sequestration by Sb species represents
the slow dynamics of free enzyme availability, whereas enzyme sequestration by Sa species dominates the fast
dynamics of phosphorylation statuses. The dynamics of chaos behavior may result from the coupling of Sa and Sb
substrates via the availability of enzymes, which varies across timescales (Figure 6).

Discussion
Reversible phosphorylation networks in living cells
To the best of our knowledge, this study is the first computational demonstration of a biochemical systemdescribed only by the Michaelis-Menten reaction scheme exhibiting chaos behavior. Our model can bediscriminated from a set of molecular models employing the synthesis and degradation of components(Decroly and Goldbeter, 1982). Although the coupling of oscillators is understood to be one of the typicalmechanisms underlying chaos generation, the model of reversible phosphorylation described by theMichaelis-Menten scheme helps connect the design principle to the possible molecular properties. We
previously showed that a similar two-site phosphorylation system could be an autonomous oscillator (Jolley
et al., 2012). Mathematically, the model of Jolley et al. (2012) is closely related to the oscillator model of
KaiC protein (Rust et al., 2007), a unique cyanobacteria protein with autokinase and autophosphatase activities. Nonetheless, the model described by the generic Michaelis-Menten scheme allows us to directlycompare parameter (k and Km) combinations, leading to the theory that several kinases, including mammalian proteins, may meet the requirements for composing autonomous oscillators. The work in the present
study may be similarly useful, and below we discuss how a cellular reversible phosphorylation network can
be a chaos generator without the need for eccentric assumptions related to the properties of enzymes andsubstrates.
We employed a model with two substrates, each of which had two modification sites modified by a sharedkinase and phosphatase pair. This setting includes some degree of complexity but should also be
commonly observed in the phosphorylation networks of cells. A previous mass spectrometry analysis revealed that at least 75% of the complete proteome is phosphorylated (Sharma et al., 2014) and thatmany substrates are phosphorylated at multiple residues. On the other hand, the number of kinases andphosphatases is much smaller than the substrate variety (Manning et al., 2002). Thus, different substrates
share a kinase/phosphatase pair for modification of phosphorylation residues. Kinases are also known to

iScience 24, 101946, January 22, 2021

11

ll

OPEN ACCESS

Figure 6. Coupling of posttranslational oscillators encoding different timescales
In a stereotypical arrangement of motif structure, the chaotic dynamics arise from a coupling of two-substrate
modification cycles communicated by the competitive availability of shared enzymes, E and F. One substrate (Sa)intrinsically drives the reciprocal sequestration of the enzymes at high frequency, and the other substrate (Sb) intrinsicallydrives the enzyme sequestration at low frequency. The coupling of two such posttranslational oscillators encoding atdifferent timescales can result in chaotic dynamics without the need for external signals or any typical allosteric regulationof enzymes.

modify a variety of non-motif substrate sequences at a different reaction rates (Miller and Turk, 2018; Ubersax and Ferrell, 2007), although kinases do have preferred substrate sequences, or consensus sequences, atwhich the kinase preferentially phosphorylates target residues. Therefore, a set of shared kinases and phosphatases is naturally observed in the phosphorylation network of cells.
Given this situation, it is reasonable to assume that the number of phosphorylation residues and the number of substrates sharing enzymes is greater than (and therefore more complex than) that in our model.
Nevertheless, our model could be applied to more complex phosphorylation networks in which more

12

iScience 24, 101946, January 22, 2021

iScience
ArticleiSciencell

ArticleOPEN ACCESS

than two phosphorylation sites and substrate species exist. By extending the complexity of multisite phosphorylation in a substrate, it may even be possible to implement the two (distinct) substrates modeled inthis study at distinct protein regions within the same substrate molecule. A similar condition can be found inphosphorylation control of the mammalian circadian protein PERIOD (PER) 2, a relatively large (100 kDa)
protein that is mostly annotated with an intrinsically disordered region (Gustafson and Partch, 2015). A previous study demonstrated that the same kinase, casein kinase I (CKI) d/ε, phosphorylates at multisite phosphorylation clusters located at two distinct protein regions within PER2 (Zhou et al., 2015). Therefore, thephosphorylation of the two regions may compete with each other through different affinities to the CKId/εkinase. Though it is currently unknown whether this circadian kinase-substrate pair exhibits autonomouschaotic phosphorylation, we expect that the phosphorylation network modeled in the present study could
be implemented in a variety of molecules.

Typical chaos motifs found in actual kinases
We found that a typical chaos motif and combination of motifs had chaotic phosphorylation dynamics. Because
the k and Km values of known kinases vary by a few orders of magnitude (Bar-Even et al., 2011) and even the samekinase can show different k and Km values depending on substrate sequences (Ubersax and Ferrell, 2007), it
should be possible to identify enzyme and substrate pairs that satisfy typical chaos motifs.
The chaotic motif shown in Figure 3A requires the formation of a tight complex between the enzyme andthe substrate/product. The arrangement of the motif structure with frequent occurrence of chaos behavior,as shown in Figure 2, suggests that it is important for an enzyme to modify the substrate in a specific order.
These structures are common to the formation of oscillators from reversible phosphorylation (Jolley et al.,2012), and the properties required of the enzyme are also common to the formation of oscillators in several
aspects. For example, CKId/ε is a candidate kinase that could achieve a chaotic motif structure. CKId/εtightly binds to PERs and other substrate proteins (Aryal et al., 2017). Moreover, the substrate/productbinding property is conserved even with short peptide substrates (Shinohara et al., 2017), which plays arole in the temperature-insensitive phosphorylation rate of this kinase. In addition, CKId/ε has a uniqueconsensus sequence of pS-X-X-S* (pS: pre-phosphorylated serine; X: any amino acid; and S*: target serine),
which allows the CKId/ε to sequentially phosphorylate the repeated S-X-X-S sequence from N- to C-terminus. Interestingly, the PER protein substrates have the repeated S-X-X-S sequence, which is known to be
involved in the frequency control of circadian clock oscillation. Recent studies found that a CKId/ε isoform
could phosphorylate serine at the N-terminal proximity of the S-X-X-S repeat and thereby prime forsequential phosphorylation (Fustin et al., 2018; Narasimamurthy et al., 2018). Therefore, CKId/ε andsubstrates with repeated S-X-X-S consensus sequences could be the building blocks for producing posttranslational chaos generators. We note that sequential phosphorylation and substrate/product complex
formation is also found in other kinase/substrate pairs involved in cellular signaling pathways such as cell
cycle control (Miller and Turk, 2018; Valk et al., 2014).

Comparison of the posttranslational oscillator and chaos generator
Our comparisons of chaos behavior and oscillation probabilities among the various combinations of motif
A and B imply that the chaos generator prefers motif A over motif B (see Figure S3B). Compared with motif
A, motif B has a higher tendency to exhibit oscillation rather than chaos behavior. This is consistent with aprevious study (Jolley et al., 2012), which indicated that for single-substrate systems with two reversible
phosphorylation sites, a single-substrate system composed from motif B [‘‘cluster 1’’ oscillator describedby Jolley et al. (2012)] has a higher probability of exhibiting oscillation than that composed from motif A[‘‘cluster 2’’ oscillator described by Jolley et al. (2012)].
Given that typical chaos behavior is generated by the coupled posttranslational oscillators in our system, a question arises: why is a higher chaos behavior probability derived from the presence of motif A rather than from a
two-substrate system composed with motif B that shows the highest probability of oscillation? The importanceof motif A for the generation of chaos behavior can be understood based on the coupling mechanism of the two
substrates, which occurs through enzymatic shearing. Our stochastic simulation indicated that enzymes are trapped and released by Sb species over a slower time course (Figure 5G). To achieve enzyme trapping by Sb, it isnecessary that the enzymes are released from Sa species. The Sa system shown in Figure 5 is composed of motif A(see Figure 4A for the schematic view). From this viewpoint, enzymes are released according to the oscillation
process driven by motif A; focusing on kinase (E) in the Sa species, Sa is mostly trapped by the Sa00 substrate
through high kinase-Sa00 affinity (low Kma2 values). The bound kinase (E) is then released from the Sa00 when

iScience 24, 101946, January 22, 2021

13

ll

OPEN ACCESS

the cyclic phosphorylation/dephosphorylation event proceeds from Sa00 to Sa10. The same scenario can also be
applied to the release of trapped phosphatase (F) by Sa11. Therefore, the trapped enzymes need to be released
through the course of the (de)phosphorylation cycle. This is not the case for an oscillator composed of motif B(i.e., the Sb system in Figure 4A); the kinase (E) trapped by Sb10 is not involved in the cyclic reaction in which Sb10 isdephosphorylated by the phosphatase, nor is the kinase released in the course of the phosphorylation/dephosphorylation cycle. This scenario is validated by the bifurcation analysis shown in Figure 4. The kinase trapping of
Sa can be very tight, as indicated by the bifurcation map in which the chaos area of Kma2 reaches the edge of thesmallest parameter boundary (0.01 mM). Figure S5A also shows that fast and slow dynamics can persist evenwith a very high affinity between Sa00-kinase and Sa11-phosphatase (low Kma2/7 values). In contrast, the kinase
trapping of Sb should be tight but to a lesser extent than that of Sa; the chaos area of Kmb4 is within the smaller
parameter range but cannot reach the smallest boundary (Figure 4). Where Kmb4 values are too low, the enzyme
may not be released from Sb10 efficiently, and the coupling between Sa and Sb might consequently be disrupted.
In summary, although motif B is potent in exhibiting autonomous oscillation, motif A is beneficial to thecoupling status of substrates with other substrate molecules through enzymatic sharing. This suggeststhat by combining different motif structures in different arrangements, it is possible to design biochemical
networks that are suitable for generating not only oscillators that induce constant rhythms but also chaos
generators that elicit spontaneous and complex dynamics in posttranslational modification.

Conclusion
Recent study demonstrated that chaos can arise from relatively simple transcriptional control models (Heltberg et al., 2019), and our study further suggests that autonomous chaos behavior can also arise from ageneric mechanism of reversible multisite phosphorylation. Therefore, multisite phosphorylation networks
could be sources of the spontaneous properties of cellular function, such as the stochastic nature of cellulardifferentiation (Furusawa and Kaneko, 2001) and the spontaneous and excitable activation of intracellularsignaling observed in unicellular organisms (Devreotes et al., 2017). Thermal fluctuations in molecular activity, as well as transduction and amplification, have previously been modeled as a source of spontaneous
activation in intracellular biochemical reactions (Oosawa, 2001, 2007; Ooyama and Shibata, 2011; Shibata
and Ueda, 2008; Ueda and Shibata, 2007). The present study suggests that spontaneous and chaotic dynamics may also be generated from the ubiquitous setting of reversible posttranslational networks. As asimilar line of research, it will also be interesting to investigate how the enzyme-substrate network avoids
chaos behavior if spontaneous action leads to uncontrollable outcomes (Gelens et al., 2015).

Limitations of the study
In this study, we did not conduct an analysis of chaos system based on the analytical solution of mathematical models. Instead, we used a set of analysis based on random parameter searching, numerical integration of ordinary differential equations, and clustering analysis of chaos parameter sets. These analyses have
limitations in the interpretation of the obtained results. First, although we could not find chaos in the systemwith single substrate having two phosphorylation sites, we cannot exclude the possibility that chaos can beobserved with certain parameter combinations in the system of single substrate. Second, the mechanism ofchaos generation explained in this study (e.g. Figure 6) is based on the stereotypical parameter combination illustrated as in Figure 4, and it is highly likely that there are certain parameter sets yielding chaos dynamics by the mechanism that are difficult to be categorized as ‘‘two-coupled fast and slow oscillators’’.
Finally, since the random parameter search was performed against reversible phosphorylation models
with the Michaelis-Menten approximation, the results of the random parameter search may show different
trends in the distribution of the chaotic parameters if similar reversible phosphorylation systems aremodeled without using the Michaelis-Menten approximation.

Resource availability
Lead contact
Further information and requests for resources should be directed to and will be fulfilled by the Lead
Contact, Hiroki R. Ueda (uedah-tky@umin.ac.jp)

Materials availability
This study did not generate new unique reagents.

14

iScience 24, 101946, January 22, 2021

iScience
ArticleiSciencell

ArticleOPEN ACCESS

Data and code availability
The code generated during this study are available at GitHub https://github.com/DSP-sleep/
PT_Chaotic_Oscillators.

Methods
All methods can be found in the accompanying Transparent Methods supplemental file.

Supplemental information
Supplemental Information can be found online at https://doi.org/10.1016/j.isci.2020.101946.

Acknowledgments
This work was supported by grants from the Brain/MINDS JP21dm0207049, Science and Technology
Directorate Platform Program for Advanced Biological Medicine JP21am0401011, AMED-CREST
JP21gm0610006 (AMED/MEXT) (H.R.U.), Grant-in-Aid for Scientific Research (S) JP25221004 (JSPS
KAKENHI) (H.R.U.) and Scientific Research (C) JP20K06576 (JSPS KAKENHI) (K.L.O.), HFSP Research Grant
Program RGP0019/2018 (HFSP) (H.R.U.), ERATO JPMJER2001 (JST) (H.R.U.) and an intramural Grant-in-Aid
from the RIKEN BDR (H.R.U.). The authors would like to thank Enago (www.enago.jp) for the English language review.

Author contributions
H.Q.Y., K.L.O., and H.R.U. designed the study. H.Q.Y. performed most of the simulation study. H.Q.Y.,
K.L.O., and H.R.U. analyzed the mechanism of chaos generation and wrote the manuscript.

Declaration of interests
H.R.U. is a member of the iScience’s editorial advisory board.

Received: October 5, 2020
Revised: December 2, 2020
Accepted: December 11, 2020
Published: January 22, 2021
References
Aihara, K., Matsumoto, G., and Ichikawa, M.
(1985). An alternating periodic-chaotic sequenceobserved in neural oscillators. Phys. Lett. A 111,
251–255.

Conradi, C., Mincheva, M., and Shiu, A. (2019).
Emergence of oscillations in a mixed-mechanism
phosphorylation system. Bull. Math. Biol. 81,
1829–1852.

Aryal, R.P., Kwak, P.B., Tamayo, A.G., Gebert, M.,
Chiu, P.L., Walz, T., and Weitz, C.J. (2017).
Macromolecular assemblies of the mammalian
circadian clock. Mol. Cell 67, 770–782 e776.

Conradi, C., and Shiu, A. (2018). Dynamics of
posttranslational modification systems: recent
progress and future directions. Biophysical J. 114,
507–515.

Bar-Even, A., Noor, E., Savir, Y., Liebermeister,W., Davidi, D., Tawfik, D.S., and Milo, R. (2011).
The moderately efficient enzyme: evolutionary
and physicochemical trends shaping enzyme
parameters. Biochemistry 50, 4402–4410.

Decroly, O., and Goldbeter, A. (1982).
Birhythmicity, chaos, and other patterns of
temporal self-organization in a multiply regulated
biochemical system. Proc. Natl. Acad. Sci. U S A
79, 6917–6921.

Bieler, J., Cannavo, R., Gustafson, K., Gobet, C.,
Gatfield, D., and Naef, F. (2014). Robust
synchronization of coupled circadian and cell
cycle oscillators in single mammalian cells. Mol.
Syst. Biol. 10, 739.

Devreotes, P.N., Bhattacharya, S., Edwards, M.,
Iglesias, P.A., Lampert, T., and Miao, Y. (2017).
Excitable signal transduction networks indirected cell migration. Annu. Rev. Cell Dev. Biol.
33, 103–125.

Chen, Z., Kibler, R.D., Hunt, A., Busch, F., Pearl, J.,
Jia, M., VanAernum, Z.L., Wicky, B.I.M., Dods, G.,
Liao, H., et al. (2020). De novo design of protein
logic gates. Science 368, 78–84.

Droin, C., Paquet, E.R., and Naef, F. (2019). Lowdimensional dynamics of two coupled biologicaloscillators. Nat. Phys. 15, 1086–1094.

Chickarmane, V., Kholodenko, B.N., and Sauro,
H.M. (2007). Oscillatory dynamics arising fromcompetitive inhibition and multisite
phosphorylation. J. Theor. Biol. 244, 68–76.

Durston, A.J. (1974). Pacemaker mutants of
Dictyostelium discoideum. Dev. Biol. 38, 308–319.
Fink, T., Lonzaric, J., Praznik, A., Plaper, T.,
Merljak, E., Leben, K., Jerala, N., Lebar, T.,

Strmsek, Z., Lapenta, F., et al. (2019). Design offast proteolysis-based signaling and logic circuits
in mammalian cells. Nat. Chem. Biol. 15, 115–122.
Furusawa, C., and Kaneko, K. (2001). Theory ofrobustness of irreversible differentiation in a stem
cell system: chaos hypothesis. J. Theor. Biol. 209,
395–416.
Fustin, J.M., Kojima, R., Itoh, K., Chang, H.Y., Ye,
S., Zhuang, B., Oji, A., Gibo, S., Narasimamurthy,
R., Virshup, D., et al. (2018). Two Ck1deltatranscripts regulated by m6A methylation codefor two antagonistic kinases in the control of thecircadian clock. Proc. Natl. Acad. Sci. U S A 115,
5980–5985.
Gao, X.J., Chong, L.S., Kim, M.S., and Elowitz,
M.B. (2018). Programmable protein circuits inliving cells. Science 361, 1252–1258.
Gelens, L., Huang, K.C., and Ferrell, J.E., Jr.
(2015). How does the Xenopus laevis embryoniccell cycle avoid spatial chaos? Cell Rep. 12,
892–900.
Gerard, C., and Goldbeter, A. (2012). Entrainmentof the mammalian cell cycle by the circadianclock: modeling two coupled cellular rhythms.
PLoS Comput. Biol. 8, e1002516.

iScience 24, 101946, January 22, 2021

15

iSciencell

ArticleOPEN ACCESS

Goldbeter, A. (2002). Computational approaches
to cellular rhythms. Nature 420, 238–245.
Gustafson, C.L., and Partch, C.L. (2015). Emerging
models for the molecular basis of mammalian
circadian timing. Biochemistry 54, 134–149.
Heltberg, M.L., Krishna, S., and Jensen, M.H.
(2019). On chaotic dynamics in transcription
factors and the associated effects in differential
gene regulation. Nat. Commun. 10, 71.
Izhikevich, E.M. (2007). Dynamical Systems in
Neuroscience (MIT Press).
Jolley, C.C., Ode, K.L., and Ueda, H.R. (2012). A
design principle for a posttranslational
biochemical oscillator. Cell Rep. 2, 938–950.
Kauffman, S.A. (1992). Origins of order in
evolution: self-organization and Selection. In
Understanding Origins: Contemporary Views on
the Origin of Life, Mind and Society, F.J. Varela
and J.-P. Dupuy, eds. (Springer Netherlands),
pp. 153–181.
Kholodenko, B.N. (2006). Cell-signalling
dynamics in time and space. Nat. Rev. Mol. Cell
Biol. 7, 165–176.
Langton, C.G. (1990). Computation at the edge ofchaos: phase transitions and emergent
computation. Physica D Nonlinear Phenomena
42, 12–37.
Liu, P., Kevrekidis, I.G., and Shvartsman, S.Y.
(2011). Substrate-dependent control of ERK
phosphorylation can lead to oscillations.
Biophysical J. 101, 2572–2581.
Manning, G., Whyte, D.B., Martinez, R., Hunter,T., and Sudarsanam, S. (2002). The protein kinase
complement of the human genome. Science 298,
1912–1934.
Markevich, N.I., Hoek, J.B., and Kholodenko, B.N.
(2004). Signaling switches and bistability arisingfrom multisite phosphorylation in protein kinasecascades. J. Cell Biol. 164, 353–359.
Markevich, N.I., Tsyganov, M.A., Hoek, J.B., andKholodenko, B.N. (2006). Long-range signalingby phosphoprotein waves arising from bistability
in protein kinase cascades. Mol. Syst. Biol. 2, 61.
Martiel, J.L., and Goldbeter, A. (1985).
Autonomous chaotic behaviour of the slimemould Dictyostelium discoideum predicted by amodel for cyclic AMP signalling. Nature 313,
590–592.
Matsuo, T., Yamaguchi, S., Mitsui, S., Emi, A.,
Shimoda, F., and Okamura, H. (2003). Control
mechanism of the circadian clock for timing of cell
division in vivo. Science 302, 255–259.
Miller, C.J., and Turk, B.E. (2018). Homing in:mechanisms of substrate targeting by protein
kinases. Trends Biochem. Sci. 43, 380–394.
Murray, A.W., and Kirschner, M.W. (1989). Cyclin
synthesis drives the early embryonic cell cycle.
Nature 339, 275–280.

16

iScience 24, 101946, January 22, 2021

Nagoshi, E., Saini, C., Bauer, C., Laroche, T., Naef,F., and Schibler, U. (2004). Circadian gene
expression in individual fibroblasts: cellautonomous and self-sustained oscillators passtime to daughter cells. Cell 119, 693–705.
Narasimamurthy, R., Hunt, S.R., Lu, Y., Fustin,
J.M., Okamura, H., Partch, C.L., Forger, D.B., Kim,
J.K., and Virshup, D.M. (2018). CK1delta/epsilon
protein kinase primes the PER2 circadianphosphoswitch. Proc. Natl. Acad. Sci. U S A 115,
5986–5991.
Olsen, L.F., and Degn, H. (1977). Chaos in an
enzyme reaction. Nature 267, 177–178.
Oosawa, F. (2001). Spontaneous signal
generation in living cells. Bull Math. Biol. 63,
643–654.
Oosawa, F. (2007). Spontaneous activity of living
cells. Bio Syst. 88, 191–201.
Ooyama, S., and Shibata, T. (2011). Hierarchicalorganization of noise generates spontaneous
signal in Paramecium cell. J. Theor. Biol. 283, 1–9.
Palmeirim, I., Henrique, D., Ish-Horowicz, D., andPourquie, O. (1997). Avian hairy gene expression
identifies a molecular clock linked to vertebrate
segmentation and somitogenesis. Cell 91,
639–648.
Pomerening, J.R., Kim, S.Y., and Ferrell, J.E., Jr.
(2005). Systems-level dissection of the cell-cycle
oscillator: bypassing positive feedback produces
damped oscillations. Cell 122, 565–578.
Pomerening, J.R., Sontag, E.D., and Ferrell, J.E.,
Jr. (2003). Building a cell cycle oscillator:
hysteresis and bistability in the activation of Cdc2.
Nat. Cell Biol. 5, 346–351.
Qiao, L., Nachbar, R.B., Kevrekidis, I.G., andShvartsman, S.Y. (2007). Bistability andoscillations in the Huang-Ferrell model of MAPK
signaling. PLoS Comput. Biol. 3, 1819–1826.
Raskin, D.M., and de Boer, P.A. (1999). Rapid
pole-to-pole oscillation of a protein required for
directing division to the middle of Escherichia
coli. Proc. Natl. Acad. Sci. U S A 96, 4971–4976.
Rasmussen, R., Jensen, M.H., and Heltberg, M.L.
(2017). Chaotic dynamics mediate Brain statetransitions, driven by changes in extracellular ion
concentrations. Cell Syst. 5, 591–603 e594.
Rust, M.J., Markson, J.S., Lane, W.S., Fisher, D.S.,
and O’Shea, E.K. (2007). Orderedphosphorylation governs oscillation of a threeprotein circadian clock. Science 318, 809–812.
Shankaran, H., Ippolito, D.L., Chrisler, W.B.,
Resat, H., Bollinger, N., Opresko, L.K., and Wiley,
H.S. (2009). Rapid and sustained nuclearcytoplasmic ERK oscillations induced by
epidermal growth factor. Mol. Syst. Biol. 5, 332.
Sharma, K., D’Souza, R.C., Tyanova, S., Schaab,C., Wisniewski, J.R., Cox, J., and Mann, M. (2014).
Ultradeep human phosphoproteome reveals a
distinct regulatory nature of Tyr and Ser/Thrbased signaling. Cell Rep. 8, 1583–1594.

Shibata, T., and Ueda, M. (2008). Noise
generation, amplification and propagation inchemotactic signaling systems of living cells. Bio
Syst. 93, 126–132.
Shinohara, Y., Koyama, Y.M., Ukai-Tadenuma, M.,
Hirokawa, T., Kikuchi, M., Yamada, R.G., Ukai, H.,
Fujishima, H., Umehara, T., Tainaka, K., et al.
(2017). Temperature-sensitive substrate andproduct binding underlie temperaturecompensated phosphorylation in the clock. Mol.
Cell 67, 783–798 e720.
Skarda, C.A., and Freeman, W.J. (1987). How
brains make chaos in order to make sense of theworld. Behav. Brain Sci. 10, 161–173.
Strogatz, S.H. (1994). Nonlinear Dynamics andChaos: With Applications to Physics, Biology,
Chemistry, and Engineering (Reading, Mass
(Addison-Wesley Pub.).
Sugai, S.S., Ode, K.L., and Ueda, H.R. (2017). A
design principle for an autonomous posttranslational pattern formation. Cell Rep. 19,
863–874.
Suwanmajo, T., and Krishnan, J. (2015). Mixedmechanisms of multi-site phosphorylation. J. R.
Soc. Interface 12, 20141405.
Thomson, M., and Gunawardena, J. (2009).
Unlimited multistability in multisite
phosphorylation systems. Nature 460, 274–277.
Ubersax, J.A., and Ferrell, J.E., Jr. (2007).
Mechanisms of specificity in protein
phosphorylation. Nat. Rev. Mol. Cell Biol. 8,
530–541.
Ueda, M., and Shibata, T. (2007). Stochastic signal
processing and transduction in chemotactic
response of eukaryotic cells. Biophysical J. 93,
11–20.
Valk, E., Venta, R., Ord, M., Faustova, I.,
Koivomagi, M., and Loog, M. (2014). Multistep
phosphorylation systems: tunable components ofbiological signaling circuits. Mol. Biol. Cell 25,
3456–3460.
Yan, J., and Goldbeter, A. (2019a). Multirhythmicity generated by coupling two cellularrhythms. J. R. Soc. Interface 16, 20180835.
Yan, J., and Goldbeter, A. (2019b). Robust
synchronization of the cell cycle and the circadianclock through bidirectional coupling. J. R. Soc.
Interface 16, 20190376.
Yang, Q., Pando, B.F., Dong, G., Golden, S.S.,
and van Oudenaarden, A. (2010). Circadiangating of the cell cycle revealed in single
cyanobacterial cells. Science 327, 1522–1526.
Yoshioka-Kobayashi, K., Matsumiya, M., Niino, Y.,
Isomura, A., Kori, H., Miyawaki, A., andKageyama, R. (2020). Coupling delay controls
synchronized oscillation in the segmentation
clock. Nature 580, 119–123.
Zhou, M., Kim, J.K., Eng, G.W., Forger, D.B., andVirshup, D.M. (2015). A Period2 phosphoswitch
regulates and temperature compensatescircadian period. Mol. Cell 60, 77–88.

iScience, Volume 24

Supplemental Information

A design principle for posttranslational
chaotic oscillators
Hiroto Q. Yamaguchi, Koji L. Ode, and Hiroki R. Ueda

Supplemental Figures

A

Reversible phosphorylationof signle substrate with
one phosphorylation site

Kinase and phosphatase
reactions

Reversible phosphorylationof signle substrate withtwo phosphorylation sites
S01
5

5
m

1

,k

K

,k

7

7

3

P
8

8

P

4

2

P

S10

m

,k

6

4

,k

K

,k

m

6

2

S0

S11

K

S00

,k

K

m

1

,k

P

Km

Km2, k2

Km2, k2

S1

Km

S1

Km1, k1

S0

F
P

Km

P

3

S1

Km

Km1, k1

S0

,k

E

P

B

Concentration (μM）

1000

500
0
1000

500
0
4880
Sa00

4940
Time (min)
Sa10
Sa01

Sa11

Sb00

Sb01

Sb11

ratio

n of

S

a00

Trajectory

25

(µM

)

e
nc

on

i
at

r

nt

0

cent

20

0
Con

of

S

a10

)

M

(µ

1
a0

C

o

State at t = 5000
started from slightly different initial states

Power (μM2*min)

10

00

50

Concentration of S

150

Power (μM2*min)

D
(µM)

C

Sb10

6
4

3

5000

× 107

2

2 1
0 0
6 3
4 2
2 1
0 0
0.00

0.02

0.04

0.06 0.08 0.10
Frequency (min-1)

0.12

Sa00

Sa01

Sa10

Sa11

Sb00

Sb01

Sb10

Sb11

Figure S1. Yamaguchi et al.

0.14

Figure S1. Oscillatory behavior of a coupled dual-phosphorylation system(A) Scheme of reversible phosphorylation models.
(B) An example time course of oscillatory phosphorylation dynamics.
(C) The example oscillatory dynamics of t = 0 min to t = 5,000 min are projected to S a00-Sa01-Sa10
space (green line). Red dots indicate the phosphorylation states at t = 5,000 min obtained from 200
runs of the simulation with slightly different initial conditions.
(D) Power spectrum of the example oscillatory trajectory analyzed by FFT.
See also Figure 1.

Cluster 2
Sa01
F
Sa11

Sb00

Sa10

Sb10

Sa

Sb

0.4

E
Sa00

Sb11

F
Sa11

Sa01

Sb00

Sa10

Sb10

Sa

Sb
0.4

0.3

0.3

0.3

0.2

0.2

0.2

0.1

0.1

0.1

0.0

0

1

2

3 0
Log10 k (min-1)

1

2

3

0.0

0

1

2

3 0
1
Log10 k (min-1)

2

3

0.0

0.4

0.4

0.4

0.3

0.3

0.3

0.2

0.2

0.2

0.1

0.1

0.1

−2 −1

0

1

2 3 −2 −1
Log10 Km (μM)

0

1

2

3

0.0
−2 −1

0

1

2 3 −2 −1
Log10 Km (μM)

Cluster 5
F
Sa01

Sb01
Sb00

Sa11

Sa00

Sb11

Sa11

E
Sb10

Sa10

Sa

Sb

Sa

3

Sb00

E

0.2
0.1
1

2

3

0
Log10 k (min-1)

1

2

3

0.2

0.1

0.1

0.0

−2 −1

0

1

2 3 −2 −1
Log10 Km (μM)

0

1

2

3

2

0

1

3 0
1
Log10 k (min-1)

2

3 −2 −1
Log10 Km (μM)

2

0

1

3

2

3

1

2

3 0
1
Log10 k (min-1)

2

3

0

Sa11

Sb00

Sb11
E

Sa10

Sb10

Sa

Sb

3

k1
k2
k3
k4
k5
k6
k7
k8

1

2

3 0
1
Log10 k (min-1)

2

3

Km1
Km2
Km3
Km4
Km5
Km6
Km7
Km8

0.3
0.2

0.1

0.1

0.0
−2 −1

0.0
0

1

2

3 −2 −1
Log10 Km (μM)

0

1

2

3

Figure S2. Yamaguchi et al.

−2 −1

0

Km1
Km2
Km3
Km4
Km5
Km6
Km7
Km8

E
Sb01

0.4

0.2

k1
k2
k3
k4
k5
k6
k7
k8

Cluster 7

0.0
0

0.3

0.1

1

0.3

0.4

0.2

Sb

0.4

0.4
0.3

Sa

Sb

0.2

0.0

Sb11
Sb10

−2 −1

Sb10

Sb00

Sa10

Sa00

Sb11

0.3

0

0.0

F
Sa11

F
Sa01

0.4

0.3

0.0

2

Sb01

Sa10

0.4

1

0

Cluster 6

E F
Sa01

Sa00

0

Sb01

E
Sa00

Sb11

0.4

0.0

Probability

Cluster 4
F
Sb01

Sa01

Sb01

E
Sa00

Percentage

Cluster 312

3 −2 −1 0
Log10 Km (μM)

12Figure S2. Parameter motifs found in chaotic parameter sets
Schematic representation and parameter histograms of the chaotic parameter clusters shown in thesame manner as in Figure 2.
See also Figure 2.

A

Probability of oscillation (x 10-3)

Motif A only

Motif B only

1x Motif A, 1x motif B

10

7.5
Cluster
1, 2, 4

5

Motif

0
Arrangement #Cluster6, 7

Cluster
5

2.5

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

Sa

Sb

Probability of oscillation (x 10-3)

2x Motif A, 1x motif B

1x Motif A, 2x motif B

2x Motif A, 2x motif B

10

7.5

Cluster
3

5

2.5
0
21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

Sa

Sb

B

0.012

0.01
Motif A only
Motif B only
1x Motif A, 1x motif B

0.008
Probability of oscillation

Motif

20

2x Motif A, 1x motif B
1x Motif A, 2x motif B
2x Motif A, 2x motif B

0.006

0.004

0.002

0

0

0.000002

0.000004

0.000006

0.000008

Probability of chaos

Figure S3. Yamaguchi et al.

0.00001

45

46

Figure S3. Specific arrangements of motif structure are important for chaos generation(A) Bar charts indicate the probability of oscillatory parameter sets found in the presence of each
imposed motif. Imposed motif arrangements are shown in the same manner as in Figure 3.
(B) A comparison of the probability of chaos behavior and oscillation for parameters in the presenceof each combination of imposed motifs. The red arrow indicates the motif arrangement with thehighest probability of chaos behavior (i.e., arrangement #46).
See also Figure 3.

A

Asymmetric purturbation
Sa parameters
3

ka1

ka2

ka3

ka4

Sb parameters

Kma1 Kma2 Kma3 Kma4 kb1

kb2

kb3

kb4

Kmb1 Kmb2 Kmb3 Kmb4

3

2
1
0

1

Log10 Km (μM)

Log10 k (min-1)

2

−1
−2

0
Chaos

B

Oscillation

Convergence

Symmetric purturbation
Sa parameters3

Sb parameters

ka1/8 ka2/7 ka3/6 ka4/5 Kma1/8 Kma2/7 Kma3/6 Kma4/5 kb1/8 kb2/7 kb3/6 kb4/5 Kmb1/8 Kmb2/7 Kmb3/6 Kmb4/5

3

2
1
0

1

−1
−2

0
Chaos

Oscillation

Convergence

Figure S4. Yamaguchi et al.

Log10 Km (μM)

Log10 k (min-1)

2

Figure S4. Bifurcation analysis of stereotypical oscillation parameter sets.
(A) Bifurcation analysis was conducted for stereotypical oscillatory parameter sets (n = 97), all ofwhich preserve the motif arrangement shown in the schematic representation in (A). The method
used is same as that reported in Figure 4B.
(B) The same analysis as that in (A) is shown, except that the indicated parameter and its
symmetrical pair (e.g., ka1 and ka8, Kmb4 and Kmb5) are simultaneously fixed to an indicated value.
See also Figure 4.

0
1

0

4.96

Concentration (μM)

0
1

0

5.00

Sb

4.96

B
1

1

Sa

Time (x1000 min)

1

1

5.00

Sb

Time (x1000 min)

S00

Sa

0
1

0

5.00

S01

0.01 x Kmb4/5

4.96

5.00

Sb

4.96

D

100 x Kma2/7

4.96

4.96

0

0

5.00

Concentration (μM)

Concentration (μM)

1

Sa

C

0.01 x Kma2/7

Concentration (μM)

A

Sa

Time (x1000 min)

5.00

100 x Kmb4/5

4.96

5.00

Sb

4.96

Time (x1000 min)

S10

Figure S5. Yamaguchi et al.

S11

5.00

Figure S5. Roles of enzyme trap efficiency in the coupling of two oscillators(A, B) Time course of oscillatory dynamics for the representative parameter set with a 0.01-folddecrease (A) or 100-fold increase (B) in the parameter values for symmetrical Kma2 and Kma7.
Increasing the enzyme sequestration efficiency by Sa substrate (A) did not severely affect chaotic
dynamics; reducing the enzyme sequestration by Sa (B) only allowed Sb-sequestration to drive therhythmic synchronization of phosphorylation states such that the entire dynamics became
dominated by slow-frequency oscillation.
(C, D) Time course of oscillatory dynamics for the representative parameter set with a 0.01-folddecrease (C) or 100-fold increase (D) in the parameter values for symmetrical Kmb4 and Kmb5.
Increasing the enzyme sequestration efficiency by S b substrate (C) masked the high-frequency
rhythmicity driven by Sa; reducing the enzyme sequestration by Sb (D) only allowed Sasequestration to drive the rhythmic synchronization of phosphorylation states such that the entire
dynamics became dominated by high-frequency oscillation.

Supplemental Tables
Table S1. Example stereotypical chaos parameter set, related to Figure 5.
Parameter

Example stereotypical chaos parameter setk (min-1), Km (μM)
ka1

47.24854606

ka2

103.9531911

ka3

10.08171512

ka4

1060.571158

ka5

1060.571158

ka6

10.08171512

ka7

103.9531911

ka8

47.24854606

Kma1

1.471248283

Kma2

0.008378714

Kma3

57.41102352

Kma4

2.597629575

Kma5

2.597629575

Kma6

57.41102352

Kma7

0.008378714

Kma8

1.471248283

kb1

209.3650037

kb2

10.6326728

kb3

19.66831864

kb4

12.42448634

kb5

12.42448634

kb6

19.66831864

kb7

10.6326728

kb8

209.3650037

Kmb1

87.77317768

Kmb2

38.41198729

Kmb3

4.314899945

Kmb4

0.04746365

Kmb5

0.04746365

Kmb6

4.314899945

Kmb7

38.41198729

Kmb8

87.77317768

Table S2. Example stereotypical chaos parameter set for stochastic simulation, related to Figure 5.
Parameterk (min-1), k

Example stereotypical chaos parameter setc(min-1 μM-1), kuc(min-1)(stochastic simulation)

ka1

47.24854606

ka2

103.9531911

ka3

10.08171512

ka4

1060.571158

ka5

1060.571158

ka6

10.08171512

ka7

103.9531911

ka8

47.24854606

kuca1 - kuca8

10

kca1

32.794

kca2

12526

kca3

0.19302kca4

408.67

kca5

408.67kca60.19302

kca7

12526

kca8

32.794kb1

209.3650037

kb2

10.6326728

kb3

19.66831864

kb4

12.42448634

kb5

12.42448634

kb6

19.66831864

kb7

10.6326728

kb8

209.3650037

kucb1 - kucb8

10

kcb1

2.3967kcb2

0.30284

kcb3

4.79

kcb4

282.84kcb5

282.84

kcb6

4.79

kcb7

0.30284kcb8

2.3967

Transparent Methods
Modeling of coupled-posttranslational oscillator system
The coupled-posttranslational oscillator system was formulated as a set of ten ordinary differential
equations that describe the temporal evolution of the concentrations of eight substrate
phosphorylation states and two unbound enzymes.
d[Sa00 ]
𝑘𝑎1
𝑘a2
𝑘a5
𝑘a6
[Sa01 ] +[S ]) [F]= −(+) [E][Sa00 ] + (d𝑡
𝐾𝑚a1 𝐾𝑚a2
𝐾𝑚a5
𝐾𝑚a6 a10(1)d[Sa01 ]
𝑘a1
𝑘a3
𝑘a5
𝑘a7
[Sa00 ] −[Sa01 ]) [E] + (−[Sa01 ] +[S ]) [F]
=(
d𝑡
𝐾𝑚a1
𝐾𝑚a3
𝐾𝑚a5
𝐾𝑚a7 a11(2)

d[Sa10 ]
𝑘a2
𝑘a4
𝑘a6
𝑘a8
[S ] −[S ]) [E] + (−[S ] +[S ]) [F]=(d𝑡𝐾𝑚a2 a00
𝐾𝑚a4 a10
𝐾𝑚a6 a10
𝐾𝑚a8 a11(3)

d[Sa11 ]
𝑘a3
𝑘a4
𝑘a6
𝑘a8
[S ] +[S ]) [E] − (=(+) [F][Sa11 ]d𝑡
𝐾𝑚a3 a01
𝐾𝑚a4 a10
𝐾𝑚a6 𝐾𝑚a8(4)d[Sb00 ]
𝑘b1
𝑘b2
𝑘b5
𝑘b6
[Sb01 ] +[S ]) [F]= −(+) [E][Sb00 ] + (d𝑡
𝐾𝑚b1 𝐾𝑚b2
𝐾𝑚b5
𝐾𝑚b6 b10(5)

d[Sb01 ]
𝑘b1
𝑘b3
𝑘b5
𝑘b7
[Sb00 ] −[Sb01 ]) [E] + (−[Sb01 ] +[S ]) [F]
=(
d𝑡
𝐾𝑚b1
𝐾𝑚b3
𝐾𝑚b5
𝐾𝑚b7 b11(6)

d[Sb10 ]
𝑘b2
𝑘b4
𝑘b6
𝑘b8
[Sb00 ] −[Sb10 ]) [E] + (−[Sb10 ] +
[S ])
=(
d𝑡
𝐾𝑚b2
𝐾𝑚b4
𝐾m𝑏6
𝐾𝑚b8 b11(7)

d[Sb11 ]
𝑘b3
𝑘b4
𝑘b6
𝑘b8
[S ] +[S ]) [E] − (=(+) [F][Sb11 ]d𝑡
𝐾𝑚b3 b01
𝐾𝑚b4 b10
𝐾mb6 𝐾𝑚b8

(8)

[E] =

[F] =

𝐸total[S ] [S ] [S ] [S ] [S ] [S ] [S ] [S ]1 + a00 + a00 + a01 + a10 + b00 + b01 + b10 + b10
𝐾𝑚a1
𝐾𝑚a2
𝐾𝑚a3
𝐾𝑚a4
𝐾𝑚b1
𝐾𝑚b2
𝐾𝑚b3
𝐾𝑚b4
𝐹total[Sa01 ] [Sa10 ] [Sa11 ] [Sa11 ] [Sb01 ] [Sb10 ] [Sb11 ] [Sb11 ]
1+
+
+
+
+
+
+
+
𝐾𝑚a5
𝐾𝑚a6
𝐾𝑚a7
𝐾𝑚a8
𝐾𝑚b5
𝐾𝑚b6
𝐾𝑚b7
𝐾𝑚b8

(9)

(10)

In these equations, [Sa00], [Sa01], [Sa10], and [Sa11] represent the concentrations of the four
phosphorylation states of substrate A; [Sb00], [Sb01], [Sb10], and [Sb11] represent the concentrations of
the four phosphorylation states of substrate B; [E] and [F] are the concentrations of free (unbound)kinase and phosphatase; ka1-ka8 and kb1-kb8 are the reaction rate constants; Kma1-Kma8 and Kmb1-Kmb8
are the Michaelis−Menten constants; and Etotal and Ftotal are the total concentrations of kinase andphosphatase.

Random parameter search for chaotic and oscillatory parameter sets
We randomly generated parameter sets that consisted of 32 constants (16 reaction rate constants
and 16 Michaelis−Menten constants) and then numerically solved the equations to find chaotic
parameter sets. Reaction rate constants ka1-ka8 and kb1-kb8 were independently generated fromexponential distributions bounded between 1 and 1000 min−1. Michaelis−Menten constants Kma1Kma8 and Kmb1-Kmb8 were independently generated from exponential distributions bounded between

0.01 and 1000 μM. The integration began from a state in which all the substrate molecules were
dephosphorylated and all the enzyme molecules were unbounded ([Sa00] = [Sb00] = 1000 μM; [E] =[F] = 20 μM). The solution for each parameter set could be convergent, oscillatory, or chaotic. We
judged a solution to be chaotic when its power spectrum was continuous rather than discrete; thespectrum was firstly numerically classified as not discrete if the proportion of the power at the peakfrequency to the total power was below 15%. Then we plotted those candidate continuous
spectrums and visually checked their continuity. We further checked the chaotic parameter sets toensure that, after sufficient time had passed (t = 5000 min), the solution was divergent because ofsmall changes in the initial state of the integration. The degree of divergence is visually inspected.
Oscillatory parameter sets were collected in the same manner as chaotic parameter sets. Solutions
were considered to be oscillatory when the power spectrum was discrete.

Clustering of chaotic parameter sets
Chaotic parameter sets collected in the random parameter search were standardized and clusteredby performing a clustering analysis. In the standardization procedure, each constant wastransformed so that the defined lower bound, upper bound, and the midpoint between them mappedto −1, 1, and 0, respectively. To achieve this, reaction rate constants were log-transformed withbase 10, 1.5 was subtracted from the value [= (0 + 3)/2], and then this value was divided by 1.5 [=(3 − 0)/2]. Michaelis−Menten constants were log-transformed with base 10, 0.5 was subtracted fromthe value [= (3 + (−2))/2], and then this value was divided by 2.5 [= (3 − (−2))/2]. In the clustering
analysis, standardized parameters were hierarchically clustered using Ward’s algorithm. We used a
special distance metric for the algorithm to consider the three symmetries of the system. Thesesymmetries existed (1) in the phosphorylation states of each substrate (i.e., exchanging the positionof Sa01 and Sa10, or exchanging the position of Sb01 and Sb10, conserves the overall structure), (2) asan enzymatic symmetry (i.e., exchanging the role of kinase and phosphatase conserves the overall
structure), and (3) as a substrate symmetry (i.e., exchanging Sa and Sb conserves the overall
structure). By considering the first symmetry, the parameter set that exchanged parameter {ka1, Kma1,
ka2, Kma2, ka3, Kma3, ka4, Kma4, ka5, Kma5, ka6, Kma6, ka7, Kma7, ka8, Kma8} with {ka2, Kma2, ka1, Kma1, ka4,
Kma4, ka3, Kma3, ka6, Kma6, ka5, Kma5, ka8, Kma8, ka7, Kma7} can be regarded as the same as the original
set. Similarly, by considering the second symmetry, the parameter set that exchanged parameter{ka1, Kma1, ka2, Kma2, ka3, Kma3, ka4, Kma4, ka5, Kma5, ka6, Kma6, ka7, Kma7, ka8, Kma8, kb1, Kmb1, kb2, Kmb2, kb3,
Kmb3, kb4, Kmb4, kb5, Kmb5, kb6, Kmb6, kb7, Kmb7, kb8, Kmb8} with {kb7, Kmb7, kb8, Kmb8, kb5, Kmb5, kb6, Kmb6,
kb3, Kmb3, kb4, Kmb4, kb1, Kmb1, kb2, Kmb2, kb7, Kmb7, kb8, Kmb8, kb5, Kmb5, kb6, Kmb6, kb3, Kmb3, kb4, Kmb4, kb1,
Kmb1, kb2, Kmb2} can be regarded as the same set. By considering the third symmetry, the parameter
set that exchanged parameter {ka1, Kma1, ka2, Kma2, ka3, Kma3, ka4, Kma4, ka5, Kma5, ka6, Kma6, ka7, Kma7,
ka8, Kma8, kb1, Kmb1, kb2, Kmb2, kb3, Kmb3, kb4, Kmb4, kb5, Kmb5, kb6, Kmb6, kb7, Kmb7, kb8, Kma8} with {kb1,Kmb1, kb2, Kmb2, kb3, Kmb3, kb4, Kmb4, kb5, Kmb5, kb6, Kmb6, kb7, Kmb7, kb8, Kmb8, ka1, Kma1, ka2, Kma2, ka3,
Kma3, ka4, Kma4, ka5, Kma5, ka6, Kma6, ka7, Kma7, ka8, Kma8} can be regarded as the same set. Giventhese symmetries, the distance between two parameter sets P1 and P2 is defined as the minimum of
the distances between P1 and all the symmetrically exchangeable parameter sets of P 2. To abstract
the motif structure, parameter histograms were calculated for each cluster. When observing the
distribution of all parameters in all clusters (Figures 2 and S2), the highest frequency at whichparameters occurred with highly skewed distributions was ~0.4. Therefore, we set a threshold valueof 0.2, i.e., half the highest frequency.

Biased parameter searchParameter sets were generated by introducing some specific biases to some of their parameters;
the solutions were then classified in the same manner as for the random parameter search basedon the power spectrum. The biases were determined as combinations of motif A bias and motif B
bias. Both motif A bias and motif B bias consisted of two biased parameters, one for a high reaction
rate constant and the other for a low Michaelis−Menten constant. When the bias for the reaction
rate constant was set, the parameter was randomly generated from an exponential distribution
bounded between 100 and 1000 min−1. When the bias for the Michaelis−Menten constant was set,
the parameter was randomly generated from an exponential distribution bounded between 0.01 and0.1 μM. In this biased parameter search, we did not apply visual inspection of the sentisitivity toinitial conditions. The lack of this confirmation would not critically deteriorate the classification
performance for a chaos parameter set: for a follow-up analysis, we randomly chose 1,000
parameter sets that were classified as chaos through all the biased conditions and inspected their
sensitivity to initial conditions. Through the inspection, only 12 parameters were found to be
misclassified (in other words, 98.8% of chaos parameters were accurately classified as chaos
without the visual inspection of initial condition sensitivity).

Collection of stereotypical chaos and oscillatory parameter sets
To obtain stereotypical chaos parameter sets with symmetrical structures, we used the results of thebiased parameter search corresponding to motif arrangement #46. Parameters were randomly
sampled from a Gaussian distribution defined by the mean and the covariance matrix of chaosparameters obtained from the biased parameter search. In this sampling, distributions were notbounded within a specific range. The sampled parameter was modified so that it had a symmetricalstructure, i.e., kai = ka(9-i), kbi = kb(9-i), Kmai = Kma(9-i), Kmbi = Kmb(9-i) for i = 1…4. The solution was thendetermined to be either chaotic or not chaotic as described in the method section “Randomparameter search for chaotic and oscillatory parameter sets”. In total, 97 chaotic parameters werecollected from 10 million random sampling repeats. Chaotic parameters were confirmed to show
initial condition sensitivity. The stereotypical oscillatory parameter sets were obtained by the sameprocedure except that the numerical solution was considered to be oscillation. Collection ofoscillatory parameter sets continued until 97 sets were found.
The example chaos parameter set was selected as follows. Each stereotypical chaosparameter set was first standardized, which was conducted using the same method described forthe clustering procedure. The example parameter set was then selected as the one closest to themean of the standardized parameter sets. The distance between parameter sets was measured by
Euclidian distance. The values of the selected stereotypical chaos parameter sets are shown in
Table S1.

Bifurcation analysis
Each parameter in the typical chaotic parameter set varied between the following ranges: 1 and1000 min−1 for the reaction rate constant; 0.01 and 0.1 μM for the binding constant. The varied
value took 30 logarithmically-spaced points between the ranges, including both sides. For eachanalysis, the solution type was determined as as described in the method section “Randomparameter search for chaotic and oscillatory parameter sets”. For asymmetric perturbation, whenbifurcation analysis on each specific parameter had been conducted, the other parameters at thesymmetrical positions were fixed to the original value, resulting in a violation of symmetry. For
symmetric perturbation, both parameters at the symmetrical positions (e.g., ka1 and ka8) varied
simultaneously. Bifurcation maps were produced based on the results of the bifurcation analysis.
The color of each cell was specified using a RGB color code: the R channel was set to the
proportion of chaos behavior, G was the proportion of oscillation, and B was the proportion of
convergence.

Stochastic simulation
The stochastic simulation was performed using Gillespie’s direct method stochastic simulation
algorithm implemented in StochPy (Maarleveld et al., 2013). We modified the simulation to track the
behavior of the full system including each kinase−substrate complex molecule. The full ordinary
differential equation system, without application of the Michaelis−Menten approximation, is
described below.
d[Sa00 ]= −(𝑘𝑐a1 + 𝑘𝑐a2 )[E][Sa00 ] + 𝑘𝑢𝑐a1 [ESa00→𝑎01 ]d𝑡+𝑘𝑢𝑐a2 [ESa00→a10 ] + 𝑘a5 [FSa01 ] + 𝑘a6 [FSa10 ](11)d[Sa01 ]= −𝑘𝑐a3 [E][Sa01 ] + 𝑘𝑢𝑐a3 [ESa01 ] − 𝑘𝑐a5 [F][Sa01 ]d𝑡+𝑘𝑢𝑐a5 [FSa01 ] + 𝑘a1 [ESa00→a01 ] + 𝑘a7 [FSa11→a01 ](12)

d[Sa10 ]= −𝑘𝑐a4 [E][Sa10 ] + 𝑘𝑢𝑐a4 [ESa10 ] − 𝑘𝑐a6 [F][Sa10 ]d𝑡+𝑘𝑢𝑐a6 [FSa10 ] + 𝑘a2 [ESa00→a10 ] + 𝑘a8 [FSa11→a10 ](13)

d[Sa11 ]= −(𝑘𝑐a7 + 𝑘𝑐a8 )[F][Sa11 ] + 𝑘𝑢𝑐a7 [FSa11→a01 ]d𝑡+𝑘𝑢𝑐a8 [FSa11→a10 ] + 𝑘a3 [ESa01 ] + 𝑘a4 [ESa10 ](14)

d[ESa00→a01 ]= 𝑘𝑐a1 [E][Sa00 ] − (𝑘𝑢𝑐a1 + 𝑘a1 )[ESa00→a01 ]
d𝑡

(15)

d[ESa00→a10 ]= 𝑘𝑐a2 [E][Sa00 ] − (𝑘𝑢𝑐a2 + 𝑘a2 )[ESa00→a10 ]d𝑡
d[ESa01 ]= 𝑘𝑐a3 [E][Sa01 ] − (𝑘𝑢𝑐a3 + 𝑘a3 )[ESa01 ]d𝑡d[ESa10 ]= 𝑘𝑐a4 [E][Sa10 ] − (𝑘𝑢𝑐a4 + 𝑘a4 )[ESa10 ]d𝑡
d[FSa01 ]= 𝑘𝑐a5 [F][Sa01 ] − (𝑘𝑢𝑐a5 + 𝑘a5 )[FSa01 ]d𝑡d[FSa10 ]= 𝑘𝑐a6 [F][Sa10 ] − (𝑘𝑢𝑐a6 + 𝑘a6 )[FSa10 ]d𝑡
d[FSa11→a01 ]= 𝑘𝑐a7 [F][Sa11 ] − (𝑘𝑢𝑐a7 + 𝑘a7 )[FSa11→a01 ]d𝑡d[FSa11→a10 ]= 𝑘𝑐a8 [F][Sa11 ] − (𝑘𝑢𝑐a8 + 𝑘a8 )[FSa11→a10 ]d𝑡
d[Sb00 ]= −(𝑘𝑐b1 + 𝑘𝑐b2 )[E][Sb00 ] + 𝑘𝑢𝑐b1 [ESb00→b01 ]d𝑡+𝑘𝑢𝑐b2 [ESb00→b10 ] + 𝑘b5 [FSb01 ] + 𝑘b6 [FSb10 ]

(16)
(17)
(18)
(19)
(20)
(21)
(22)

(23)

d[Sb01 ]= −𝑘𝑐b3 [E][Sb01 ] + 𝑘𝑢𝑐b3 [ESb01 ] − 𝑘𝑐b5 [F][Sb01 ]d𝑡+𝑘𝑢𝑐b5 [FSb01 ] + 𝑘b1 [ESb00→b01 ] + 𝑘b7 [FSb11→b01 ](24)

d[Sb10 ]= −𝑘𝑐b4 [E][Sb10 ] + 𝑘𝑢𝑐b4 [ESb10 ] − 𝑘𝑐b6 [F][Sb10 ]d𝑡+𝑘𝑢𝑐b6 [FSb10 ] + 𝑘b2 [ESb00→b10 ] + 𝑘b8 [FSb11→b10 ](25)

d[Sb11 ]= −(𝑘𝑐b7 + 𝑘𝑐b8 )[F][Sb11 ] + 𝑘𝑢𝑐b7 [FSb11→b01 ]d𝑡+𝑘𝑢𝑐b8 [FSb11→b10 ] + 𝑘b3 [ESb01 ] + 𝑘b4 [ESb10 ](26)

d[ESb00→b01 ]= 𝑘𝑐b1 [E][Sb00 ] − (𝑘𝑢𝑐b1 + 𝑘b1 )[ESb00→b01 ]d𝑡
d[ESb00→b10 ]= 𝑘𝑐b2 [E][Sb00 ] − (𝑘𝑢𝑐b2 + 𝑘b2 )[ESb00−b10 ]d𝑡
d[ESb01 ]= 𝑘𝑐b3 [E][Sb01 ] − (𝑘𝑢𝑐b3 + 𝑘b3 )[ESb01 ]d𝑡
d[ESb10 ]= 𝑘𝑐b4 [E][Sb10 ] − (𝑘𝑢𝑐b4 + 𝑘b4 )[ESb10 ]d𝑡
d[FSb01 ]= 𝑘𝑐b5 [F][Sb01 ] − (𝑘𝑢𝑐b5 + 𝑘b5 )[FSb01 ]d𝑡
d[FSb10 ]= 𝑘𝑐b6 [F][Sb10 ] − (𝑘𝑢𝑐b6 + 𝑘b6 )[FSb10 ]d𝑡
d[FSb11→b01 ]= 𝑘𝑐b7 [F][Sb11 ] − (𝑘𝑢𝑐b7 + 𝑘b7 )[FSb11→b01 ]d𝑡
d[FSb11→b10 ]= 𝑘𝑐b8 [F][Sb11 ] − (𝑘𝑢𝑐b8 + 𝑘b8 )[FSb11→b10 ]
d𝑡

(28)
(29)
(30)
(31)
(32)
(33)
(34)
(35)

d[E]
= −((𝑘𝑐a1 + 𝑘𝑐a2 )[Sa00 ] + 𝑘𝑐a3 [Sa01 ] + 𝑘𝑐a4 [Sa10 ])[E] + (𝑘𝑢𝑐a1 + 𝑘a1 )[ESa00→a10 ]d𝑡+(𝑘𝑢𝑐a2 + 𝑘a2 )[ESa00→a10 ] + (𝑘𝑢𝑐a3 + 𝑘a3 )[ESa01 ] + (𝑘𝑢𝑐a4 + 𝑘a4 )[ESa10 ]−((𝑘𝑐b1 + 𝑘𝑐b2 )[Sb00 ] + 𝑘𝑐b3 [Sb01 ] + 𝑘𝑐b4 [Sb10 ])[E] + (𝑘𝑢𝑐b1 + 𝑘b1 )[ESb00→b01 ]+(𝑘𝑢𝑐b2 + 𝑘b2 )[ESb00→b10 ] + (𝑘𝑢𝑐b3 + 𝑘b3 )[ESb01 ] + (𝑘𝑢𝑐b4 + 𝑘b4 )[ESb10 ]

(36)

d[F]
= −(𝑘𝑐a5 [Sa00 ] + 𝑘𝑐a6 [Sa10 ] + (𝑘𝑐a7 + 𝑘𝑐a8 )[Sa11 ])[F] + (𝑘𝑢𝑐a5 + 𝑘a5 )[FS01 ]d𝑡+(𝑘𝑢𝑐a6 + 𝑘a6 )[FSa10 ] + (𝑘𝑢𝑐a7 + 𝑘a7 )[FSa11→a01 ] + (𝑘𝑢𝑐a8 + 𝑘a8 )[FSa11→a10 ]−(𝑘𝑐b5 [Sb00 ] + 𝑘𝑐b6 [Sb10 ] + (𝑘𝑐b7 + 𝑘𝑐b8 )[Sb11 ])[F] + (𝑘𝑢𝑐b5 + 𝑘b5 )[FS01 ]+(𝑘𝑢𝑐b6 + 𝑘b6 )[FSb10 ] + (𝑘𝑢𝑐b7 + 𝑘b7 )[FSb11→b01 ] + (𝑘𝑢𝑐b8 + 𝑘b8 )[FSb11→b10 ](37)

In these equations, [Sa00], [Sa01], [Sa10], and [Sa11] represent the concentrations of the four
phosphorylation states of substrate A; [Sb00], [Sb01], [Sb10], and [Sb11] represent the concentrations of
the four phosphorylation states of substrate B; [E] and [F] are the concentrations of free (unbound)kinase and phosphatase; [ES] and [FS] are the enzyme-substrate complex concentrations with thearrow in the subscript indicating which reaction the complex is involved in (e.g., ESa00→a01 denotesthat the ES complex is involved in the kinase reaction converting S a00 to Sa01); ka1-ka8 and kb1-kb8 are
the catalytic rate constants; kca1-kca8 and kcb1-kcb8 are the binding constants; and kuca1-kuca8 and kucb1kucb8 are the unbinding constants. The values of the parameter sets are shown in Table S2. The
relationship among the Michaelis−-Menten constants, binding constants, and unbinding constants is
described as follows:
𝐾𝑚i =𝑘𝑢𝑐i + 𝑘i, for i = 1 … 8
𝑘𝑐i

(38)

Software for computer simulation
Numerical integration and clustering were carried out by Python 3.6.1 with the libraries of numpy
1.12.1 and scipy 0.19.1, or by Mathematica 12.0 (Wolfram Research). Stochastic simulation was
carried out by StochPy 2.3.