Article

A Design Principle for an Autonomous Posttranslational Pattern Formation
Graphical Abstract

Authors
Shuhei S. Sugai, Koji L. Ode,Hiroki R. Ueda

Correspondence
uedah-tky@umin.ac.jp

In Brief
Using computer simulations, Sugai et al.
showed that a generic Michaelis-Menten
scheme of two-site reversible
phosphorylation can produce Turing
patterns in the spatial distribution of a
substrate’s modification states. A
random parameter search found typical
combinations of reaction parameters
accounting for the pattern formation andtuning the pattern shapes.

HighlightsdReversible two-site phosphorylation of a substrate canproduce spatial patternsdThis design principle is mass conserved, non-autocatalytic,
and non-allostericdCyclic reactions with biased diffusion and enzyme
sequestration are design motifsdStochastic simulation revealed two reaction-diffusion cycles
that shape patterns

Sugai et al., 2017, Cell Reports 19, 863–874
April 25, 2017 ª 2017 The Author(s).
http://dx.doi.org/10.1016/j.celrep.2017.03.081

Cell Reports

ArticleA Design Principle for an Autonomous
Post-translational Pattern Formation
Shuhei S. Sugai,1,3 Koji L. Ode,1,2,3 and Hiroki R. Ueda1,2,4,*
1Department of Systems Pharmacology, Graduate School of Medicine, The University of Tokyo, Hongo 7-3-1, Bunkyo-ku,
Tokyo 113-0033, Japan
2Laboratory for Synthetic Biology, RIKEN Quantitative Biology Center, 1-3 Yamadaoka, Suita, Osaka 565-0871, Japan
3These authors contributed equally
4Lead Contact*Correspondence: uedah-tky@umin.ac.jp
http://dx.doi.org/10.1016/j.celrep.2017.03.081

SUMMARY

Previous autonomous pattern-formation models
often assumed complex molecular and cellular networks. This theoretical study, however, shows thata system composed of one substrate with multisite
phosphorylation and a pair of kinase and phosphatase can generate autonomous spatial information,
including complex stripe patterns. All (de-)phosphorylation reactions are described with a genericMichaelis-Menten scheme, and all species freely
diffuse without pre-existing gradients. Computational simulation upon >23,000,000 randomly generated parameter sets revealed the design motifs ofcyclic reaction and enzyme sequestration by slowdiffusing substrates. These motifs constitute shortrange positive and long-range negative feedback
loops to induce Turing instability. The width andheight of spatial patterns can be controlled independently by distinct reaction-diffusion processes.
Therefore, multisite reversible post-translational
modification can be a ubiquitous source for various
patterns without requiring other complex regulations
such as autocatalytic regulation of enzymes and isapplicable to molecular mechanisms for inducing
subcellular localization of proteins driven by posttranslational modifications.
INTRODUCTION
Theoretical modeling has been a powerful tool to analyze
the mechanism of autonomous pattern formation in biologicalsystems (Karsenti, 2008; Kondo and Miura, 2010). One of thebest-known models is the Turing model deploying reactiondiffusion equations (Turing, 1952). The Turing instability can arise
from a system with two components, an activator with a slow
diffusion rate and an inhibitor with a fast diffusion rate (Gierer
and Meinhardt, 1972). The activator increases the concentrations of both the activator and the inhibitor within a short rangeof space, while the inhibitor represses the concentration of theactivator over a long range. The detailed mechanism to realize
the activation/inhibition can be varied; for example, an activator-depleted substrate scheme assumes that the inhibitory
pathway is passively mediated by the insufficient supply to produce the activator (Meinhardt, 2008).
The key features underlying the Turing model have been
discovered in the molecular and cellular mechanisms of biological pattern formation. At a cellular-circuit level, biological examples of the Turing model were reported in pattern formationof animal skin (Asai et al., 1999; Kondo and Asal, 1995) and in
the process of vertebrate morphogenesis (Economou et al.,€ller et al., 2012; Sheth et al., 2012). At a molecular2012; Mu
network level, a system consisting of a group of proteins called
MinC, MinD, and MinE is known to show typical spatiotemporal
patterns explained based on Turing instability (Raskin andde Boer, 1999; Zieske and Schwille, 2013; Loose et al., 2008). Eukaryotic cells also employ reaction-diffusion systems, including
GTPase Cdc42 and kinase-substrate network of PAR proteins
and atypical protein kinase C (aPKC) (Etienne-Manneville andHall, 2002; Hoege and Hyman, 2013) as an underlying mechanism for cellular pattern formation.
Theoretical models to simulate such molecular mechanisms ofpattern formation typically employed complex reaction networks(e.g., oligomerization and/or mutual inhibition among PAR
proteins [Dawes and Munro, 2011; Tostevin and Howard,
2008] or cell-compartment specific activation of enzymatic activities [Alonso, 2016; Alonso and Bär, 2010, 2014; Halatek andFrey, 2012; Otsuji et al., 2007]). Alternatively, a more genericset of enzymes and substrates might be sufficient for autonomous pattern formation. Analysis using the generic but not too
abstractive reaction scheme such as Michaelis-Menten scheme(Michaelis et al., 2011) allows us to compare the model parameters with the biochemically measurable values. A generic set ofcomponents to achieve reversible phosphorylation is one substrate and a kinase and a phosphatase (Kholodenko, 2006).
Despite its simple setting, it has been found that reversible phosphorylation at a single substrate site can produce an ultrasensitive response in the phosphorylation state of the substrate alongwith a linear change of kinase/phosphatase activity (Goldbeter
and Koshland, 1981). If the number of phosphorylation sites isincreased, the phosphorylation status can have two or more
distinct steady states (Markevich et al., 2004; Thomson and Gunawardena, 2009). It was also discovered that a traveling wave

Cell Reports 19, 863–874, April 25, 2017 ª 2017 The Author(s). 863
This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).

A

B

Figure 1. Autonomous Spatial Pattern Formation Arises from Reversible Dual Phosphorylation

C(A) A scheme for reversible two-site phosphorylation. The enzymatic reaction follows a MichaelisMenten reaction scheme with two parameters; kdenotes the reaction speed and Km corresponds tothe inverse of enzyme-substrate affinity. Green
arrow indicates the phosphorylation reaction. Bluearrow indicates the de-phosphorylation reaction.
(B) An example of pattern formation using thereaction scheme in (A). Horizontal dotted line indicates the time section at t = 1,000 min shown in(C). Vertical dotted line indicates the spatial section
at the position of 2.5 and 7.5 mm shown in (D). The
parameter set used in this figure is ‘‘representative
parameter set for cluster 1’’ shown in Table S1.
(C) A section of (B) at t = 1,000 min shows thespatial pattern at the steady state. Yellow-shaded
area corresponds to the S01 peak area, where S01
concentration is high.
(D) Time course of the concentration change at theposition of 2.5 and 7.5 mm.
See also Figure S1 and Table S1.

D(Markevich et al., 2006) and self-sustained oscillation (Jolley
et al., 2012) can arise from a reversible multisite phosphorylation
system. Importantly, these models of reversible multisite phosphorylation system did not assume any typical allosteric regulation of enzymes (i.e., an enzyme is regulated by a molecule that is
not a direct substrate of the reaction but acts on the enzyme
through a site different from a catalytic pocket or a substrate
recognition site), changes in the amount of protein through protein synthesis and degradation, or any special function of enzymes and substrate (e.g., autocatalytic activity), suggesting
that, under the proper combination of reaction parameters, a
simple reversible phosphorylation reaction can be the core ofcomplex and autonomous behaviors in biological systems.
A series of studies has investigated that the reversible andmultisite post-translational scheme can amplify and maintainthe spatial information in the presence of pre-existing patterning
such as spatial compartmentalization/gradient of preferential
reaction, components, and external signals (Alam-Nazki andKrishnan, 2012, 2013, 2015), or with the ability of similar reversible post-translational systems to process spatially inhomogeneous stimulus (Krishnan, 2009; Seaton and Krishnan, 2011).
Nonetheless, it had been unclear whether a stable spatial
pattern can autonomously arise from reversible multisite phosphorylation in the absence of any pre-existing patterning. In
this study, we show that reversible two-site phosphorylation
described as Michaelis-Menten-type reaction processes is sufficient for autonomous and stable spatial pattern formation. A
random parameter search revealed the core mechanisms forautonomous pattern formation; local enzymatic sequestrationby the interaction between enzymes and slow-diffusing substrates inhibits the production of fast-diffusing substrates. The
sequestration creates an inflow of fast-diffusion substrates toward a specific area, maintaining spatial heterogeneity. Stochastic simulation further revealed the presence of two types
864 Cell Reports 19, 863–874, April 25, 2017of reaction-diffusion reaction cycles that operate to maintain alocal accumulation of substrate with specific phosphorylation
states.
RESULTS
Autonomous Spatial Pattern Formation Arises from
Reversible Dual Phosphorylation
To construct a simple pattern-formation system, we consider aminimal model that comprises a kinase, a phosphatase, and theirsingle substrate. It was analytically shown that, if the substrate
has only one phosphorylation site, this system will always
converge to a temporally and spatially uniform steady state,
and a spatial pattern will not emerge autonomously (Supplemental Experimental Procedures). Therefore, we considered adual-phosphorylation system where the two phosphorylationsites on the substrate are phosphorylated by a common kinaseand de-phosphorylated by a common phosphatase (Figure 1A).
In this model, the four possible modification states will be called
as S00, S10, S01, and S11, respectively; S00 accounts for dephosphorylated substrate, S10 and S01 account for mono-phosphorylated substrate, and S11 accounts for di-phosphorylated
substrate. Each chemical species is (de-)phosphorylated anddiffuses with certain reaction and diffusion rates (Figure 1A).
Mass-action kinetics is applied to one-to-one binding of enzyme
and substrate and the phosphate transfer reactions. To simplify
the simulation process, a Michaelis-Menten approximation isapplied, but a more detailed model without this approximation
will be considered in the later sections. Each substrate diffusesin one-dimensional (1D) space with circular boundary condition
following Fick’s law of diffusion. The 1D system size was set to
10 mm corresponding to cellular/subcellular scale. Consistent
with our previous study (Jolley et al., 2012), particular parameter
sets showed oscillatory behavior (Figures S1A and S1B). WeACFigure 2. Three Design Motifs in Parameters
for Spatial Pattern Formation

B

D

E

further found that reversible two-site phosphorylation can induce
autonomous pattern formation under certain parameter combinations (Figures 1B–1D).
Three Design Motifs in Parameters for Spatial Pattern
Formation
To elucidate the underlying design motifs for spatial pattern formation, a comprehensive random parameter search was conducted (Figure 2A). Reaction rate constants (k), Michaelis constants (Km), and diffusion coefficients (D) were randomly chosenfrom 1–1,000 min–1, 0.01–1,000 mM, and 0.001–100 mm2 min–1(1.67 3 10 5-1.67 mm2 s–1), respectively. These ranges include
biochemically probable reaction and diffusion velocities based(A) Workflow for parameter motif identification.
(B) Clustering of the spatial parameter sets with
Ward’s algorithm.
(C–E) Parameter histograms of the spatial parameter sets. For each cluster, schematic representation is drawn at the top. Fast reactions are represented with thick arrows. Substrates with highest
affinity to kinase (E) and phosphatase (F) are indicated by attaching the corresponding enzymes tothe edge. Diffusion coefficients are represented as
the labels ‘‘fast’’ or ‘‘slow.’’
See also Figure S2.

on the previous studies (Arrio-Dupont
et al., 2000; Jolley et al., 2012; Markevich
et al., 2004). As a result, about 0.05% and0.08% of the total parameter sets yielded
spatial patterns and oscillatory behaviors,
respectively. These parameter sets werethen clustered by Ward’s algorithm. The
reaction-diffusion scheme in our model isunchanged if the first and second phosphorylation sites are exchanged (e.g.,exchanging S01 and S10) and if the kinaseand phosphatase are exchanged (e.g.,exchanging S00 and S11). To take such
symmetries into account and extract qualitatively different clusters, the distance matrix for the clustering was made to incorporate symmetrical reflection and rotation.
The resulting clustering dendrogram suggested the existence of two clusters foroscillatory parameter sets (Figure S2A)
and three clusters for spatial parameter
sets (Figure 2B).
Two clusters in oscillatory parameter
sets commonly have two design motifs
consistent with the previous study (Jolley
et al., 2012) (Figures S2B and S2C). In
brief, the biochemical system in oscillatory parameter sets tends to have clockwise cyclic reactions (S00/S01/S11/
S10). In addition, a particular substrate
tends to sequestrate an enzyme due to its strong substrateenzyme affinity, leading to synchronization of substrate states.
For example, low Km5 results in strong binding of S01 to phosphatase in cluster 1. Therefore, S01 sequestrates phosphatase andinhibits the S11/S10 step, leading to synchronization of substrate states at S11. Cluster 2 can be interpreted as a 90-degree
rotation of cluster 1. These ‘‘reaction’’ and ‘‘binding’’ motifs
found in the oscillatory clusters are similarly conserved in spatial
clusters 1 and 2 (Figures 2C and 2D). Notably, a bias of diffusion
rates was found in spatial parameter sets (‘‘diffusion’’ motif) but
not in oscillatory parameter sets. To further investigate the rolesof the three design motifs in spatial pattern formation, we focus
on cluster 1 in the following sections without much loss of
Cell Reports 19, 863–874, April 25, 2017 865

A

B

C

D

E

F

generality because cluster 1 is the largest cluster, and because
cluster 1 and its rotationally symmetric cluster (cluster 2)
together account for more than 85% of all spatial parameter
sets. Mechanism of pattern formation in cluster 3 (Figure 2E)
will be revisited in the later section.
Necessity of Three Design Motifs
Three design motifs in cluster 1 can be represented as follows(Figure 3A). For the ‘‘reaction’’ motif, reaction rates (k1 and k8)
are biased toward the clockwise direction. For the ‘‘binding’’
motif, Michaelis constants (Km4 and Km5) are biased to promotekinase sequestration by S10 and phosphatase sequestration by
S01, preventing the other species from binding the enzymes.
For the ‘‘diffusion’’ motif, diffusion rates (D) of the enzyme-sequestrating substrates (i.e., DS01 and DS10) are extremely slow,
866 Cell Reports 19, 863–874, April 25, 2017

Figure 3. Inner and Outer Reaction-Diffusion Cycles Account for Peak Shape(A) Schematic representation of reaction motif (k)for cyclic reaction direction, binding motif (Km) for
enzymatic sequestration, and diffusion motif (D) for
biased diffusion velocity.
(B) Applying a constraint on the condition(s) ofrandom parameter search corresponding to thethree design motifs improves the rate of finding
parameter sets for spatial pattern formation.
(C) Schematic representation of the ‘‘inner’’ and‘‘outer’’ cycles.
(D) Spatial distribution of kinase-driven reaction
rates (upper half) and phosphatase-driven reaction
rates (lower half) at the steady state of spatial
pattern formation. The simulation is conductedusing the representative typical parameter set for
spatial cluster 1 (Table S1). The yellow-shaded
area corresponds to the S01 peak area.
(E) Three significant pattern-determining factors(k3, k5, and DS01) were altered one by one in bifurcation analyses from the original values in the
representative parameter set shown by the red
dotted lines. For each value on the horizontal axis,
the two dots represent the peak and trough concentration of S01 in the resultant spatial pattern.
The color of the dots stands for the width of thepeak. Black dot means spatially uniform steadystate. ‘‘Oscillation’’ indicates the presence oftemporal oscillation. In the other parameter areas
with no dot or indication, the calculation failed(diverged) probably because the parameter value
is close to the bifurcation point.
(F) Example distributions of substrates arising frommodified parameter sets. The indicated parameter
was either increased 3-fold or decreased to
one-third compared with original representative
parameter set.
See also Figure S3 and Table S1.

while those of the other substrates (i.e.,DS00 and DS11) are fast. Biased parametersearches with constraints reflecting thesemotifs resulted in an improved hit rate
for spatial patterns (Figure 3B). When all
three types of bias were implemented,
the hit rate synergistically increased to as much as 1.86%, amore than a 30-fold increase from the original rate, suggestingthat these three design motifs work cooperatively in spatial
pattern formation.
Inner and Outer Reaction-Diffusion Cycles Account forPeak ShapeTo elucidate the underlying reaction-diffusion mechanism of
spatial pattern formation, we investigated a representative stereotypical parameter set belonging to spatial cluster 1 (Table S1).
There are two characteristic areas in the steady-state distribution,
where S01 and S10 are abundant (i.e., S01 and S10 peak areas) (seeFigure 1C). As expected from the distribution pattern, certain reactions tend to occur in a limited range of space at the steady state(Figure 3D). By inspecting such spatially restricted reactions, wehypothesized that there are two types of reaction-diffusion cycles,
namely, the ‘‘inner’’ and ‘‘outer’’ cycles (Figure 3C). The ‘‘inner’’ cycle can be defined by S01/S00/S01 (occurring at the yellow
shaded area in Figure 3D). Once an S01 substrate is converted
to S00, S00 eventually returns to S01 and closes the ‘‘inner’’ cycle,
since there is almost no reaction in the direction of S00/S10
throughout 1D space. We also note that another ‘‘inner’’ cycle
can be similarly defined by S10/S11/S10.
On the other hand, the ‘‘outer’’ cycle can be defined by S01/
S11/S10/S00/S01. This ‘‘outer’’ cycle starts from the reaction
S01/S11 that predominantly occurs at the S01 peak area (Figure 3D, S01/S11). Note that there is almost no reaction of thedirect conversion of S11/S01 throughout 1D space. Then, S11
likely undergoes a transportation from S01 peak area to S10
peak area, because S11 is abundant at the S01 peak area anddepleted at the S10 peak area (Figure 1C). The subsequent reaction S11/S10 occurs at the S10 peak area (Figure 3D, S11/S10).
A similar explanation is applicable for the reaction-diffusion process S10/S00/S01.
To extract the decisive parameters for peak height and width,
we conducted comprehensive bifurcation analyses (Figures 3E
and S3). We found three parameters (k3, k5, and DS01) (Figure 3E)
and their symmetrical counterparts (k6, k4, and DS10) (Figure S3A)
are markedly important for determining peak height and width.
For example, a reaction parameter (k3) of the first reaction inthe ‘‘outer’’ cycle largely affects peak width in a negative manner(Figure 3F, left). On the other hand, parameters (k5, and DS01)
involved in the ‘‘inner’’ cycle largely affects peak height: diffusionof S01 (DS01) attenuates the ‘‘inner’’ cycle by facilitating theoutflow of S01 from its peak, and its decrease leads to a higher
peak (Figure 3F, right). Similarly, peak growth is also facilitatedby attenuating the outflow of S01 from its peak by accelerating
k5 (Figure 3F, middle). In summary, stimulation of the ‘‘outer’’ cycle carries substrates away from S01 peak to the other area, leading to a narrower peak. By contrast, stimulation of the ‘‘inner’’ cycle keeps substrates inside the peak and leads to a higher peak.
Stochastic Simulations Directly Track the ‘‘Inner’’ and‘‘Outer’’ Reaction-Diffusion Cycles
To validate the existence of ‘‘inner’’ and ‘‘outer’’ cycles, we next
aimed to directly observe the biochemical mechanisms underlying
the spatial pattern formation by conducting a simulation using astochastic simulation algorithm (SSA) (Gillespie, 1977), which
can track individual molecules and their reactions at each position
and time. Note that substrate-free enzymes and substrateenzyme complexes are not explicitly considered in the presenceof Michaelis-Menten approximation. To conduct SSA simulations,representative parameter set for cluster 1 was converted tomass-action parameters without Michaelis-Menten approximation (Table S2). For the sake of simplicity, all diffusion rates for substrate-free enzymes and all types of enzyme-substrate complexes
are chosen to be the center of the parameter range used for
random parameter search (i.e., 0.01667 mm2 s–1).
With this setting, the total number of molecules was changed
to search for the minimal number of molecules that can yield aspatial pattern (Figure 4A). A simulation with 100 enzyme and150,000 substrate molecules results in a stable spatial patternthat remains at the same spatial position for nearly 1,000 min.

The distribution of each substrate averaged throughout
100 min of simulation time with SSA (Figure S4A, left column)
is mostly similar to those of the deterministic partial differential
equation (PDE) (see Figure 1C). In addition, diffusion rates ofsubstrate-free enzymes and enzyme-substrate complexes do
not qualitatively change the steady-state spatial pattern at least
for the representative parameter set, suggesting that the spatial
pattern formation observed in the PDE model is not an artifact
caused by Michaelis-Menten approximation (Figure S4B).
Surprisingly, spatial pattern is still visible in the simulation witha single enzyme and 100 substrate molecules. Although thepeak position of S01 fluctuates, the distribution pattern across
1D space at defined time shows a spatial heterogeneity qualitatively different from random noise arising from a parameter set
without any design motifs for spatial pattern formation (Figures
4A and 4B). Because many reaction events are required to calculate the histogram of reaction rates, we analyzed the simulation
with 100 enzyme and 150,000 substrate molecules in thefollowing analysis. Nevertheless, the observed pattern formation
with a small number of enzyme and substrate moleculessuggests that pattern formation by reversible multisite phosphorylation can effectively work in subcellular compartmentswith a limited number of molecules.
To verify the existence of the ‘‘inner’’ and ‘‘outer’’ cycles, wenext calculated the distribution of each enzymatic reaction anddiffusion events in the SSA simulation (Figure S4A). Detailed inspection of these distributions confirmed the existence of the ‘‘inner’’ and ‘‘outer’’ cycles (Figure 4C). For example, ES00/01 andFS01/00 are both enriched within the S01 peak area (Figure S4A,
left column), which corresponds to the ‘‘inner’’ cycle (S01/S00/
S01) predominant in S01 peak area. On the other hand, the first halfof the ‘‘outer’’ cycle (S01/S11/S10) starts at S01 peak area withthe formation of ES01/11 to produce S11 (Figure S4A, left column). The fast-diffusing S11 is then transported from the S01
peak area to the S10 peak area (Figure S4A, right column) and results in FS11/10 formation, leading to S10 production (Figure S4A,
left column). Similarly, the second half of the ‘‘outer’’ cycle (S10/S00/S01) continues at the S10 peak area with the formation of
FS10/00 to produce S00 (Figure S4A, left column). The fastdiffusing S00 is then transported from S10 peak area to S01 peakarea (Figure S4A, right column) and results in ES00/01 formation,
leading to S01 production (Figure S4A, left column).
To more directly observe the ‘‘inner’’ and ‘‘outer’’ cycles, we
tracked a single molecule over a long trajectory (Figure 4D).
When an S01 molecule is converted to S00, the subsequent enzymatic reaction is almost always the conversion to S01 inside S01
peak area (‘‘inner’’ cycle). On the other hand, when the S01 molecule is converted to S11, the subsequent S10 production tends tooccur outside the S01 peak area (the first half of the ‘‘outer’’ cycle).
Three Design Motifs Cooperatively Contribute to the‘‘Inner’’ and ‘‘Outer’’ Cycles to Generate Spatial Pattern
We then sampled 200 trajectories in the representative parameter set. In the case of trajectories for S11 and S00 molecules,
there are four possible trajectories as shown in Figure 5A. The
trajectories other than the ‘‘inner’’ and ‘‘outer’’ cycles are termedas ‘‘other pathways.’’ The re-production of S01 through the ‘‘inner’’ cycle (S01/S00/S01) mostly occurs within the S01 peak
Cell Reports 19, 863–874, April 25, 2017 867

Total S = 150000, E = F = 100

with spatial parameter set

Number of S01 molecules

Number of S01 molecules

2000

500

1000
500

0
0

2

4
6
8
Position (μm)

0

0

4
6
8
Position (μm)

210Total S = 100, E = F = 1

with spatial parameter set

parameter set without motifs
Number of S01 molecules
60

40
500
20

0

2

4
6
8
Position (μm)

1000

0

2

4
6
8
Position (μm)

C

10

4
2
0

0

0

3
2
1
0

0

S00
S01

P

E
S10

8
4
6
Position (μm)

S00

S10

S01

S11

10

(A) Time evolution of S01 concentration based on
the stochastic simulation. Total number of molecules altered in three ways as shown in the blackbox (top left, top right, and bottom left). A simulation with a parameter set that does not retain any
motifs for spatial parameter set is also shown in thebottom right.
(B) Distribution of each substrate at t = 100 min with
the indicated condition for the total amount of
substrate and enzymes and the parameter set.
(C) Overall event following the ‘‘inner’’ and ‘‘outer’’cycles. The cyclic reaction bias and enzyme
sequestration underlie the mutually exclusive
relationship between S01 and S10. Phosphatasesequestration by S01 inhibits the production of S10
and S00 in the S01 peak area. Absence of S10 in thisarea allows the use of kinase to produce S01 from
S00, to further promote the phosphatase sequestration. The efficient consumption of S00 creates aconcentration gradient that promotes the influxof fast-diffusing S00 from the surrounding area. A
similar story is valid for the symmetrical S10 peak.
(D) Single-molecule trajectory based on the stochastic simulation with spatial parameter set.
See also Figure S4 and Table S2.

400

S111

P

Time (min)

F P
S01

2

Figure 4. Stochastic Simulations Directly
Track Inner and Outer Reaction-Diffusion
CyclesS11
P

se Kin
qu as
es e
tra
tio
n

Ph
se ospqu ha
es tas
tra e
tio
n

S11

Fast diffusion

10

4

D

P

4
6
8
Position (μm)

parameter set without motifs

S10

P

2

Total S = 100, E = F = 1

1

0

6

3

500

0

10

with spatial parameter set

2

0

S11

S01

Total S = 100, E = F = 1

105

Total S = 1500, E = F = 1

1000
Time (min)

500

0

10

Number of S01 molecules

0Time (min)

1500

1000

S10

S00

Number of moleculeswith spatial parameter set

Time (min)Time (min)

1000

B

Total S = 100, E = F = 1

Number of molecules

A

200

S00
Fast diffusion
ff

S00

0
0

area (Figure 5B). On the other hand, the first half of the ‘‘outer’’cycle (S01/S11/S10) starts in the S01 peak area but ends in
the S10 peak area (Figure 5B).
The essential roles of ‘‘reaction,’’ ‘‘binding,’’ and ‘‘diffusion’’motifs were examined by removing each one of them from thestochastic simulation model (Table S2). The model without the‘‘reaction’’ motif failed to generate spatial patterns, with the position of each unit reaction having become totally random (Figure 5C). Furthermore, all reactions were restricted to the ‘‘inner’’
cycle in this simulation. A similar result was observed whenthe ‘‘diffusion’’ motif was removed (Figure 5D). These results
indicate that both ‘‘reaction’’ and ‘‘diffusion’’ motifs are required
for the establishment of the ‘‘outer’’ cycle (S01/S11/S10). The‘‘reaction’’ motif is important for the ‘‘outer’’ cycle, because the‘‘outer’’ cycle requires high concentration of S01 molecule thatis supported by a high reaction rate of S00/S01 in the ‘‘reaction’’motif. The ‘‘diffusion’’ motif is also important for the ‘‘outer’’ cycle(S01/S11/S10) because the ‘‘outer’’ cycle requires free kinasemolecules, which are easily trapped by S10 without a low diffusion rate of S10 in the ‘‘diffusion’’ motif.
868 Cell Reports 19, 863–874, April 25, 2017

In contrast, the model without the‘‘binding’’ motif preserved the ‘‘outer’’cycle but still failed to generate spatial
patterns (Figure 5E). The model withoutthe ‘‘binding’’ motif allows ‘‘other pathways’’ (S01/S11/S01 or S10/S11/
S01), which are prohibited with spatial
2
4
6
8
10
parameter sets, suggesting that the
Position (μm)‘‘binding’’ motif can exclude other reactions than ‘‘inner’’ or ‘‘outer’’ cycles.
Note that the model without three motifs
altogether led to the complete collapse of reaction order, witha similar number of reactions in every direction (Figure 5F).
Three Design Motifs Contribute to Short-Range Positive
Feedback and Long-Range Negative Feedback Loops toGenerate Spatial Pattern
The previous theoretical studies propose that autonomous spatial
pattern formation requires a short-range positive feedback anda long-range negative feedback (Gierer and Meinhardt, 1972;
Turing, 1952). In our model, the ‘‘inner’’ and ‘‘outer’’ cycles seems
to act as short-range positive feedback and long-range negative
feedback loops, respectively (Figure 6A). In the case of shortrange positive feedback loop, S01 indirectly activates its own production through the ‘‘inner’’ cycle. In detail, the reaction S00/S01
in the ‘‘inner’’ cycle is mediated by the kinase that is strongly trapped by S10. Alternatively, S01 strongly binds to the phosphataseand effectively represses the production of S10 from S11. Thus,
S01 and S10 are mutually exclusive each other, creating a positive
feedback. Because S01 is the slow-diffusing molecule, thispositive feedback loop acts within a short range. The efficientFigure 5. Three Design Motifs Cooperatively
Contribute to Inner and Outer Cycles toGenerate Spatial Pattern

A(A) Diagrams explaining the reaction track classified as the ‘‘inner’’ cycle, the ‘‘outer’’ cycle, and‘‘other pathways.’’(B) Transition diagram of individual molecules
derived from cluster 1-representative parameter
set. The tracks are extracted from the reactions
that occurred during 50 min in the steady state. To
highlight the exclusive existence of ‘‘inner’’ and‘‘outer’’ cycle and the absence of any other transitions, only the transitions following the triplet order shown in the box are selected for drawing. For
each pathway, the fast-diffusing species’ movements are represented as the lines, whose startingand ending points are shown as dots and triangles,
respectively. In this condition, no reaction was
found to follow the ‘‘other pathways.’’(C–F) The parameter sets are changed to the ones
with collapsed motifs. For (C) and (D), there are notransitions in the direction of the ‘‘outer’’ cycle and‘‘other pathways.’’See also Table S2.

B

C

D

E

F

S00/S01 reaction promotes the consumption of S00 at this area.
The S00 gradient results in the inflow of fast-diffusing S00 toward
the S01 peak area, further promoting the S00/S01 reaction. In
the case of a long-range negative feedback loop, we can assume
two types of antagonizing effect to the production of S01. The first
type is the depletion of fast-diffusing substrate molecule S00.

Indeed, the level of S01 production is suppressed at the center of S01 (Figure 3D),
and, when the width of S01 increased, theS01 maxima started to split apart (Figure 3F). This behavior is typically observedin an activator-depleted substrate model(Meinhardt, 2008). The second type is theindirect suppression of S01 production
through the ‘‘outer’’ cycle. In detail, S01
indirectly activates the production of S10
through the ‘‘outer’’ cycle, in which S01 is
first converted to S11 in the S01 peak area.
S11 then diffuses rapidly from S01 peakarea to S10 peak area. In the S10 peak
area, S11 is eventually converted into S10.
Since S10 strongly binds to a kinase andeffectively represses S01 production from
S00, the ‘‘outer’’ cycle can form the negative feedback loop for S01. Because S11 is
the fast-diffusing molecule, this negative
feedback loop acts over a long range.
A Reversible Two-Site Modification
Model Can Generate Complex
Pattern Formation in TwoDimensional SpaceThe Turing mechanism is well known toelicit various spatial patterns when it operates in two-dimensional (2D) space. To verify the potential of areversible two-site modification model, we extended our simulation in 2D space with periodic boundary condition. The area was
set to 900 mm2 corresponding to the order of eukaryotic cell
membrane (Kholodenko et al., 2000). Because a representative
parameter set in cluster 1 has a symmetric and mutually
Cell Reports 19, 863–874, April 25, 2017 869

A

B

C

D

E

exclusive relationship between S01 and S10, it exhibited a stripe
pattern when it operates in 2D space, where S01 and S10 complementarily occupied the space with almost identical stripe widthand peak concentration (Figure 6B). When we changed the valueof k3, the critical parameter for the ‘‘outer’’ cycle and hence forpeak width, the spatial pattern in 2D space was greatly affected:
the increase of k3 resulted in the decreased stripe width andended up with the fragmented dots of high S01 regions (Figure 6C). Along the 100-fold change in k3 value, the peak valueof S01 was less affected. On the other hand, changes in k5 andDS01, parameters responsible for the ‘‘inner’’ cycle, and hence
the peak height had more impact on the peak concentration of
stripe. The 100-fold increase in the k5 value and 100-fold
870 Cell Reports 19, 863–874, April 25, 2017

Figure 6. A Reversible Two-Site Modification Model Can Generate Complex Pattern
Formation in 2D Space(A) The ‘‘inner’’ cycle with phosphatase sequestration contributes to the short-range positive
feedback for the production of S01. The ‘‘outer’’cycle with kinase sequestration contributes tolong-range negative feedback for the productionof S01.
(B) 2D patterns from cluster 1-representative
parameter set.
(C–E) Modifying one of the parameters in thecluster 1-representative parameter set resulted invarious pattern shapes. Only the distributions of
S01 are shown. Changes in k3 predominantly affectthe peak width. The maximum concentration of
S01 is altered significantly with k5 and Ds01 (see theconcentration shown below each 2D plot).
See also Figure S5 and Table S1.

decrease in the DS01 value changed thepeak concentration almost 10-fold (Figures 6D and 6E), while it did not affectthe stripe width much. A similar mechanism of pattern formation is likely to takeplace for the parameter sets of cluster 2(Figure S5).
Breakdown to the Minimal
Requirement for Spatial Pattern
Formation
The remaining cluster 3 had horizontal
reflection symmetry (i.e., the overall
parameter distribution is unchanged byexchanging the position of S01 and S10)
(Figure 2E). Even when one of the two
mono-phosphorylated substrates wasomitted from the scheme (Figure 7A), themodel with three modification states (i.e.,
S0, S1, and S2) exhibited spatial patternformation. A random parameter search revealed one cluster for the spatial parameter set (Figure 7B). Parameter distributionof the cluster was almost the same aswhat was obtained by eliminating one of
the mono-phosphorylated forms from spatial cluster 3 in thefour-modification-state model (cf. Figures 2E and 7C). Therefore,
spatial cluster 3 in the four-state model seems a superimposedversion of two three-state models.
A core feature of design motifs as well as mechanism ofpattern formation found in the four-state model is conserved in
the three-state model: there is a clear reaction bias to produce
the slow-diffusion species (S0) from the fast-diffusion species(S2). Based on the parameter distribution (Figure 7C), we createda representative parameter set for this three-state model (Table
S3). A steady-state distribution of substrate concentration andreaction rates (Figure 7D) indicated that the slow-diffusing
species S0 amplifies its own synthesis via the inflow of theA

B
Km3, k3

S0

Km1, k1

DS0

S1

Km4, k4 P

P

Km2, k2

DS1

S2

Figure 7. Breakdown to the Minimal
Requirement for Spatial Pattern Formation

Random parameter sets
5,407,000

P

Spatial parameter sets
1,600

DS2

C
E P
S1

S0

P

3

0.0

Km4

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

E

4
2
0

S0
S1

-2

S2

10

S0 to S1
S1 to S00130

Overall production

20

Number of events per minute

Reaction rate
(x100 μM/min)

3

−5
−3
−1
1
Log10 D (μm2 sec-1)

4
2
0
-2

30

Reaction rate
(x100 μM/min)

0.4

Km3

Log10 of theaveraged count

1
2
Log10 k (min-1)

DS0
DS1
DS2

0.8

0.4

Km2

Rightward movement

Log10 of theconcentration (μM)

0

D

Km1

0.0

0.0

Probability
0.2

k1
k2
k3
k4

P

(Fast)

0.2

0.4(Slow)

S2

4
0
-4
-8

S0
S1
S2

(A) Three-state reaction scheme.
(B) Workflow to derive the condition for spatial
pattern formation.
(C) Parameter histograms for the spatial parameter
set and schematic representation of spatial
parameter set in three-state model.
(D) Spatial pattern arising from the representative
parameter set for three-state model. Concentrationof each substance (top), reaction rate involving
kinase (middle), and phosphatase (bottom) areplotted for each position. The shaded area corresponds to the S0 peak area. See Table S3 for the
representative parameter values.
(E) Result from stochastic simulation with thethree-state scheme and the total initial number
S = 150,000, E = F = 100. Averaged distribution(top), overall production rate (middle), and overall
rightward movement (bottom) are shown for eachbin position. Distribution was averaged for 100 min
after the establishment of spatial pattern. In thesame time frame, overall production rate andoverall movements were calculated in the samemanner as in Figure S4A. See Table S4 for theparameter values.
(F) Schematic expression of the spatial patternmaintenance. Kinase sequestration by S0 at thecenter of the peak and its release at the peripheralarea, and overall absorption of substance in theform of fast-diffusing S2 are depicted.
(G) Parameter sets that met algebraic conditions
for Turing instability (red) and that produced spatial
patterns in the simulation (blue) are projected on
the same principle component (PC) space.
See also Figure S6 and Tables S3 and S4.

60
40
20
0
-20
-40
-60

PC 2

2

0

S1 to S2
S2 to S1

4case of four-state model, the long-range
negative feedback in the three-state
110
model can be interpreted as activatordepleted scheme (Meinhardt, 2008):
0
2
4
6
8
10
0
2
4
6
8
10
Position (μm)
slow-diffusing S0 is produced from fastPosition (μm)
FGdiffusing S2 and production rate of S2 atthe center of S0 is suppressed.
Enzymesequestration
An apparent difference in the parameterS0
inhibits S2
distribution between the three-state
production
model and the four-state model lies in
Algebraic
the binding affinity of substrate andSimulatedP
PPP
enzyme. In the four-state model, theS2
S2highest substrate-enzyme affinity lies in a
Fast diffusion
slow-diffusion substrate making the peak
4
6
2
0
-2
(Figures 2C and 2D). In the three-state
PC 1
model, S1, the second slowest diffusion
S0 peak area
species has the highest affinity with kinase. This represents the restriction of
fast-diffusing species S2 (Figure 7E). At the center of S0 peak, the three-state model; n the three-state model, reproduction of S2
overwhelming amount of S0 sequestrates the kinase, inhibiting requires kinase-driven phosphorylation of S1. To this end, S1
S1/S2. This makes the overall reaction inclined toward S2/ needs to take back kinase from S0 at some point of reaction
S1/S0, and S2 concentration tends to be lower. Following this cascade, thus needs to have a high kinase-substrate affinity.
The simplicity of this three-state model enabled an analytgradient, fast-diffusing S2 flows into the center of S0 peak thatfurther stimulates the production of S0 (Figure 7F). Like the ical approach using linear stability analysis (Supplemental
120

Cell Reports 19, 863–874, April 25, 2017 871

Experimental Procedures). One requirement for the emergenceof spatial pattern is that the phosphorylation and de-phosphorylation of two sites are carried out by the same kinase and phosphatase, indicating the importance of enzymatic sequestration
to affect other reaction steps sharing the same enzyme. Furthermore, analytical examination revealed the requirements of theparameter sets to have a potential to generate spatial patterns
by amplifying small fluctuations in an almost homogeneous 1D
space. The parameter space satisfying such analytic requirements are closely overlapped with the parameter space that resulted in spatial pattern formation in numerical simulation,
implying the validity of this analytical approach (Figure 7G).
Unlike clusters 1 and 2, cluster 3 does not have a mutually
exclusive relationship between two species (e.g., S01 and S10in cluster 1). Therefore, there is no constraint to partition the
2D space for the two antagonizing species. This relaxed
constraint leads to a single peak pattern even in the 2D space.
The width of resultant peak was tunable from a small single
dot to a more blurred concentration gradient (Figure S6).
DISCUSSION
Comparison with Other Mass-Conserved ReactionDiffusion Models for Scale-free Pattern Formation
When the total amount of substances is kept constant byassuming that no de novo synthesis and degradation of components happen, the reaction-diffusion system results in the emergence of spatial pattern with single peak in 1D space (Goryachev
and Pokhilko, 2008; Ishihara and Kaneko, 2006; Otsuji et al.,2007). Under the local positive feedback for peak emergence,
the highest peak consumes the resource more effectively thanlower peaks, further stimulating the growth of the highest peak.
The winner-take-all situation allows only one peak to survive at
steady state in an entire space where the competition of limited
resources effectively occurs. Consistent with the previous
finding (Otsuji et al., 2007), most (>80%) of the spatial parameters found in our random parameter search in 1D space resulted
in the single peak for substrate with a specific modification state.
The single peak pattern, which is resilient to the changes in thesystem size, was also observed in a reaction-diffusion systemwith a mechanism called wave-pinning (Jilkine and EdelsteinKeshet, 2011; Mori et al., 2008). In the wave-pinning system,
spatially homogeneous steady state is stable, whereas in the
mass-conserved Turing system, spatially homogeneous steadystate is unstable and an infinitesimally small perturbation to break
the homogeneity leads to the emergence of spatial pattern. The
spatial pattern formation in our model belongs to the massconserved Turing system and not to the wave-pinning system,
because part of our criteria for the screening of spatial parameter set requires Turing instability. Nonetheless, it is expectedthat the scheme of reversible-multisite phosphorylation can
exhibit wave-pinning mechanism because the scheme is massconserved and can have a bistability (Markevich et al., 2004).
General Model of Pattern Formation by Reversible
Two-Site Phosphorylation
Our study proposes that the reversible two-site modification of asingle substrate by two competing enzymes could be the direct
872 Cell Reports 19, 863–874, April 25, 2017source for spatial pattern formation. For example, it was shownthat a two-phosphorylation mitogen-activated protein kinase(MAPK) signaling forms a focus in a unicellular organism fission
yeast (Dudin et al., 2016). The focalization process is assumed
to involve several feedback pathways including actin stabilization, but our model suggests that the reversible phosphorylationof MAPK itself may play a role to evoke a local accumulation of
the signaling pathway.
A part of ‘‘reaction,’’ ‘‘binding,’’ and ‘‘diffusion’’ motifs can befound in enzymes/substrates including mitogen-activated protein kinase kinase (MAPKK) and CKId/ε. These kinases phosphorylate substrates at multiple sites following a preferred order,
corresponding to the cyclic bias revealed in our modeling study(Ferrell and Bhatt, 1997; Flotow et al., 1990; Vanselow et al.,2006). In these systems, previous studies have suggested thatenzyme sequestration plays an important role to elicit non-linear
responses (Jolley et al., 2012; Markevich et al., 2004). In these
enzymes, the phosphorylation status of the substrate affects
the substrate affinity through a substrate recognition site, whose
interaction with the substrate is dependent on the electrostatic
status of the substrate. The catalytic constant is also affected
because the substrate recognition can define the conformationof the substrate toward the catalytic pocket and hence determine its catalytic constant. In this sense, our biochemical model
can naturally accommodate the characteristics of these enzymes, which do not have a typical allosteric regulation. Anotherspecific requirement for the formation of spatial patterns is adifference in diffusion rates depending on the phosphorylation
state of the substrate. Although the change in molecular weightby the addition of one or two phosphate groups itself will not alter
the diffusion rate of modified substrate significantly, the diffusion
rate can be dramatically altered by the altered affinity of modified
substrate with scaffold structures such as cytoskeletal proteins,cell membranes, and chromosomes.
Although our modeling study started without constraints
based on previously identified molecular networks, the overallmechanism encompasses several key features found in otherspecific biological models including sequestration-based directional flow and mutual exclusion of two key species (Dawes andMunro, 2011; Halatek and Frey, 2012; Otsuji et al., 2007; Tostevin and Howard, 2008). Furthermore, our model can provide a
unified view for the role of protein localization; when the substrate-enzyme complex is localized in the specific area, therole of localization may be enzyme sequestration for autonomous spatial formation rather than only enrichment of enzymatic
activity. Taken together, the simplicity and generality of ourreversible phosphorylation model provides a framework toinvestigate a biochemical source of autonomous pattern formation and ultimately, to synthesize a de novo molecular networkthat inherently encodes spatial information.
EXPERIMENTAL PROCEDURES
All the details of equations and simulation processes are presented in Supplemental Experimental Procedures.
Modeling the Reversible Phosphorylation with Diffusion Process
The phosphorylation and de-phosphorylation processes were describedas the Michaelis-Menten scheme for enzymatic reaction, consisting of thebinding of enzyme and substrate, the dissociation of enzyme-substrate complex, and the reaction of phosphate transfer/removal. Quasi-steady-state
approximation on the formation of enzyme-substrate complex (i.e., Michaelis-Menten approximation) was applied for deterministic simulations.
Random and Biased Parameter Search
For random parameter search, each reaction rate constant (k), Michaelis constant (Km), and diffusion coefficient (D) were randomly chosen from the exponentially distributed parameters within the range of 1–1,000 min–1, 0.01–
1,000 mM, and 0.001–100 mm2 min–1 (1.67 3 10 5–1.67 mm2 s–1), respectively.
The obtained spatial and oscillatory parameter sets were then clusteredthrough Ward’s algorithm. The four-state model has a horizontal symmetry(e.g., exchanging the position of S01 and S10) and a point symmetry (e.g.,exchanging the role of kinase and phosphatase). A distance matrix between
the parameter sets was made to cancel those inherent symmetries of thereaction scheme.
For each biased parameter search reflecting one or more of the three design
motifs, randomly generated parameters that satisfied the following constraints
were used for PDE calculation: (1) cyclic reaction motif: k5 < k1 and k4 < k8, (2)binding motif: Max {Km4, Km5} < Min {Km1, Km2, Km3, Km6, Km7, Km8}, and (3)biased diffusion motif: Max {DS01, DS10} < Min {DS00, DS11}.
Representative Parameter Set
To analyze the mechanisms of spatial pattern formation from each cluster, a
representative parameter set was created to represent the tendency of thegiven cluster’s parameter distribution. For monotonously increasing or
decreasing histogram, its extreme value was adopted. For a histogram with
an apparent single peak value, the peak value was adopted. For the parameters that had no significant tendency, a value around its average was adopted.
Bifurcation Analysis
The representative parameter set (four-state model, spatial cluster 1) was used
for the bifurcation analysis shown in Figures 3 and S3. Each parameter wasvaried among the indicated range. After the system reaches steady state,peak and trough of S01 concentration and spatial width between two points
that give half the concentration of S01 peak were calculated.
Stochastic Simulation
The representative parameter sets for clusters 1 and 3 and the representative
parameter set for the three-state model were all converted to the full mass-action parameter sets before conducting SSA calculations. For the single-molecule tracking, 150,000 substrate molecules and only one ‘‘tagged’’ substratemolecule were simulated in a common environment but recorded separately.

KAKENHI grant number 25221004, H.R.U.), a Grant-in-Aid for Scientific
Research on Innovative Areas (JSPS KAKENHI grant number 23115006,
H.R.U.), the strategic programs for R&D (President’s discretionary fund) of
RIKEN (H.R.U.), an intramural Grant-in-Aid from the RIKEN QBiC (H.R.U.), a
grant from Takeda Science Foundation, a Grant-in-Aid for Young Scientists(B) (JSPS KAKENHI grant number 25830146, K.L.O.), and a Grant-in-Aid for
Challenging Exploratory Research (JSPS KAKENHI grant number 16K14653,
K.L.O.).
Received: October 1, 2016
Revised: January 22, 2017
Accepted: March 30, 2017
Published: April 25, 2017; corrected online: May 4, 2017
REFERENCES
Alam-Nazki, A., and Krishnan, J. (2012). An investigation of spatial signal transduction in cellular networks. BMC Syst. Biol. 6, 83.
Alam-Nazki, A., and Krishnan, J. (2013). Covalent modification cycles throughthe spatial prism. Biophys. J. 105, 1720–1731.
Alam-Nazki, A., and Krishnan, J. (2015). Spatial control of biochemical modification cascades and pathways. Biophys. J. 108, 2912–2924.
Alonso, S. (2016). Pattern formation at cellular membranes by phosphorylationand dephosphorylation of proteins. In Nonlinear Dynamics in Biological Systems, J. Carballido-Landeira and B. Escribano, eds. (Springer International
Publishing), pp. 63–82.
Alonso, S., and Bär, M. (2010). Phase separation and bistability in a threedimensional model for protein domain formation at biomembranes. Phys.
Biol. 7, 046012.
Alonso, S., and Bär, M. (2014). Modeling domain formation of MARCKS andprotein kinase C at cellular membranes. EPJ Nonlinear Biomedical Physics
2, 1–18.
Arrio-Dupont, M., Foucault, G., Vacher, M., Devaux, P.F., and Cribier, S.
(2000). Translational diffusion of globular proteins in the cytoplasm of cultured
muscle cells. Biophys. J. 78, 901–907.
Asai, R., Taguchi, E., Kume, Y., Saito, M., and Kondo, S. (1999). Zebrafish
leopard gene as a component of the putative reaction-diffusion system.
Mech. Dev. 89, 87–92.
Dawes, A.T., and Munro, E.M. (2011). PAR-3 oligomerization may provide an
actin-independent mechanism to maintain distinct par protein domains in
the early Caenorhabditis elegans embryo. Biophys. J. 101, 1412–1422.

SUPPLEMENTAL INFORMATION

Dudin, O., Merlini, L., and Martin, S.G. (2016). Spatial focalization of pheromone/MAPK signaling triggers commitment to cell-cell fusion. Genes Dev.
30, 2226–2239.

Supplemental Information includes Supplemental Experimental Procedures,
six figures, and four tables and can be found with this article online at http://dx.doi.org/10.1016/j.celrep.2017.03.081.

Economou, A.D., Ohazama, A., Porntaveetus, T., Sharpe, P.T., Kondo, S.,
Basson, M.A., Gritli-Linde, A., Cobourne, M.T., and Green, J.B. (2012). Periodic stripe formation by a Turing mechanism operating at growth zones in
the mammalian palate. Nat. Genet. 44, 348–351.

AUTHOR CONTRIBUTIONS

Etienne-Manneville, S., and Hall, A. (2002). Rho GTPases in cell biology.
Nature 420, 629–635.

S.S.S., K.L.O., and H.R.U. designed the study. S.S.S. performed most of thesimulation study. S.S.S., K.L.O., and H.R.U. analyzed the mechanism ofpattern formation and wrote the manuscript.

Ferrell, J.E., Jr., and Bhatt, R.R. (1997). Mechanistic studies of the dual phosphorylation of mitogen-activated protein kinase. J. Biol. Chem. 272, 19008–
19016.

ACKNOWLEDGMENTS
We thank all the lab members at RIKEN QBiC and the University of Tokyo, in
particular Dr. Craig C. Jolley for assisting the computer simulation and Dr.
Arthur Millius for critical reading of this manuscript. We are grateful to Drs. Shigeru Kondo, Shuji Ishihara, and Ryokichi Tanaka for valuable comments. This
work was supported by a grant from AMED-CREST (AMED/MEXT, H.R.U.),
CREST (JST/MEXT, H.R.U.), Brain/MINDS (AMED/MEXT, H.R.U.), Basic Science and Platform Technology Program for Innovative Biological Medicine(AMED/MEXT, H.R.U.), a Grant-in-Aid for Scientific Research (S) (JSPS

Flotow, H., Graves, P.R., Wang, A.Q., Fiol, C.J., Roeske, R.W., and Roach, P.J.
(1990). Phosphate groups as substrate determinants for casein kinase I action.
J. Biol. Chem. 265, 14264–14269.
Gierer, A., and Meinhardt, H. (1972). A theory of biological pattern formation.
Kybernetik 12, 30–39.
Gillespie, D.T. (1977). Exact stochastic simulation of coupled chemical reactions. J. Phys. Chem. 81, 2340–2361.
Goldbeter, A., and Koshland, D.E., Jr. (1981). An amplified sensitivity arisingfrom covalent modification in biological systems. Proc. Natl. Acad. Sci. USA
78, 6840–6844.

Cell Reports 19, 863–874, April 25, 2017 873

Goryachev, A.B., and Pokhilko, A.V. (2008). Dynamics of Cdc42 network embodies a Turing-type mechanism of yeast cell polarity. FEBS Lett. 582, 1437–
1443.
Halatek, J., and Frey, E. (2012). Highly canalized MinD transfer and MinE
sequestration explain the origin of robust MinCDE-protein dynamics. Cell
Rep. 1, 741–752.
Hoege, C., and Hyman, A.A. (2013). Principles of PAR polarity in Caenorhabditis elegans embryos. Nat. Rev. Mol. Cell Biol. 14, 315–322.
Ishihara, S., and Kaneko, K. (2006). Turing pattern with proportion preservation. J. Theor. Biol. 238, 683–693.
Jilkine, A., and Edelstein-Keshet, L. (2011). A comparison of mathematicalmodels for polarization of single eukaryotic cells in response to guided cues.
PLoS Comput. Biol. 7, e1001121.
Jolley, C.C., Ode, K.L., and Ueda, H.R. (2012). A design principle for a posttranslational biochemical oscillator. Cell Rep. 2, 938–950.

Meinhardt, H. (2008). Models of biological pattern formation: From elementary
steps to the organization of embryonic axes. Curr. Top. Dev. Biol. 81, 1–63.
Michaelis, L., Menten, M.L., Johnson, K.A., and Goody, R.S. (2011). The
original Michaelis constant: Translation of the 1913 Michaelis-Menten paper.
Biochemistry 50, 8264–8269.
Mori, Y., Jilkine, A., and Edelstein-Keshet, L. (2008). Wave-pinning and cell
polarity from a bistable reaction-diffusion system. Biophys. J. 94, 3684–3697.
€ller, P., Rogers, K.W., Jordan, B.M., Lee, J.S., Robson, D., Ramanathan, S.,
Mu
and Schier, A.F. (2012). Differential diffusivity of Nodal and Lefty underlies areaction-diffusion patterning system. Science 336, 721–724.
Otsuji, M., Ishihara, S., Co, C., Kaibuchi, K., Mochizuki, A., and Kuroda, S.
(2007). A mass conserved reaction-diffusion system captures properties ofcell polarity. PLoS Comput. Biol. 3, e108.

Karsenti, E. (2008). Self-organization in cell biology: A brief history. Nat. Rev.
Mol. Cell Biol. 9, 255–262.

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. USA 96, 4971–4976.

Kholodenko, B.N. (2006). Cell-signalling dynamics in time and space. Nat. Rev.
Mol. Cell Biol. 7, 165–176.

Seaton, D.D., and Krishnan, J. (2011). The coupling of pathways and processes through shared components. BMC Syst. Biol. 5, 103.

Kholodenko, B.N., Hoek, J.B., and Westerhoff, H.V. (2000). Why cytoplasmic
signalling proteins should be recruited to cell membranes. Trends Cell Biol.
10, 173–178.

Sheth, R., Marcon, L., Bastida, M.F., Junco, M., Quintana, L., Dahn, R., Kmita,
M., Sharpe, J., and Ros, M.A. (2012). Hox genes regulate digit patterning bycontrolling the wavelength of a Turing-type mechanism. Science 338, 1476–
1480.

Kondo, S., and Asal, R. (1995). A reaction-diffusion wave on the skin of the marine angelfish Pomacanthus. Nature 376, 765–768.
Kondo, S., and Miura, T. (2010). Reaction-diffusion model as a framework forunderstanding biological pattern formation. Science 329, 1616–1620.
Krishnan, J. (2009). Signal processing through a generalized module of adaptation and spatial sensing. J. Theor. Biol. 259, 31–43.
Loose, M., Fischer-Friedrich, E., Ries, J., Kruse, K., and Schwille, P. (2008).
Spatial regulators for bacterial cell division self-organize into surface waves
in vitro. Science 320, 789–792.

Thomson, M., and Gunawardena, J. (2009). Unlimited multistability in multisite
phosphorylation systems. Nature 460, 274–277.
Tostevin, F., and Howard, M. (2008). Modeling the establishment of PAR protein polarity in the one-cell C. elegans embryo. Biophys. J. 95, 4512–4522.
Turing, A.M. (1952). The chemical basis of morphogenesis. Philos. Trans. R.
Soc. Lond. B Biol. Sci. 237, 37–72.

Markevich, N.I., Hoek, J.B., and Kholodenko, B.N. (2004). Signaling switches
and bistability arising from multisite phosphorylation in protein kinase cascades. J. Cell Biol. 164, 353–359.

Vanselow, K., Vanselow, J.T., Westermark, P.O., Reischl, S., Maier, B., Korte,
T., Herrmann, A., Herzel, H., Schlosser, A., and Kramer, A. (2006). Differential
effects of PER2 phosphorylation: Molecular basis for the human familial
advanced sleep phase syndrome (FASPS). Genes Dev. 20, 2660–2672.

Markevich, N.I., Tsyganov, M.A., Hoek, J.B., and Kholodenko, B.N. (2006).
Long-range signaling by phosphoprotein waves arising from bistability in protein kinase cascades. Mol. Syst. Biol. 2, 61.

Zieske, K., and Schwille, P. (2013). Reconstitution of pole-to-pole oscillationsof min proteins in microengineered polydimethylsiloxane compartments.
Angew. Chem. Int. Ed. Engl. 52, 459–462.

874 Cell Reports 19, 863–874, April 25, 2017

Cell Reports, Volume 19

Supplemental Information

A Design Principle for an Autonomous
Post-translational Pattern Formation
Shuhei S. Sugai, Koji L. Ode, and Hiroki R. Ueda

A

B

Time (min)

0.5

500

0

1

0

2

4
6
8
Position (μm)

100Concentration (mM)Concentration of S01 (mM)

1000

1
0.8
0.6
0.4
0.2
0

S00
S01
S10
S11
0

500
Time (min)1000

Figure S1. Two-site reversible phosphorylation can generate oscillation, related to Figure 1.
(A) An example of oscillation, showing the concentration of S 01 along position and time. Vertical dotted line indicates the spatial
section at the position of 5 µm shown in (B).
(B) A section of (A) at the position of 5 µm shows the limit cycle oscillation. Concentration of each substrate is plotted along with
time.

A

Figure S2. Random parameter search for oscillatory
parameter sets, related to Figure 2.
(A) Oscillatory parameter sets found Figure 2A wereclustered with Ward’s algorithm.
(B, C) For each cluster of the oscillatory parameter sets,
parameter histograms and a pictorial representation for
tendencies of the parameters are depicted in the samemanner as in Figure 2. No significant trend is found in
the parameter distribution of diffusion coefficients.

0Distance

10000

Oscillatory parameter sets

1

B

2

ClusterC

Cluster 1

Cluster 2

F P
S01

S01

P

S00

S11

P

E

P F P
S11

S00
P

●
●
●
●

●
●
●
●

●
●
●
●

●
●
●
●

●
●
●

●
●

●
●
●

●
●
●
●

●
●
●
●

−5
−3
−1
1
Log10 D (μm2 sec-1)

Ds00
Ds01
Ds10
Ds11

0.0

Probability
0.2

0.4

S10

0.0

Probability
0.2

0.4

1
2
3
Log10 k (min-1)-2 -1 0 1 2
Log10 Km (μM)

30.2

Probability
0.25
0.0
0.0

Probability
0.1

0.2

-2 -1 0 1 2
Log10 Km (μM)

0Probability
0.1

3

Km1
Km2
Km3
Km4
Km5
Km6
Km7
Km8

Probability
0.1
0.2

3

0.0
1
2
Log10 k (min-1)

0.5

0

0.0

0.3

P E
S10
k1
k2
k3
k4
k5
k6
k7
k8

P
P

●
●
●

●
●
●
●

●
●

●
●

●
●
●

●
●
●

●
●
●

●
●
●
●

●
●
●

●

−5
−3
−1
1
Log10 D (μm2 sec-1)

k1
k2
k3
k4
k5
k6
k7
k8
Km1
Km2
Km3
Km4
Km5
Km6
Km7
Km8

Ds00
Ds01
Ds10
Ds11

1
2
Log10 k1 (min-1)

Km1

6
4
2
0

8

-2

-1

0
2
1
Log10 Km1 (μM)

Km5

6
4
2

8

-2

-1

1
0
2
Log10 Km5 (μM)

DS00

6
4
2
0

-2
-1
0
-3
Log10 DS00 (μm2 sec-1)012Log10 k2 (min-1)

2
-2

-1

1
0
2
Log10 Km2 (μM)

428S01 concentration (mM)

Km6

-2

-1

1
0
2
Log10 Km6 (μM)

3

0
-2
-1
-3
Log10 DS10 (μm2 sec-1)
k7

6
4
2

8

0

1
2
Log10 k7 (min-1)

Km3

4
2

8

-2

-1

0
2
1
Log10 Km3 (μM)

Km7

4
2
0

-2

-1

0
2
1
Log10 Km7 (μM)

3

k8

6
4
2

8

0

1
2
Log10 k8 (min-1)

3

-1

0
1
2
Log10 Km4 (μM)

3

-1

0
1
2
Log10 Km8 (μM)

3

Km4

6
4
2
0

3

6

8

0

3

6

0

3

6

8

0

Km2

4

8

0

3

6

0

3

S01 concentration (mM)

0

2

8

2

k2

4

0

3

3

6

0

3

S01 concentration (mM)

S01 concentration (mM)

8

0

1
2
Log10 k6 (min-1)

4

S01 concentration (mM)

2

8

0

6

S01 concentration (mM)

40

S01 concentration (mM)

k1

Oscillation

2
0

3

6

4

DS10

S01 concentration (mM)

8

1
2
Log10 k4 (min-1)

S01 concentration (mM)

S01 concentration (mM)

B

0

6

8

S01 concentration (mM)2

k6

S01 concentration (mM)

48

S01 concentration (mM)

60

S01 concentration (mM)k4
S01 concentration (mM)

8

S01 concentration (mM)

S01 concentration (mM)

A

8

-2
Km8

6
4
2
0

-2

DS11

6
4

0

2
0

5
10
Half-peak width (µm)

-2
-1
0
-3
Log10 DS11 (μm2 sec-1)

Figure S3. Bifurcation analysis to find the critical parameter for peak height and width, related to Figure 3.
The bifurcation analyses were conducted as shown in Figure 3E, and all the parameters except for the parameters shown in Figure
3E are shown in this figure. The symmetrical counterparts of parameters shown in Figure 3E are shown in (A), and the others areshown in (B).
See also Table S1.

Average distribution
throughout 100 min

Number of overall
production
100

1000

2

50

500

0

0

0

-2

-50

-500

-4

-100

-1000

Averaged count

1.2
0.8
0.4
0
Averaged count

1.2

B

5
0
-5
-10
10

0

0.4
0

10

5

0.8

0

Number of events per minute4

Number of events per minuteLog10 of averaged count

A

-5
2

4
6
Position (μm)

8

10

-10

0

2

4
6
8
Position (μm)

10Number of overall
rightward movements
S00
S01
S10
S11

30
20
10
0
-10
-20
-30
30
20
10
0
-10
-20
-30
0

ES00→01
ES00→10ES01→11
ES10→11E
FS01→00
FS10→00
FS11→01
FS11→10F2

468
Position (μm)

10

Log10 of averaged countDiffusion rates of enzymes and enzyme-substrate complexes
× 1/100

4
2

2

0

0

-2

-2

-4

0

2

4
6
8
Position (μm)

× 100

4

10

-4

S00
S01
S10
S11
0

2

4
6
Position (μm)

810Figure S4. Stochastic simulation, related to Figure 4.
(A) The stochastic simulation was conducted with a parameter set corresponding to the representative parameter set for cluster 1(see Table S2). For the final spatial pattern arising in the simulation, averaged distribution (left), overall production rate (middle),
and overall rightward movements (right) are shown for each bin position. Distribution was averaged for 400 minutes after theestablishment of spatial pattern. In the same time frame, overall production rate was evaluated by deducing depletion of thespecies of interest from its production, and overall movements were calculated for each bin by deducing leftward movements (fromthe next bin at right to the bin of interest) from the rightward movements (toward the next bin at right from the bin of interest).
The distribution of enzyme-substrate complexes is explicitly calculated in SSA simulation. As these complexes are converted into
products at a constant rate, their distributions are proportionate with the reaction rates calculated with PDE (see Figure 3D).
(B) Changing the diffusion rates of enzymes and enzyme-substrate complexes did not qualitatively alter the shape of the spatial
pattern. Of the parameter values shown in Table S2, the diffusion rates of enzymes and enzyme-substrate complexes were alteredeither to one hundred times (1.667 μm2 sec-1) or to one-hundredth (1.667×10-4 μm2 sec-1) of the original value (0.01667 μm2 sec-1).
In both cases, distribution of each species was averaged for 100 minutes after the establishment of spatial pattern.
See also Table S2.

A

B

DS11 ↓

Inner P F P
cycle S11

E

Outer
cycle

Y-position (μm)(k8 ↑)(k7 ↑)

S00 concentration

30

S00

20

20

10

10

0

0

S11
0

Y-position (μm)

C

0.8 × k80.82 × k8

0

1
2
3
Concentration (mM)

30

20

20

20

20

10

10

10

10

10
20
X-position (μm)

30

0

D

0

0

10
20
X-position (μm)

30

0

1.2
2.4
S11 Concentration (mM)

1/10 × k7

0

0

1/3 × k7

10
20
X-position (μm)

30

0

1.8
3.6
S11 Concentration (mM)

30

20

20

20

20

10

10

10

10

30

0

0
1.8
3.6
S11 Concentration (mM)

E

0

10
20
X-position (μm)

0

1.8
3.6
S11 Concentration (mM)

1/10 × DS11

30

0

0

0

1/3 × DS11

10
20
X-position (μm)

30

0

1.8
3.6
S11 Concentration (mM)

30

20

20

20

20

10

10

10

10

30

0
3
6
S11 Concentration (mM)

0

0

10
20
X-position (μm)

0

3
1.5
4.5
S11 Concentration (mM)

30

0

0

0

10
20
X-position (μm)

30

0

10
20
X-position (μm)

30

10 × DS11

30

10
20
X-position (μm)

1020X-position (μm)

0.4
0.8
S11 Concentration (mM)

3 × DS11

30

0

0

0
1.8
3.6
S11 Concentration (mM)

30

0

30

10 × k730

10
20
X-position (μm)

123Concentration (mM)

3 × k7

30

0

0

0

30

0

10
20
X-position (μm)

1.3 × k830

001.22 × k8

30

0
0.3
0.6
0.9
S11 Concentration (mM)Y-position (μm)

30

30

0

Y-position (μm)

1020X-position (μm)

S11 concentration

30

30

1.5
0.5
1
S11 Concentration (mM)

0

0

10
20
X-position (μm)

30

0
0.2
0.4
S11 Concentration (mM)

Figure S5. 2D pattern for the cluster 2-represenatative parameter set, related to Figure 6.
(A) Overall mechanism of spatial pattern formation and maintenance with cluster 2-representative parameter set similar to that withcluster 1-representative parameter set.
(B) 2D patterns emerged from cluster 2-representative parameter set.
(C-E) Modification of one of the parameter in the cluster 2-representative parameter set. Only the distributions of S 11 are shown.
See also Table S1.

A

Cluster 3-representativeB
30
Y-position (μm)

Y-position (μm)

30
20
10
0

100 × DS00

0

10
20
X-position (μm)

30

15
45
0
30
S00 Concentration (mM)

20
10
0

0

0

10
20
X-position (μm)

30

2
4
6
S00 Concentration (mM)

Figure S6. 2D pattern for the cluster 3-representative parameter set, related to Figure 7.
(A) Cluster 3-representative parameter set exhibits a small dot-shaped pattern, with most of the substance taking form of S 00 andaccumulating in the middle of the simulated area. The scale-independent emergence of single peak might be useful in circumstances where creating a single signal gradient in cell and subcellular structure is important, for example asymmetric cell division.
(B) Increasing DS00 caused broadening of the single peak, making a smooth gradient.
See also Table S1.

Table S1. Representative parameter sets of the four-component Michaelis-Menten model, related to
Figures 1, 3, 6, S3, S5, and S6.
Parameter

Cluster 1

Cluster 2

Cluster 3

1000

30

30

k2 (min-1)

30

30

30

k3 (min-1)

30

1000

300

k4 (min-1)

30

3

300

k5 (min-1)

30

3

1000

k6 (min-1)

30

1000

1000

k7 (min-1)

30

30

300

k8 (min-1)

1000

30

300

Km1 (μM)

0.3

0.01

20

Km2 (μM)

300

300

20

Km3 (μM)

300

0.3

0.02

Km4 (μM)

0.01

300

0.02

Km5 (μM)

0.01

300

30

Km6 (μM)

300

0.3

30

Km7 (μM)

300

300

1000

Km8 (μM)

0.3

0.01

1000

DS00 (μm2 sec-1)

1.667

5×10-4

1.667×10-5

DS01 (μm2 sec-1)

5×10-41.667

1.667×10-3

DS10 (μm2 sec-1)

5×10-41.667

1.667×10-3

DS11 (μm2 sec-1)

1.667

5×10-4

1.667-1

k1 (min )Table S2. Representative parameter sets and collapsed parameter sets of the four-component model,
related to Figures 4, 5, and S4.
Parameter
kb1 (min-1 μM-1)
kb2 (min-1 μM-1)
kb3 (min-1 μM-1)
kb4 (min-1 μM-1)
kb5 (min-1 μM-1)
kb6 (min-1 μM-1)
kb7 (min-1 μM-1)
kb8 (min-1 μM-1)
kub1-kub8 (min-1)
k1 (min-1)
k2- k7 (min-1)
k8 (min-1)
DS00 (μm2 sec-1)
DS01 (μm2 sec-1)
DS10 (μm2 sec-1)
DS11 (μm2 sec-1)
DES0001, DES0010,
DES0111, DES1011,
DFS0100, DFS1000,
DFS1101, DFS1110,
D E, D F(μm2 sec-1)

3367
0.13330.13334000
4000
0.13330.1333
3367
10
1000
30
1000
1.667
5×10-4
5×10-4
1.667

Reactionmotifcollapsed
274
0.274
0.274
8210
8210
0.274
0.274
274
10
72
72
72
1.667
5×10-4
5×10-41.667

0.01667

0.01667

Representative

Binding
Diffusion
motif
motif
collapsed
collapsed
249
3367
9.87
0.13339.87
0.1333
9.87
4000
9.87
4000
9.87
0.13339.87
0.1333
249
3367
10
10
1000
1000
30
30
1000
1000
1.667
0.02
5×10-4
0.02
-4
5×10
0.02
1.667
0.02

0.01667

0.01667

All motifs
collapsed
22
22
22
22
22
22
22
22
10
72
72
72
0.02
0.02
0.02
0.02

0.01667

Table S3. Representative parameter set of the three-component Michaelis-Menten model, related to
Figure 7.
Parameter

Representative-1

k1 (min )20

k2- k4 (min-1)

1000

Km1 (μM)

16

Km2 (μM)

0.01

Km3 (μM)

16

Km4 (μM)

103

DS0 (μm2 sec-1)

1.667×10-5

DS1 (μm2 sec-1)

1.667×10-3

DS2 (μm2 sec-1)

1.667

Table S4. Representative parameter set of the three-component model converted to the full mass-action
parameter, related to Figure 7.
Parameter

Representative-1

kb1 (min μM )1.875

kb2 (min-1 μM-1)

1.01×105kb3 (min-1 μM-1)

63.13

kb4 (min-1 μM-1)

1.01

kub1-kub4 (min-1)

10

k1 (min-1)

20

k2-k4 (min-1)

103

DS0 (μm2 sec-1)

1.667×10-5

DS1 (μm2 sec-1)

1.667×10-3

DS2 (μm2 sec-1)

1.667

DES0, DES1, DFS1, DFS2, DE, DF(μm2 sec-1)

0.01667-1

SUPPLEMENTAL EXPERIMENTAL PROCEDURES
Software for computer simulation
Numerical integration of systems with PDE was carried out by using Mathematica software version 9.0 and
10.0 (Wolfram Research). Clustering of parameter sets was carried out by R software (The R project). The
stochastic simulation was carried out by SSA (Gillespie, 1977) implemented in Cain software version 1.9(http://cain.sourceforge.net/).
Reaction-diffusion equations for the single-site reversible phosphorylation
The reaction-diffusion equations for the reversible mono-phosphorylation scheme are described as below.
𝜕

(1.1)

[𝑆 ],
𝜕𝑥 2 1𝜕[𝐸] = −𝑅ES0 + 𝐷E 2 [𝐸],
𝜕𝑡
𝜕𝑥

(1.2)

𝜕𝑡
𝜕
𝜕
𝜕𝑡
𝜕
𝜕𝑡
𝜕
𝜕𝑡

[𝑆0 ] = 𝑅S0 + 𝐷S0

𝜕

[𝑆0 ],𝜕𝑡𝜕

[𝑆1 ] = 𝑅S1 + 𝐷S1

𝜕𝑥 2
𝜕[𝐹] = −𝑅FS1 + 𝐷F𝜕𝜕𝑥 2

[𝐹𝑆1 ] = 𝑅FS1 + 𝐷FS1

[𝐹],
𝜕

[𝐸𝑆0 ] = 𝑅ES0 + 𝐷ES0

(1.3)

𝜕𝑥 2
𝜕
𝜕𝑥 2

(1.4)

[𝐸𝑆0 ],(1.5)[𝐹𝑆1 ],(1.6)

𝑅S0 = −𝑘𝑏1 [𝐸][𝑆0 ] + 𝑘𝑢𝑏1 [𝐸𝑆0 ] + 𝑘2 [𝐹𝑆1 ],(1.7)

𝑅S1 = −𝑘𝑏2 [𝐹][𝑆1 ] + 𝑘𝑢𝑏2 [𝐹𝑆1 ] + 𝑘1 [𝐸𝑆0 ],(1.8)

𝑅ES0 = 𝑘𝑏1 [𝐸][𝑆0 ] − 𝑘𝑢𝑏1 [𝐸𝑆0 ] − 𝑘1 [𝐸𝑆0 ],(1.9)

𝑅FS1 = 𝑘𝑏2 [𝐹][𝑆1 ] − 𝑘𝑢𝑏2 [𝐹𝑆1 ] − 𝑘2 [𝐹𝑆1 ].

(1.10)For any given substance S, Rs stands for the “reaction term” denoting the phosphorylation and dephosphorylation reactions. The “diffusion term,” denoting each species’ own free-diffusion process withthe diffusion coefficient Ds, is described as following;
𝐷S

𝜕
𝜕𝑥 2

[𝑆].

(1.11)[S0] and [S1] denote the concentrations of de- and mono-phosphorylated substrate. [E] and [F] denote theconcentrations of kinase and phosphatase. [ES0] and [FS1] denote the concentrations of the enzymesubstrate complexes. Mass-action kinetics parameters, kb, kub, and k, stand for the velocity coefficients ofsubstrate-enzyme binding, unbinding, and phospho-group transfer. Properties of phosphorylation reactions
are determined by kb1, kub1, and k1, while de-phosphorylation reactions are described by kb2, kub2, and k2.
Linear stability analysis
A reaction-diffusion system is said to have Turing instability or diffusion-induced instability when theuniform steady state is stable against a spatially homogeneous perturbation and unstable against aheterogeneous perturbation. To determine whether a given system can display Turing instability, linear
stability analysis was conducted on the Jacobian matrix at the uniform steady state. For n-species system,
the uniform steady state S* = ([S1*], [S2*],,, [Sn*]) was evaluated either by solving Rs = 0 for every chemical
species, or by running the simulation without the diffusion term until the system converges to S*. Supposethat the system resting at S* receives a small perturbation approximated as σm eikx (m = 1, 2,,, n), a sinefunction with wavenumber k ≥ 0 and a small amplitude σm in spatial position x. Then the reaction-diffusion
equations are reduced to;
𝟎 = (𝑹 − 𝑘 2 𝑫 − 𝜆𝑬)𝜮,

(2.1)

where
𝑹 = {𝑟𝑝𝑞 | 1 ≤ 𝑝 ≤ 𝑛, 1 ≤ 𝑞 ≤ 𝑛, 𝑟𝑝𝑞 =𝜕𝜕[𝑆𝑝 ]𝑅S𝑞 (𝑺∗ )},𝑫 = {𝑑𝑝𝑞 |1 ≤ 𝑝 ≤ 𝑛, 1 ≤ 𝑞 ≤ 𝑛, 𝑑𝑝𝑞 = 𝛿𝑝𝑞 𝐷s𝑝 }.

(2.2)
(2.3)

E is an identity matrix and Σ is a vector with n members denoting a fluctuation from the steady state. Let
λmax be one of the solutions of the equation below, with the largest value of real part|𝑹 − 𝑘 2 𝑫 − 𝜆𝑬| = 0.

(2.4)

For each k ≥ 0, Re (λmax) represents growth rate of a perturbation of sine-wave shape with wavenumber k.
Then the condition for Turing instability is as follows;
Re(𝜆max ) < 0 for 𝑘 = 0 and Re(𝜆max ) > 0 for some 𝑘 > 0.
(2.5)
Absence of Turing instability in the single-site system
For the single-site reversible phosphorylation scheme described as the equations (1.1)-(1.10), thecharacteristic equation (2.4) becomes a sixth-degree polynomial equation of λ. From Routh-Hurwitz
criterion (Edelstein-Keshet, 2005), this equation only yields solutions whose real part is negative, regardlessof the parameter set. From the condition (2.5), possibility of Turing instability in this system is ruled out.
Reaction-diffusion equations for the four-state model
The four-state model shown in Figure 1A is described as below;
𝜕
𝜕𝑡
𝜕
𝜕𝑡
𝜕
𝜕𝑡
𝜕
𝜕𝑡
𝜕
𝜕𝑡
𝜕
𝜕𝑡
𝜕
𝜕𝑡
𝜕
𝜕𝑡
𝜕
𝜕𝑡
𝜕
𝜕𝑡
𝜕
𝜕𝑡
𝜕
𝜕𝑡
𝜕
𝜕𝑡
𝜕
𝜕𝑡

[𝑆00 ] = 𝑅S00 + 𝐷S00
[𝑆01 ] = 𝑅S01 + 𝐷S01
[𝑆10 ] = 𝑅S10 + 𝐷S10
[𝑆11 ] = 𝑅S11 + 𝐷S11

𝜕
𝜕𝑥 2
𝜕
𝜕𝑥 2
𝜕
𝜕𝑥 2
𝜕
𝜕𝑥 2

[𝑆00 ],(3.1)[𝑆01 ],(3.2)[𝑆10 ],(3.3)[𝑆11 ],(3.4)[𝐸] = −𝑅ES0001 −𝑅ES0010 −𝑅ES0111 −𝑅ES1011 + 𝐷E
[𝐹] = −𝑅FS0100 −𝑅FS1000−𝑅FS1101 −𝑅FS1110 + 𝐷F
[𝐸𝑆0001 ] = 𝑅ES0001 + 𝐷ES0001[𝐸𝑆0010 ] = 𝑅ES0010 + 𝐷ES0010[𝐸𝑆0111 ] = 𝑅ES0111 + 𝐷ES0111[𝐸𝑆1011 ] = 𝑅ES1011 + 𝐷ES1011[𝐹𝑆0100 ] = 𝑅FS0100 + 𝐷FS0100[𝐹𝑆1000 ] = 𝑅FS1000 + 𝐷FS1000[𝐹𝑆1101 ] = 𝑅FS1101 + 𝐷FS1101[𝐹𝑆1110 ] = 𝑅FS1110 + 𝐷FS1110

𝜕
𝜕𝑥 2
𝜕
𝜕𝑥 2
𝜕
𝜕𝑥 2
𝜕
𝜕𝑥 2
𝜕
𝜕𝑥 2
𝜕
𝜕𝑥 2
𝜕
𝜕𝑥 2
𝜕
𝜕𝑥 2

𝜕
𝜕𝑥 2
𝜕
𝜕𝑥 2

[𝐸],

(3.5)

[𝐹],

(3.6)

[𝐸𝑆0001 ],(3.7)[𝐸𝑆0010 ],(3.8)[𝐸𝑆0111 ],(3.9)[𝐸𝑆1011 ],(3.10)[𝐹𝑆0100 ],(3.11)[𝐹𝑆1000 ],(3.12)[𝐹𝑆1101 ],(3.13)[𝐹𝑆1110 ],(3.14)

𝑅S00 = −(𝑘𝑏1 + 𝑘𝑏2 )[𝐸][𝑆00 ] + 𝑘𝑢𝑏1 [𝐸𝑆0001 ] + 𝑘𝑢𝑏2 [𝐸𝑆0010 ] + 𝑘5 [𝐹𝑆0100 ] + 𝑘6 [𝐹𝑆1000 ],(3.15)
𝑅S01 = −𝑘𝑏3 [𝐸][𝑆01 ] − 𝑘𝑏5 [𝐹][𝑆01 ] + 𝑘𝑢𝑏3 [𝐸𝑆0111 ] + 𝑘𝑢𝑏5 [𝐹𝑆0100 ] + 𝑘1 [𝐸𝑆0001 ] + 𝑘7 [𝐹𝑆1101 ],(3.16)
𝑅S10 = −𝑘𝑏4 [𝐸][𝑆10 ] − 𝑘𝑏6 [𝐹][𝑆10 ] + 𝑘𝑢𝑏4 [𝐸𝑆1011 ] + 𝑘𝑢𝑏6 [𝐹𝑆1000 ] + 𝑘2 [𝐸𝑆0010 ] + 𝑘8 [𝐹𝑆1110 ],(3.17)
𝑅S11 = −(𝑘𝑏7 + 𝑘𝑏8 )[𝐹][𝑆11 ] + 𝑘𝑢𝑏7 [𝐹𝑆1101 ] + 𝑘𝑢𝑏8 [𝐹𝑆1110 ] + 𝑘3 [𝐸𝑆0111 ] + 𝑘4 [𝐸𝑆1011 ],(3.18)
𝑅ES0001 = 𝑘𝑏1 [𝐸][𝑆00 ] − 𝑘𝑢𝑏1 [𝐸𝑆0001 ] − 𝑘1 [𝐸𝑆0001 ],(3.19)

𝑅ES0010 = 𝑘𝑏2 [𝐸][𝑆00 ] − 𝑘𝑢𝑏2 [𝐸𝑆0010 ] − 𝑘2 [𝐸𝑆0010 ],(3.20)

𝑅ES0111 = 𝑘𝑏3 [𝐸][𝑆01 ] − 𝑘𝑢𝑏3 [𝐸𝑆0111 ] − 𝑘3 [𝐸𝑆0111 ],(3.21)

𝑅ES1011 = 𝑘𝑏4 [𝐸][𝑆10 ] − 𝑘𝑢𝑏4 [𝐸𝑆1011 ] − 𝑘4 [𝐸𝑆1011 ],(3.22)

𝑅FS0100 = 𝑘𝑏5 [𝐹][𝑆01 ] − 𝑘𝑢𝑏5 [𝐹𝑆0100 ] − 𝑘5 [𝐹𝑆0100 ],(3.23)

𝑅FS1000 = 𝑘𝑏6 [𝐹][𝑆10 ] − 𝑘𝑢𝑏6 [𝐹𝑆1000 ] − 𝑘6 [𝐹𝑆1000 ],(3.24)

𝑅FS1101 = 𝑘𝑏7 [𝐹][𝑆11 ] − 𝑘𝑢𝑏7 [𝐹𝑆1101 ] − 𝑘7 [𝐹𝑆1101 ],(3.25)

𝑅FS1110 = 𝑘𝑏8 [𝐹][𝑆11 ] − 𝑘𝑢𝑏8 [𝐹𝑆1110 ] − 𝑘8 [𝐹𝑆1110 ].

(3.26)For the massive random parameter search, Michaelis-Menten approximation was applied toreduce the number of variables. It assumes the enzyme-substrate binding and unbinding as having reachedequilibrium, and reduces the set of reaction rate parameters {kb, kub, k} into {Km, k}. Km is called Michaelis
constant given by𝐾𝑚 =𝑘𝑢𝑏 +𝑘𝑘𝑏.

(3.27)In the Michaelis-Menten -approximated version, the entire system is described by four equations,exchanging the reaction terms of (3.1)-(3.4) to the following;
𝑅S00 = − (𝑅S01 = (𝑅S10 = (𝑅S11 = (
[𝐸] =

𝑘1
𝐾𝑚1

𝑘1

𝐾𝑚1
𝑘2
𝐾𝑚2
𝑘3
𝐾𝑚3

𝑘2
𝐾𝑚2

[𝑆00 ] −[𝑆00 ] −[𝑆01 ] +[𝐸tot ]

) [𝐸][𝑆00 ] + (
𝑘3

𝐾𝑚3
𝑘4
𝐾𝑚4
𝑘4
𝐾𝑚4

𝑘5
𝐾𝑚5

[𝑆01 ] +[𝑆01 ]) [𝐸] + (−[𝑆10 ]) [𝐸] + (−[𝑆10 ]) [𝐸] − (

𝑘5
𝐾𝑚5
𝑘6
𝐾𝑚6

𝑘7

𝐾𝑚7

+

𝑘6
𝐾𝑚6

[𝑆10 ]) [𝐹],[𝑆01 ] +[𝑆10 ] +
𝑘8
𝐾𝑚8

𝑘7
𝐾𝑚7
𝑘8
𝐾𝑚8

(3.28)[𝑆11 ]) [𝐹],(3.29)[𝑆11 ]) [𝐹],(3.30)) [𝐹][𝑆11 ],(3.31),(3.32)[𝐹tot ][𝑆 ] [𝑆 ] [𝑆 ] [𝑆 ].
1+ 01 + 10 + 11 + 11

(3.33)[𝑆 ] [𝑆 ] [𝑆 ] [𝑆 ]1+ 00 + 00 + 01 + 10
𝐾𝑚1

[𝐹] =

+

𝐾𝑚5

𝐾𝑚2

𝐾𝑚6

𝐾𝑚3

𝐾𝑚7

𝐾𝑚4

𝐾𝑚8

In our model, the total amount of enzymes are [Etot] = [Ftot] = 1 μM and total amount of substrate was 1,500
μM unless otherwise indicated.
Random parameter search in the four-state model
The random parameter search in this paper aimed to find Turing instability-based spatial pattern formationand spatially homogeneous oscillation from reversible two-site phosphorylation. The 1D system size was
set to 10 µm with a periodic boundary condition. Spatial pattern is defined by three criteria; 1) Turing
instability, 2) temporal stability (i.e. temporal change of any species’ concentration throughout a certain
time is sufficiently small), and 3) spatial heterogeneity (at least one species shows a significant variance inits concentration distribution). Oscillation is characterized by two criteria; 4) temporal cyclicity, and 5)
spatial homogeneity (i.e. spatial variance of any species’ concentration throughout the simulated space issufficiently small) of the chemical concentration.
For the four-state model, the following two steps were conducted to reduce calculation time.
First, for each parameter set, a system without D was simulated under ordinary differential equations(ODEs). In the initial state, all the substrate exists as non-phosphorylated form (S00) with the concentrationof 1,500 µM. Kinase and phosphatase exist at the concentration of 1 µM. All the other substances have zero
initial values. For each parameter set, calculation was continued at least until t = 1,000 min and at mostuntil t = 3 × 108 min. Parameter sets that showed either oscillatory behaviors or Turing instability were thenallowed to proceed to the second step. The other parameter sets, converging to Turing-stable uniform steadystate, or failing to converge within the calculation limit, were discarded. The concrete criteria of oscillation
and Turing instability are as follows. For oscillation, the temporal trajectory of any of the simulated species
must cross the same value three times. The two parts of trajectories sectioned by those three crossing points
must display a high similarity with correlation coefficient greater than 0.99, and the peak amplitude of thetwo trajectories must not differ by greater than 1%. For Turing instability, the system must have converged
to a steady state, and the determinant equation of Jacobian matrix at the steady state must have at least one
solution with positive real part (see Linear stability analysis). This satisfies the criterion 1).
The second step uses the PDE with both reaction and diffusion term, and calculates the evolutionof the system both along time and space. The initial condition of the PDE simulation was set as the final
state of the ODE, hence mostly the same as the steady state of ODE. 0.1% fluctuation was added to theconcentration of S00 for each spatial position. To prevent the memory from expiring, each simulation was
split into several epochs with duration T. The first epoch has T = 2,500 min. Each epoch was extended totwice the duration of the previous epoch if the previous epoch was calculated without errors (e.g. probably
the parameter set was near a bifurcation point and the behavior of the system was too stiff). If there was anerror in the previous epoch, the same epoch was re-calculated with the duration shortened to half theduration of the previous epoch. If the error did not resolve even when the epoch size reaches T = 100 min
or shorter, then the simulation was discarded.
The presence of stable spatial pattern is judged according to the following steps. As criterion 1)

is already established in the first step, the second step aims to establish criteria 2) and 3). After calculatingeach epoch, the concentration difference of two time points, T and 0.99T, were evaluated at each position
for each species. If the maximum of the concentration difference was less than 0.01 μM, then the simulation
result was judged as having reached the equilibrium and met criterion 2). Then, for each of S00, S01, S10,
and S11, the difference between its maximum and minimum concentration was calculated. If the largest oneof the four values, named Cmax, satisfies Cmax > 100 μM then it was judged as having met criterion 3), hence
the presence of stable spatial pattern. For oscillation, the time evolution at an arbitrary position in the space(in our simulation, x = 2.5 μm) must meet the same criteria as in the first step and be therefore considered
as having met criterion 4). The system must also meet criterion 5) with Cmax < 0.01 μM.
Clustering with collapsed symmetry
We used Ward’s algorithm when clustering the spatial and oscillatory parameter sets. The distance matrices
for the clustering were made with a metric to cancel out the inherent symmetries of the reaction scheme.
The four-state model has a horizontal symmetry (e.g. exchanging the position of S 01 and S10 conserves theoverall structure) and a point-symmetry (e.g. exchanging the role of kinase and phosphatase conserves theoverall structure). The distance between the two parameter sets P 1 and P2 is defined as the minimum of the
distances between P1 and the four symmetrical images of P2 – Min {P1P2, P1P3, P1P4, P1P5} – where P3 andP4 are the horizontal and point-symmetrical transpositions, and P5 is the combination of the twotranspositions. P3 is defined by transposing the parameters of P 2 as k1 ↔ k2, k3 ↔ k4, k5 ↔ k6, k7 ↔ k8, Km1↔ Km2, Km3 ↔ Km4, Km5 ↔ Km6, Km7 ↔ Km8, and DS01 ↔ DS10. P4 is defined by converting the parametersof P2 as k1 ↔ k8, k2 ↔ k7, k3 ↔ k6, k4 ↔ k5, Km1 ↔ Km8, Km2 ↔ Km7, Km3 ↔ Km6, Km4 ↔ Km5, DS00 ↔ DS11,
and DS01 ↔ DS10. The combination of two transpositions, P 5, is defined by the product of the above two
steps. The order of the two transpositions does not affect the result, both being same as k1 ↔ k7, k2 ↔ k8, k3↔ k5, k4 ↔ k6, Km1 ↔ Km7, Km2 ↔ Km8, Km3 ↔ Km5, Km4 ↔ Km6, and DS00 ↔ DS11.
Stochastic simulation
In Gillespie algorithm (Gillespie, 1977), two random numbers are generated for each iteration step; the firstis to determine when the next process takes place, and the second is to determine what the process will be.
Each random number is generated from a set of probability density functions that reflect the number ofpossible next processes and their ratio. In this study, the 1D spatial range (0-10 μM) with a periodic
boundary condition was split into 100 identical bins. It is assumed that each bin is homogeneous andexchanges substances with the neighboring bins with a defined set of diffusion coefficients. By assuming
diffusion as conversion process of a species in one bin into a species in the neighboring bins, diffusion canbe treated in the same way as reaction.
For the simulations, the representative parameter sets in Michaelis-Menten-approximatedversion were converted to full mass-action parameter sets by fixing kub = 10 min-1 and solving the equation(3.27). The set of two parameters {Km, k} was thus converted to the three parameters {kub, kb, k} for eachreaction. For the diffusion term, enzymes and enzyme-substrate complexes were all given the value of 1
µm2min-1 (0.01667 µm2sec-1). See Table S2 for the values.
Analysis of spatial pattern formation in the three-state model
When we consider kinase/phosphatase actions in two-site phosphorylation, we can choose from random orordered models. The random model corresponds to the four-state model, while the ordered model
corresponds to the three-state model. Although the ordered three-state model applies a strict constraint on
the action of enzymes, it is mathematically simpler than the random model. On the other hand, the random
four-state model is a biochemically more relaxed form and thus a more general model of two-site reversible
modification. The three-state model with progressive two-site phosphorylation and the shared kinase andphosphatase is described as the equations below;
𝜕
𝜕𝑡
𝜕
𝜕𝑡
𝜕
𝜕𝑡
𝜕
𝜕𝑡
𝜕
𝜕𝑡
𝜕
𝜕𝑡
𝜕
𝜕𝑡
𝜕
𝜕𝑡
𝜕
𝜕𝑡

[𝑆0 ] = 𝑅S0 + 𝐷S0[𝑆1 ] = 𝑅S1 + 𝐷S1[𝑆2 ] = 𝑅S2 + 𝐷S2

𝜕
𝜕𝑥 2
𝜕
𝜕𝑥 2
𝜕
𝜕𝑥 2

[𝑆0 ],(4.1)[𝑆1 ],(4.2)[𝑆2 ],[𝐸] = −𝑅ES0 −𝑅ES1 + 𝐷E

[𝐸],

(4.4)

[𝐹],

(4.5)

[𝐸𝑆0 ],(4.6)[𝐸𝑆1 ],(4.7)[𝐹𝑆1 ],(4.8)[𝐹𝑆2 ],(4.9)[𝐹] = −𝑅FS1 −𝑅FS2 + 𝐷F
[𝐸𝑆0 ] = 𝑅ES0 + 𝐷ES0
[𝐸𝑆1 ] = 𝑅ES1 + 𝐷ES1
[𝐹𝑆1 ] = 𝑅FS1 + 𝐷FS1
[𝐹𝑆2 ] = 𝑅FS2 + 𝐷FS2

𝜕
𝜕𝑥 2
𝜕
𝜕𝑥 2
𝜕
𝜕𝑥 2
𝜕
𝜕𝑥 2

(4.3)
𝜕

𝜕𝑥 2
𝜕
𝜕𝑥 2

𝑅S0 = −𝑘𝑏1 [𝐸][𝑆0 ] + 𝑘𝑢𝑏1 [𝐸𝑆0 ] + 𝑘3 [𝐹𝑆1 ],(4.10)

𝑅S1 = −𝑘𝑏3 [𝐹][𝑆1 ] + 𝑘𝑢𝑏3 [𝐹𝑆1 ] + 𝑘1 [𝐸𝑆0 ]−𝑘𝑏2 [𝐸][𝑆1 ] + 𝑘𝑢𝑏2 [𝐸𝑆1 ] + 𝑘4 [𝐹𝑆2 ],(4.11)

𝑅S2 = −𝑘𝑏4 [𝐹][𝑆2 ] + 𝑘𝑢𝑏4 [𝐹𝑆2 ] + 𝑘2 [𝐸𝑆1 ],(4.12)

𝑅ES0 = 𝑘𝑏1 [𝐸][𝑆0 ] − 𝑘𝑢𝑏1 [𝐸𝑆0 ] − 𝑘1 [𝐸𝑆0 ],(4.13)

𝑅ES1 = 𝑘𝑏2 [𝐸][𝑆1 ] − 𝑘𝑢𝑏2 [𝐸𝑆1 ] − 𝑘2 [𝐸𝑆1 ],(4.14)

𝑅FS1 = 𝑘𝑏3 [𝐹][𝑆1 ] − 𝑘𝑢𝑏3 [𝐹𝑆1 ] − 𝑘3 [𝐹𝑆1 ],(4.15)

𝑅FS2 = 𝑘𝑏4 [𝐹][𝑆2 ] − 𝑘𝑢𝑏4 [𝐹𝑆2 ] − 𝑘4 [𝐹𝑆2 ].

(4.16)The names of the species and the indices of reaction parameters in the equations correspond to the ones in
Figure 7A. The above set of nine PDEs was used in the stochastic simulation. With Michaelis-Menten
approximation, the entire system is written by three variables for S0, S1, and S2, exchanging the reaction
terms of (4.1)-(4.3) to the following;𝑘1

𝑅S0 = −

𝐾𝑚1
𝑘1

𝑅S1 = (
𝑅S2 =

𝐾𝑚1
𝑘2

𝐾𝑚2

[𝐸][𝑆0 ] +[𝑆0 ] −

𝑘2

𝑘3
𝐾𝑚3

[𝐹][𝑆1 ],[𝑆 ]) [𝐸] + (

𝐾𝑚2 1
𝑘4

[𝐸][𝑆1 ] −𝐾𝑚4(4.17)𝑘4
𝐾𝑚4[𝑆2 ] −𝑘3
𝐾𝑚3

[𝑆1 ]) [𝐹],(4.18)[𝐹][𝑆2 ],(4.19)

where
[𝐸] =
[𝐹] =

[𝐸tot ]
[𝑆0 ] [𝑆1 ]+𝐾𝑚1 𝐾𝑚2,(4.20)[𝐹tot ]
[𝑆1 ] [𝑆2 ]+𝐾𝑚3 𝐾𝑚4.

(4.21)

1+

1+

In the random parameter search with the three-state model, the same range for reaction rate
constants (k), Michaelis constants (Km), and diffusion coefficients (D) were used as with the four-statemodel. In the initial state, all the substrate exists as non-phosphorylated form (S0) with the concentration of
1,500 µM. The kinase and phosphatase exist at the concentration of 1 µM. All the other substances havezero initial values. Each simulation for a given parameter set was composed of a sequence of epochs withthe duration of T = 10,000 min. The calculation was continued at least for 1 epoch and at most for 20 epochs.
If the result of simulation of an epoch met the following two conditions, it was determined that thesimulation reached a stable spatial pattern; first, standard deviation of spatial distribution of at least one
species’ concentration must exceed 100 μM at the end of the epoch. Second, standard deviation of thetemporal change of every species’ concentration throughout the epoch must not exceed 0.1 μM at two points(x = 0 and x = 0.23).
Absence of Turing pattern in the two-site reversible phosphorylation model without enzyme sharing
With Michaelis-Menten approximation, it was shown that Turing pattern could never occur from the fourstate scheme without enzyme sharing (i.e. all the reactions are carried out by different enzymes). The
reaction-diffusion equations of the four-state model without enzyme sharing are described as below insteadof (3.28) through (3.33);
𝑅S00 = −𝑘1
𝐾𝑚1 +[𝑆00 ][𝐸tot1 ][𝑆00 ] −𝑘2
𝐾𝑚2 +[𝑆00 ][𝐸tot2 ][𝑆00 ] +𝑘5
𝐾𝑚5 +[𝑆01 ][𝐹tot1 ][𝑆01 ] +𝑘6
𝐾𝑚6 +[𝑆10 ][𝐹tot2 ][𝑆10 ],(5.1)

𝑅S01 =

𝑘1
𝐾𝑚1 +[𝑆00 ][𝐸tot1 ][𝑆00 ] −𝑘3
𝐾𝑚3 +[𝑆01 ][𝐸tot3 ][𝑆01 ] −𝑘5
𝐾𝑚5 +[𝑆01 ][𝐹tot1 ][𝑆01 ] +𝑘7
𝐾𝑚7 +[𝑆01 ][𝐹tot3 ][𝑆11 ],(5.2)

𝑅S10 =𝑘2
𝐾𝑚2 +[𝑆00 ][𝐸tot2 ][𝑆00 ] −𝑘4
𝐾𝑚4 +[𝑆10 ][𝐸tot4 ][𝑆10 ] −𝑘6
𝐾𝑚6 +[𝑆10 ][𝐹tot2 ][𝑆10 ] +𝑘8
𝐾𝑚8 +[𝑆10 ][𝐹tot4 ][𝑆11 ],(5.3)

𝑅S11 =𝑘3
𝐾𝑚3 +[𝑆01 ][𝐸tot3 ][𝑆01 ] +𝑘4
𝐾𝑚4 +[𝑆10 ][𝐸tot4 ][𝑆10 ] −𝑘7
𝐾𝑚7 +[𝑆11 ][𝐹tot3 ][𝑆11 ] −𝑘8
𝐾𝑚8 +[𝑆11 ][𝐹tot4 ][𝑆11 ].

(5.4)

For this set of equations, possibility of Turing instability was algebraically ruled out with any combination

of reaction parameters and total amounts of enzymes; all the Routh-Hurwitz coefficients are positive withno change of signs. Possibility of Turing instability in the three-state model without enzyme sharing was
also ruled out by setting k2 = k4 = k6 = k8 = 0. For the three-state model, we attempted the same analysis for
all possible patterns of enzyme sharing, using Routh-Hurwitz criterion. It was shown that Turing pattern ispossible when two phosphorylating reactions share the kinase, or when two de-phosphorylating reactions
share the phosphatase, or when the both are satisfied at the same time; in all the other cases, spatial pattern
emergence was ruled out. The results of analyses in Michaelis-Menten-approximated schemes are notdirectly applicable to the full mass-action scheme, especially where the fundamental premise for theapproximation fails (e.g. where the total amount of enzymes and substrates are at the similar order). We did
not conduct the same analysis on the full mass-action scheme, as it became unfeasible within our limit ofcomputation power.
Algebraic condition for Turing instability in the three-state model with shared kinase andphosphatase
In the case of three-state reversible phosphorylation scheme with a shared kinase and a shared phosphatase,the value of each rpq in definition (2.2) are as follows;
𝑟11 = −

𝑘1 𝐾𝑚1 𝐾𝑚2 [𝐸tot ](𝐾𝑚2 +[𝑆1 ∗ ])2

(𝐾𝑚1 𝐾𝑚2 +𝐾𝑚2 [𝑆0 ∗ ]+𝐾𝑚1 [𝑆1 ∗ ])

𝑘 𝐾

𝐾

[𝐹

𝑟12 = 3 𝑚3 𝑚4 tot

(6.1)

< 0,
2

](𝐾𝑚2 [𝑆0 ∗ ]+𝐾𝑚1 (𝐾𝑚2 +[𝑆1 ∗ ])) (𝐾𝑚4 +[𝑆2 ∗ ])+𝑘1 𝐾𝑚1 𝐾𝑚2 [𝐸tot ][𝑆0 ∗ ](𝐾𝑚4 [𝑆1 ∗ ]+𝐾𝑚3 (𝐾𝑚4 +[𝑆2 ∗ ]))2∗

(𝐾𝑚2 [𝑆0 ]+𝐾𝑚1 (𝐾𝑚2 +[𝑆1 ∗ ]))2 (𝐾𝑚4 [𝑆1 ∗ ]+𝐾𝑚3 (𝐾𝑚4 +[𝑆2 ∗ ]))2

> 0,
(6.2)

𝑟13 = −
𝑟31 = −

𝑘3 𝐾𝑚3 𝐾𝑚4 [Ftot ][𝑆1 ∗ ](𝐾𝑚4 [𝑆1 ∗ ]+𝐾𝑚3 (𝐾𝑚4 +[𝑆2 ∗ ]))

2

< 0,

(6.3)

2

< 0,

(6.4)

𝑘2 𝐾𝑚1 𝐾𝑚2 [𝐸tot ][𝑆1 ∗ ](𝐾𝑚2 [𝑆0 ∗ ]+𝐾𝑚1 (𝐾𝑚2 +[𝑆1 ∗ ]))

𝑘 𝐾

𝐾

[𝐹

𝑟32 = 4 𝑚3 𝑚4 tot][𝑆2 ∗ ](𝐾𝑚2 [𝑆0 ∗ ]+𝐾𝑚1 (𝐾𝑚2 +[𝑆1 ∗ ]))2 +𝑘2 𝐾𝑚1 𝐾𝑚2 [𝐸tot ](𝐾𝑚1 +[𝑆0 ∗ ])(𝐾𝑚4 [𝑆1 ∗ ]+𝐾𝑚3 (𝐾𝑚4 +[𝑆2 ∗ ]))2
(𝐾𝑚2 [𝑆0 ∗ ]+𝐾𝑚1 (𝐾𝑚2 +[𝑆1 ∗ ]))2 (𝐾𝑚4 [𝑆1 ∗ ]+𝐾𝑚3 (𝐾𝑚4 +[𝑆2 ∗ ]))2

> 0,
(6.5)

𝑟33 = −

𝑘4 𝐾𝑚3 𝐾𝑚4 [𝐹tot ](𝐾𝑚3 +[𝑆1 ∗ ])(𝐾𝑚4 [𝑆1 ∗ ]+𝐾𝑚3 (𝐾𝑚4 +[𝑆2 ∗ ]))

2

< 0.

(6.6)

The Jacobian matrix can be rewritten using positive numbers,
𝐹11 = −𝑟11 , 𝐹12 = 𝑟12 , 𝐹13 = −𝑟13 , 𝐹31 = −𝑟31 , 𝐹32 = 𝑟32 , 𝐹33 = −𝑟33 ,

(6.7)

Then,
𝑹 − 𝑘 2 𝑫 − 𝜆𝑬 =−𝐹11 − 𝐷S0 ∗ 𝑘 2 − 𝜆(𝐹11 + 𝐹31
−𝐹31
Then (2.4) becomes:

𝐹12−𝐹12 − 𝐹32 − 𝐷S1 ∗ 𝑘 2 − 𝜆
𝐹32

−𝐹13
)
𝐹33 + 𝐹13
−𝐹33 − 𝐷S2 ∗ 𝑘 2 − 𝜆𝑎0 𝜆3 + 𝑎1 𝜆2 + 𝑎2 𝜆 + 𝑎3 = 0,

(6.8)

(6.9)

where
𝑎0 = −1 < 0,

(6.10)

𝑎1 = 𝑘 2 (−𝐷S0 − 𝐷S1 − 𝐷S2 ) − 𝐹11 − 𝐹12 − 𝐹32 − 𝐹33 < 0,

(6.11)

𝑎2 = 𝑘 4 (−𝐷S0 𝐷S1 − 𝐷S1 𝐷S2 − 𝐷S0 𝐷S2 ) + 𝑘 2 (−𝐷S2 (𝐹11 + 𝐹12 + 𝐹32 ) − 𝐷S1 (𝐹11 + 𝐹33 ) − 𝐷S0 (𝐹12 +𝐹32 + 𝐹33 )) + 𝐹12 𝐹31 + 𝐹13 𝐹31 + 𝐹13 𝐹32 − 𝐹11 𝐹32 − 𝐹11 𝐹33 − 𝐹12 𝐹33,
(6.12)
𝑎3 =𝑘 2 (−𝑘 4 𝐷S0 𝐷S1 𝐷S2 − 𝑘 2 (𝐷S0 𝐷S2 (𝐹12 + 𝐹32 ) + 𝐷S1 (𝐷S2 𝐹11 + 𝐷S0 𝐹33 )) + 𝐷S0 (𝐹13 𝐹32 − 𝐹12 𝐹33 ) +𝐷S1 (𝐹13 𝐹31 − 𝐹11 𝐹33 ) + 𝐷S2 (𝐹12 𝐹31 − 𝐹11 𝐹32 )).
(6.13)As the first step for solving the analytic condition (2.5) for the system to have Turing instability,
the condition for stability at k = 0 is investigated. Since a3 becomes zero for k = 0, equation (6.9) becomes𝜆(𝑎0 𝜆2 + 𝑎1 𝜆 + 𝑎2 ) = 0,

(6.14)

where λ = 0 is a trivial solution arising from mass-conservation. Its corresponding eigenvector points to thedirection that all the substance increases at the same time, breaching the mass-conservation itself. This isbiochemically meaningless and is omitted from consideration. The equation (6.14) is then reduced to;𝑎0 𝜆2 + 𝑎1 𝜆 + 𝑎2 = 0.

(6.15)For both of the two solutions of this equation to have negative real parts, a2 < 0 (k = 0) must be met. Thisis equivalent to;
𝐹12 𝐹31 + 𝐹13 𝐹31 + 𝐹13 𝐹32 − 𝐹11 𝐹32 − 𝐹11 𝐹33 − 𝐹12 𝐹33 < 0.

(6.16)

Next, the condition for instability for some k > 0 is investigated. From Routh-Hurwitz criterion,
the signs of the 4 Routh-Hurwitz coefficients must not be the same. These coefficients are a0, a1, a3 and b1,
which is written as𝑎 𝑎𝑏1 = 𝑎2 − 0 3.
𝑎1

(6.17)

Substitution of equation (6.17) by k, DS0, DS1, DS2, F11, F12, F13, F31, F32, and F33 yields a polynomial with
only negative coefficients. This proves that b1 is always negative. Combined with equations (6.10) and(6.11), indicating that both a0 and a1 are negative, Turing instability can be made possible only when a3 >0. We investigated the condition to satisfy
𝑎3 ′ = 𝑎3 /𝑘 2 = −𝑘 4 𝐷S0 𝐷S1 𝐷S2 − 𝑘 2 (𝐷S0 𝐷S2 (𝐹12 + 𝐹32 ) + 𝐷S1 (𝐷S2 𝐹11 + 𝐷S0 𝐹33 )) + 𝐷S0 (𝐹13 𝐹32 −𝐹12 𝐹33 ) + 𝐷S1 (𝐹13 𝐹31 − 𝐹11 𝐹33 ) + 𝐷S2 (𝐹12 𝐹31 − 𝐹11 𝐹32 ) > 0.
(6.18)
The coefficients of k2 and k4 are both negative, making a3’ monotonically decrease. Therefore, the only case
we have to consider in the periodic boundary condition for the spatial range is k = 2π.
From the above analysis, a given parameter set must satisfy all the conditions in the following
steps in order to result in a spatial pattern formation, all of which can be checked algebraically;(1) Evaluate S* = ([S0*], [S1*], [S2*]) by solving the zero point of reaction terms of reaction-diffusion
equations. In the three-state model with Michaelis-Menten approximation, this yields a cubic equationand thus algebraically solvable. From S*, extract the solutions that are within a chemically possible
range (i.e. real and positive). This yields either one or three different solutions.
(2) Ensure that S* is a stable steady state without diffusion. This can be done by solving |R - λE| = 0 andchecking that Re (λmax) < 0. If S* is unstable, it must be omitted from the later analysis.
(3) For k = 2π, if any of the following two tautologous conditions are met, then the given parameter set
would result in a spatial pattern: i) a3 > 0, or ii) Re (λmax) > 0 where λ is given by solving a cubic
equation |R – k2D - λE| = 0.
A detailed analysis on the pattern formation in two-dimensional space for cluster 2 and 3
Though the cluster 1 and 2 are highly similar with 90-degree rotation of the reaction scheme (see Figure
2C and 2D) and similar regulation of pattern shape in 2D space was observed with the representative
parameter set for cluster 2, they are not exactly the same. For example, the peak height is resilient against
the changes in one parameter responsible for the “inner” cycle in cluster 2 (i.e. k7). The difference can be
attributed to the different parameter topology between cluster 1 and 2 because these clusters are not
perfectly symmetric: for example, in the cluster 1, high affinity between S 01 and phosphatase promotes thereaction S01 → S00 that counteracts the overall clockwise bias in the reaction rate. On the other hand, in thecluster 2, high affinity between S00 and kinase promoters the reaction from S00 → S01 that supports theoverall clockwise reaction bias.
In the following explanation, the parameters of cluster 2 are compared with those of cluster 1with clockwise 90-degree rotation – for example, S00 of cluster 1 is compared with S01 of cluster 2 and k1of cluster 1 is compared with k3 of cluster 2. With the representative parameter set of cluster 1, the onset ofthe “inner cycle” (i.e. S01 → S00) is mediated by phosphatase, while the onset of the “outer cycle” (i.e. S 01
→ S11) is mediated by kinase. S01 can thus proceed either to the “inner” cycle or the “outer” cycle at asimilar order of reaction rate. On the other hand, with cluster 2-representative parameter set, the two
possible reactions for S11 (i.e. S11 → S01 and S11 → S10) are both catalyzed by phosphatase. The directionin which S11 proceeds is determined by the ratio of Km7 and Km8. As Km8 is overwhelmingly smaller thanKm7, S11 predominantly proceeds to the “outer cycle.” Thus, k8 in cluster 2, the first component of the “outer”
cycle, plays a much more crucial role to the maintenance of spatial pattern than k3 in cluster 1 does. A small
alteration of k8 significantly changes the pattern width; only 18% decrease results in an S11-dominant pattern
and 22% increase produces an S00-dominant pattern (Figure S5C). The heavy dependence to the “outer”
cycle makes the pattern collapse easily with only 20% decrease or 30% increase (Figure S5C). On the other
hand, k7 in cluster 2, the first component of the “inner” cycle, does not contribute so much to themaintenance of spatial pattern formation as k5 does in cluster 1, thus altering k7 has little effect on the patternshape (Figure S5D). As diffusion process is independent from the chemical reaction, alteration of DS11
produced mostly the same result as altering DS01 in cluster 1 (Figure S5E); decrease of DS11 resulted inhigher peaks with a narrower pitch, while its increase produced an opposite effect with lower peaks with awider pitch.
With the representative parameter of cluster 3, almost all the substances accumulate at the single
position in the form of S00. Increasing DS00, the diffusion rate of the accumulating species, resulted in a
wider peak with a lower peak value (Figure S6A and B).

SUPPLEMENTAL REFERENCE
Edelstein-Keshet, L. (2005). Mathematical models in biology, Classics edn (Philadelphia: Society for
Industrial and Applied Mathematics).