Spontaneous emergence of catalytic cycles with colloidal spheres

Edited by Alexei V. Tkachenko, Brookhaven National Laboratory, Upton, NY, and accepted by Editorial Board Member John D. Weeks March 10, 2017 (received for review July 20, 2016)
April 10, 2017
114 (17) 4342-4347


Biological systems inspire a new paradigm for material synthesis, aiming to design materials that emulate living systems, providing both functionality and self-replication. Colloidal particles endowed with specific interactions provide a particularly promising approach for realizing this vision. Here, we consider a set of colloidal spheres with specific, time-dependent interactions and demonstrate that the interactions can be designed so that clusters of particles create other clusters through templating reactions. Surprisingly, simple templating rules generically give rise to the production of a sea of clusters of various sizes. The sea of clusters grows exponentially in a catalytic cycle. This is a specific realization of Dyson’s notion of an exponentially growing metabolism, emerging naturally from a simple scheme.


Colloidal particles endowed with specific time-dependent interactions are a promising route for realizing artificial materials that have the properties of living ones. Previous work has demonstrated how this system can give rise to self-replication. Here, we introduce the process of colloidal catalysis, in which clusters of particles catalyze the creation of other clusters through templating reactions. Surprisingly, we find that simple templating rules generically lead to the production of huge numbers of clusters. The templating reactions among this sea of clusters give rise to an exponentially growing catalytic cycle, a specific realization of Dyson’s notion of an exponentially growing metabolism. We demonstrate this behavior with a fixed set of interactions between particles chosen to allow a catalysis of a specific six-particle cluster from a specific seven-particle cluster, yet giving rise to the catalytic production of a sea of clusters of sizes between 2 and 11 particles. The fact that an exponentially growing cycle emerges naturally from such a simple scheme demonstrates that the emergence of exponentially growing metabolisms could be simpler than previously imagined.
The origin of life is usually associated with self-replication, the ability of an entity to create copies of itself (1). Thus inspired, there have been many efforts in recent years aimed at creating artificial systems with self-replicative capacity. These are typically inspired by the linear chain mechanism of DNA. However, living systems are more than individual entities that can replicate themselves. Dyson (2) and Oparin (3) argued that a more critical aspect of living systems is the creation of a metabolism, a complex cascade of chemical reactions that together are able to accomplish more than any single chemical reaction can do on its own. Using a simple mathematical model, the so-called “garbage bag model,” Dyson presented a scenario through which such a complex metabolism could arise spontaneously. His idea is that a random set of catalysts will catalyze arbitrary chemical reactions in a nonsynergistic fashion. However, if each catalyst is more likely to function when there are others that are synergistic with it, then there is a critical amount of cooperativity above which metabolic cycles spontaneously emerge. This idea has been explored through abstract simulations (4). In a similar spirit, Kauffman and coworkers have shown that if the probability of one species catalyzing the formation of another is above a threshold, then catalytic cycles naturally emerge (57).
Whether this scenario can naturally occur in practice remains obscure. The fundamental question is to determine how likely it is for a soup of interacting catalysts to self-organize into an exponentially growing catalytic cycle. How special do the interactions between the different components need to be for spontaneous exponentially growing catalytic cycles to emerge? In this work, we present an explicit demonstration of a set of interactions in which exponentially growing catalytic cycles naturally emerge. This set of interactions models a realizable physical system, colloidal spheres interacting with time-dependent specific interactions, in which particle clusters can template each other. The emergence of exponentially growing catalytic cycles depends on particle valences in the templating process.
We consider colloidal particles having stickers distributed uniformly over their surface. These stickers cause short-ranged, specific interactions between particles, with a programmable binding strength between pairs of particles. We have recently shown that this system has the potential of allowing high-yield equilibrium self-assembly of large structures, up to 1,000 particles (8). Classical experiments showed that DNA-based specific interactions allow assembly of crystals of nanoparticles (912). More recently, efforts for self-assembly at the colloidal scale have begun to bear fruit (1218). To achieve the exponentially growing catalytic cycles of the current study, we find that it is necessary to use time-dependent interactions between the particles, whose strength can be modulated over a programmable timescale. Schemes for time-dependent interactions based on DNA hybridization exist (19), and there is already an experimental realization of a form of time-dependent interactions (20), with more under development. We will see in what follows that the timescale of the interactions plays an essential role in determining the nature of the catalytic cycles that are produced in the reactions.


We define a “cluster” as a rigid structure in which any continuous deformation costs energy and in which N particles interact through specific short-range interactions. On the one hand, a cluster is characterized by its geometry: It has been shown that for N<6, for each N, there is only one rigid structure, and for N6, there are multiple rigid structures (2 for N=6, 5 for N=7, 13 for N=8, and so on) (21, 22), with nearly 106 at N=14 (23) (Fig. 1A). On the other hand, a cluster is characterized by the species of particles that comprise it, where the species of two particles determines the specific interaction between them. The number of different species of particles is not necessarily the same as N. It has been previously established how specific interactions can be chosen for any given geometry (8, 16, 24, 25) and how self-assembly yield of clusters depends on the chosen interactions (8, 24). In the present work, we do not consider self-assembly of clusters in a bath of particles (8, 16, 25, 26); however, our ability to enumerate cluster geometries and choose interactions to stabilize any desired cluster is a critical ingredient for the discoveries reported here. Instead, here we focus exclusively on formation of a cluster through templating from another cluster, which is a form of catalysis. Template-based catalysis is common in biology (27, 28) and offers a different route to increasing the production rate of a structure in exchange for energy expenditure.
Fig. 1.
(A) Landscape of different rigid structures that can be formed out of N particles. (B) Templating scheme for making an octahedron (III) from an N0=7 cluster (I). A catalyst particle can bind a monomer of complementary species, which are differently colored. (C1) Interaction matrix for the N0 = 7 cluster. Blue entries denote attraction between particle species. The set of complementary species has the exact same interaction rules. (C2) Contact matrix for the octahedron. Blue entries denote particles in contact. (C3) How the templating in B is possible: The contacts in the octahedron (C2) have to be allowed by the interactions of particle species in the N0 = 7 cluster (C1). This can only be realized (shaded red with particle assignments on bottom) if one of the particles (species 7) has valence 2.

Templating Catalysis.

To define our model system, we select one specific cluster of size N=N0, termed here the “parent catalyst,” and choose from the earlier enumerations (8, 24) the particle species and interactions that make it the only ground state of N0 particles. In this work, however, the particle number is not fixed, and instead of the self-assembly process, we consider exclusively templating catalytic reactions as means of creating clusters. The choice of interactions, therefore, only guarantees that the parent catalyst is indeed a rigid configuration of particles. The templating process, as will be defined shortly, is based on monomers attaching to clusters, which implies that particle valence is important in catalysis. The choice of allowed valences is, in principle, independent of the already introduced interaction rules.
We will present in detail the catalytic system arising when the parent catalyst is a specific N0=7 cluster (Fig. 1 B and C1), and then show that the presented behavior is generic in our model, occurring in a multitude of parent catalysts with sizes N0=6 to N0=8. To define the templating process and illustrate the significance of particle valence, we proceed by describing a specific templating catalysis reaction in which the parent catalyst templates an octahedron cluster (N=6).
A templating reaction occurs in three distinct steps, with our example illustrated in Fig. 1B: First, monomers from the bath bind to the surface of the catalyst cluster, with the bound monomers forming bonds with each other. In the second step, the network of bound monomers dissociates from the original cluster. In the third step, the dissociated network folds into a new cluster. To catalyze formation of a specific cluster, it is necessary to specify the rules for binding monomers to a given particle of the parent catalyst (i.e., specify the particle valence). Our previous work (29) showed that efficient autocatalysis of the octahedron, with at most one monomer binding to a parent particle (valence 1), required including an additional cluster in the reaction. Here, we consider templating from a single catalyst, without limiting the particle valence. To catalyze the formation of the octahedron from the N0=7 parent catalyst, it turns out that we must adopt the following rule—although particle species labeled 16 can only bind one monomer from the bath, particle species 7 can bind two monomers from the bath (Fig. 1C).
The reason for this can be seen in our algorithm for finding possible catalysts for a given cluster: In Fig. 1 C1, the “interaction matrix” shows which particle species in the catalyst can contact each other through their specific interactions. In Fig. 1 C2, the matrix shows which particles within the octahedron are in contact. For the N0=7 parent cluster to catalyze the octahedron, we must find a submatrix of the interaction matrix that matches the contacts in the octahedron. Fig. 1 C3 shows that this can only happen if particle species 7 has valence 2, while the other particles have valence 1. Using this algorithm, we searched through all geometries and compositions in our landscape (8, 21, 22, 24) and found multiple parent clusters that can efficiently catalyze the octahedron: The discussed parent catalyst is the only one with N0=7, and multiple parents with N0>7 can carry out the reaction efficiently (SI Text and Fig. S1). There is nothing special about the reaction in which the octahedron is the target; however, the choice of having particle valence larger than one strongly impacts the totality of catalytic reactions for a given parent catalyst, as we will discuss below.
Fig. S1.
(A–D) Parent clusters with N0=7 particles that can catalyze formation of the octahedron. For each catalyst we show: Cluster with numbered particles and colored particle species (Upper Left); the interaction matrix between particle species (Upper Right); and different ways the octahedron is catalyzed: numbers on particles denote the valence toward attaching monomers (zero means the particle is not involved in catalysis), and the ordered list above shows the numbered particle species (in Upper Left) that attach monomers that become the octahedron particles (Lower).
The reactions of template-driven catalysis require time-dependent bond strengths: The bond strengths of the monomers binding to the catalytic cluster must weaken in time, whereas the bonds between bound monomers must strengthen in time. The dissociation timescale of the network of bound monomers is set by the timescale of the weakening bonds. By varying this timescale, it is possible to change the size of the dissociated network of monomers and hence the size of the catalyzed cluster. This process is in spirit similar to the time-dependent type of enzyme catalysis, which has been established as a common type of catalysis in biology (28). In our simulations, we implement this varying timescale by triggering the dissociation event based on the number of bonds within the monomer network. We choose the necessary number of bonds for triggering dissociation from a distribution, where moving the distribution peak toward, e.g., larger numbers corresponds to a longer dissociation timescale (SI Text and Fig. S2).
Fig. S2.
Folding pathways of the octahedron, starting from all structures formed out of 6 particles (top row). The particle colors represent their interactions without detailing their possible species. Particles colored differently attract, whereas particles that repel are colored the same. Rows are labeled by the total number of particle bonds. Blue lines mark possible foldings of structures that can only end in the fully formed octahedron as the number of bonds increases. Red lines mark possible foldings that can end up in the local minimum state of the octahedron. Black disc color (legend) marks if monomers with nc=8 number of mutual contacts can form the given structure as they dissociate from the parent catalyst cluster with N0=7 particles. There are four possibilities, and for each there are folding pathways leading to the fully formed octahedron. The first such network from the left appears the most in the simulations and always leads to the ground state.
When the interactions are realized by DNA stickers, the monomer–catalyst binding involves complementary strands, so that for each particle species in the parent, there is a complementary particle species that can bind to it. The complementary species interact among themselves in the same way as the original species do. The catalytic system of the N0=7 parent catalyst (Fig. 1 C1) is therefore defined to have double the number of species, i.e., 14 species (Fig. 2A). We now proceed to test in detail the performance of the N0=7 parent catalyst, using dissipative particle dynamics (DPD) techniques introduced previously for colloidal particles with stickers (8, 29). Our simulation contains colloidal spheres of diameter d, with an interaction range of 1.05d. The particles are immersed into a DPD solvent of smaller particles. We simulated systems with 512 colloidal particles, out of which 7 comprise the parent catalyst. Simulations are run at a fixed temperature with a volume fraction of colloids ϕcoll=0.01. Fig. 2B shows a particular event of the parent cluster producing an octahedra in a simulation.
Fig. 2.
(A) The interaction matrix shows that the mutual interactions among the catalyst cluster particles (the first seven rows/columns) are identical to the ones among their complementary species (the last seven rows/columns). Each particle species can bind its complementary species (off-diagonal seven by seven blocks). Orange entries mark the only interactions with valence 2. (B) Simulation of catalysis of octahedra from the N0=7 parent catalyst. Particles in the catalyst are colored according to their interaction matrix, with the particles bound from the bath colored white (unbound particles and solvent particles are omitted for clarity). (B1–B4) Snapshots of an octahedron being formed. (B, Inset) Plot of the number of octahedra as function of time obtained in a simulation. (C) Snapshot of simulation started with a single N0=7 catalyst (labeled 1), where all generated clusters were allowed to template other clusters. A few selected clusters (labeled 2–9) are marked, with the sequence of them getting templated shown on the right. Only clusters 3 and 5 are nonrigid structures. Clusters 3 and 6 have eight particles, cluster 5 has nine particles, and the rest have seven particles.

Emergent Catalytic Reactions.

Although Fig. 2B demonstrates that the N0=7 parent catalyst indeed effectively templates octahedra, a close examination of the simulations demonstrates that the cluster is also capable of templating other clusters. In turn, the clusters that are templated from the parent N0=7 catalyst can themselves template even more clusters. The rich templating capability becomes apparent as we increase the number of monomer bonds necessary for dissociation event and thereby favor larger catalysis products. Fig. 2C shows a snapshot from a simulation, highlighting a few examples of clusters and nonrigid structures.
The root cause of the degenerate structures that are produced is the valence-2 particle in the N0=7 parent catalyst. Without this valence-2 particle, only some clusters with < 7 particles can be templated (SI Text and Fig. S3). However, the degeneracy causes an enormous zoo of clusters to be produced by the templating reactions. To enumerate the set of possible rigid clusters that can be produced starting from the N0=7 cluster, we take our database of cluster geometries up to 14 particles (2123) and ask which of these geometries can be realized by using the interaction matrix of the parent N0=7 catalyst (Fig. 1 C1), where multiple copies of the valence-2 particle are allowed (particle species 7). There are 100 different geometries that are consistent with this interaction matrix, from N=2 to N=11 particles. These give rise to a total of 176 different clusters, because some of the geometries have multiple colorings (i.e., some geometries can be constructed using different combinations of particle species). Of the 176 clusters, 135 have >5 particles, and we focus on these.
Fig. S3.
Directed graph showing the templating connections between all clusters realized by using the interaction matrix of the N0=7 parent cluster, where the valence of each particle species is 1. In contrast to Fig. 4, there are no cycles in this graph.
We next ask which of the 135 clusters are capable of templating each other. Note that from here on, for simplicity of presentation, we consider clusters with complementary colorings to be one and the same; i.e., a pair of complementary particle species is regarded as a single species. This simplification does not affect any of our results and conclusions (SI Text). Fig. 3 shows a directed graph, where each of the 135 nodes corresponds to a rigid cluster. A (directed) edge between two nodes indicates that the initial cluster can template the target one. There are 8,057 pairs of connected clusters in the graph, so the probability of a connection between any two clusters is 90%. This percentage is higher than the critical value predicted by Kauffman above which catalytic cycles emerge from a reactive graph (5, 6). Indeed, there are many cycles within our graph, an example of which is shown in Fig. 4A.
Fig. 3.
Directed graph showing the templating connections between all clusters with >5 particles, realized by using the interaction matrix of the N0=7 parent catalyst. There are 135 clusters in total.
Fig. 4.
(A) An example of emergent catalytic cycle that includes the N0=7 parent catalyst. (B) The reorganization of the directed graph in Fig. 3, where clusters are grouped by their number of particles (N) and by the species of each of their particles. A directed edge from a group of clusters (a circle) to another group means that all clusters in the former can template all clusters in the latter group. When the edge has an endpoint on a horizontal line that connects several circles, it represents multiple edges to each of those circles. The group to which the N0=7 cluster belongs is marked by a red circle.
To make the structure of the cycles visible, Fig. 4B reorganizes the graph in Fig. 3, grouping clusters by their size and coloring. For example, to generate a cycle with a N=11 cluster starting from the parent N0=7 catalyst, we need to template a type of N=8 cluster (8 possibilities out of 24 clusters), which then templates a type of N=9 cluster (31 possibilities out of 37 clusters). There are >1.7×106,1.0×107,6.2×108 cycles of length 3, 4, 5, respectively. The parent N0=7 catalyst is involved in 15% of the average number of cycles per cluster. The largest cycle we found has size 83 and involves the parent catalyst.

Exponential Growth of Catalytic Cycles.

Given this web of templating possibilities between the different clusters (Figs. 3 and 4B), what is the outcome when a N0=7 catalyst is placed in a monomer bath? To answer this, let 𝐜(t) be a vector of probabilities, where ci(t) is the normalized concentration of the ith cluster. Then,
where 𝐌 is the matrix of first-order production rates between the different clusters, with Mij representing the rate at which the jth cluster templates the ith cluster.
Now, Mij can only be nonzero if the corresponding edge in the directed graph (Figs. 3 and 4B) is nonzero. The magnitude of Mij is set by kinetic considerations: A templating reaction requires that all of the templated particles diffuse to the surface of the catalyst and then dissociate together once enough bonds between the particles are formed. The critical parameter is the timescale of dissociation, τdis, the length of time the templating particles are allowed to stay bound to the catalyst before they dissociate. The number of particles that template before dissociation increases with increasing τdis, leading to a preference for larger clusters. In contrast, if τdis is short, templating larger clusters is impossible, and the catalysts only produce small clusters.
A reasonable model chooses τdis so that templated clusters of a given size are favored, and then both larger and smaller clusters are disfavored. We model this by choosing the value of each matrix element of 𝐌, i.e., the value of kinetic rate for each templating reaction, which is inversely proportional to the corresponding τdis. Given a matrix 𝐌, what is the dynamics of Eq. 1, assuming that we start from a single instance of the N0=7 catalyst? If 𝐜 has K components, then Eq. 1 has K eigenmodes, namely, solutions Ciα that grow/shrink exponentially when Re[λα] is positive/negative, where λα is the corresponding eigenvalue. When we choose the kinetic rates to vary weakly (SI Text), the spectrum corresponding to our model is shown in Fig. 5A. We sort the eigenvalues so that Re[λ1] is the smallest and Re[λK] is the largest. The same figure shows the spectrum for a different choice of kinetic rates, drawn from an exponential distribution, where kinetic rates for production of the largest clusters happen to be relatively increased (SI Text). In both cases, the λK is real and λKRe[λα] for α<K, and we find this to be a robust feature under variation of kinetic rates. This is consistent with the Perron–Frobenius theorem, which guarantees a maximal real eigenvalue for a matrix representing a strongly connected graph (30). In our example, the 83 clusters from the largest cycle form the largest strongly connected component and lead to the maximal real eigenvalue. At long times, the solution will be dominated by the eigenmode corresponding to the largest exponential term,
given that the projection, β, of the initial condition onto this eigenmode is nonzero. We indeed find that for the initial condition with only the N0=7 cluster in the bath, there is robustly a significant overlap with the 𝐂K eigenmode.
Fig. 5.
(A) The spectrum corresponding to the kinetic equations for templating, with eigenvalues λα, α=1135. The two spectra belong to two different choices of kinetic rates for templating: slightly varying rates (squares), and rates chosen from an exponential distribution (circles) (Exponential Growth of Catalytic Cycles). The fastest growing eigenmode for each choice of kinetic rates corresponds to the marked eigenvalue. (B1) The fastest growing mode for slightly varying kinetic rates of templating. The size of the cluster image is proportional to the concentration of the cluster. Directed connections show the templating relations. (B2) The fastest growing mode for kinetic rates chosen from an exponential distribution (Exponential Growth of Catalytic Cycles). The size of the cluster image is proportional to the concentration of the cluster. We show the mutual catalytic relations only for five clusters.
The structure of the eigenmode 𝐂K is sensitive to the choice of kinetic rates, even though the existence of a dominating λK is robust. The clusters involved in the longtime solution are shown in Fig. 5B for our two different choices of kinetic rates. When all kinetic rates are similar, the solution involves clusters of all sizes, preferring slightly the smaller ones because these can be templated by both small and large clusters, whereas the converse does not hold. When kinetic rates are uneven, as in our second choice, only a subset of clusters is involved in the long-term solution, representing a smaller set of catalytic cycles, as shown in Fig. 5 B2.

Generality of Emergent Cycles in Colloidal Catalysis.

The question arises of how general is the behavior of the catalytic system based on the N0=7 parent catalyst? To answer this, we analyzed the catalytic systems arising from each cluster between N0=6 and N0=8 as the parent catalyst. Allowing a single (any) particle of the fixed parent (a total of 105 parents) to have valence 2 results in 716 different catalytic systems. We generate the directed graph for each, just as in Fig. 3, checking for geometries from N=6 to N=10 (a total of 333 geometries). Compared with the described N0=7 parent, we find that dozens of other parent catalysts generate catalytic systems that are at least as large and have as big a maximal cycle (Fig. 6, blue bars). We note that throughout the entire histogram (Fig. 6, blue bars), one can find catalytic systems based on any number of particle species; the minimum is 4 species (when the parent has 2), and the maximum is 16.
Fig. 6.
Histograms for the 716 different valence-2 catalytic systems (blue) and 105 different valence-1 catalytic systems (gray), according to their largest catalytic cycle. The number above the bar counts the number of systems having the maximal cycle size within the range of the bar. Each catalytic system is generated by parent cluster of size N0=6,7 or 8, with a single (any) particle species allowed valence 2 (blue) or all having valence 1 (gray). In each system we are considering only different geometries of clusters with N between 6 and 10. Systems without cycles >1 (i.e., the only cycles are self-catalysis) were omitted from the histogram (125 and 47 for valence 2 and 1, respectively). The arrow marks the bin in which lies the N0=7 catalytic parent discussed in detail.
In contrast, in the case when all particles have valence 1, the catalytic systems generated in the same way as above are drastically smaller; there are very few cycles, and all of them contain <6 clusters (Fig. 6, gray bars). Because of valence constraint, a cluster cannot template a larger one, so that cycles can occur only among clusters with the same N. Our example of N0=7 parent catalyst illustrates this: With valence set to 1, we find no catalytic cycles (SI Text), in stark contrast to the case with valence 2.


Our exponentially growing catalytic cycles appear to be an explicit realization of Dyson’s classical mechanism (2) for the emergence of catalytic cycles. Instead of a single proliferating entity, our system results in exponentially growing catalytic cycles, which is reminiscent of a metabolism. Were the reactions to occur in a cell volume that itself divided in response to, e.g. the amount of one of the clusters in our cycle, we would have a natural mechanism for an exponentially growing catalytic cell cycle.
Whereas previous models realizing Dyson’s metabolism are abstract (46), without specific physical realization, the present construction demonstrates that it is surprisingly easy for a templated set of reactions to give rise to an exponentially growing catalytic cycle. Our only caveat is that the specific nature and properties of the catalytic cycle are emergent properties of the programmed interactions and cannot easily be directly designed. The only special aspect of the system studied here is that we have a complete characterization of the energy landscape of the colloidal particles, and therefore are able to easily make the design choices for a templated catalytic reaction to emerge. However, we see no reason why similar phenomena could not emerge from other soft-matter systems, in particular, polymers or gels (31). As but one example, it is interesting to contemplate whether a system of RNA strands might be designed to lead to an exponential catalytic cycle.
We note that the basic ingredients required to experimentally realize our system are becoming available. Controlled valence of isotropic mesoscale particles has been demonstrated (18, 32). The interaction among monomers needs to strengthen in time, so that they do not aggregate in solution, but do bind to each other during time they are attached to the template. We believe that the dissociation from the template could be induced by the very process of bonding among attached monomers (33) and/or by weakening in time of the monomer–template bonds (19), both of which require consumption of energy. The weakening and breaking of monomer–template bonds might be controlled by an externally imposed (i.e., day–night) cycle for these bonds (31, 34, 35), which could set the disassociation timescale. A first step toward time-dependent interactions has been realized between nanoparticles by using complex strand-displacement reactions that rely on a DNA fuel source (20, 36).
SI Text

Possible Catalysts for the Octahedron

We present the analysis of templating of an octahedron from all clusters of size N0=7 from our database (Fig. S1), which is representative of our analysis of all N0=7,8 possible catalysts. The first step is the algorithm (described in Templating Catalysis) of identifying the contact matrix of the octahedron within the interaction matrix of the potential catalyst, which reveals the 4 possible catalysts out of the total 5 clusters with N0=7. All except one of the catalysts must have at least one particle with valence 2. Having completed this step, we are guaranteed that 6 monomers attached to a particular set of particles on any one of these catalysts (Fig. S1) will mutually interact in a way that octahedron is the only rigid structure they can form.
The next step in considering the catalysis of an octahedron is the question of geometry: In which way can the monomers attached to a catalyst geometrically bond with each other (we assume exactly 6 monomers attached to the proper cluster particles)? We find by direct inspection that on every catalyst in Fig. S1 the attached monomers can geometrically form a continuous network, i.e., a network that upon dissociation from the catalyst is in a single piece. All of the possible 6-monomer networks are listed in Fig. S2. The single catalyst in which all particles have valence 1 (Fig. S1D) is at a disadvantage: The monomers can only connect into a single chain, and that only when all reach a special spatial configuration, making the chain unlikely to form.
In the final step, we consider the folding of the dissociated monomer network into an octahedron (Fig. S2). For a successful folding, it is important to note that the octahedron ground state with 12 particle bonds has one excited state, a local minimum (LM) with two fewer bonds. The portion of monomer networks that can only fold into the octahedron, and not into the LM, grows as the number of bonds in the network is increased. On the other hand, smaller monomer networks assemble more quickly on the surface of the catalyst. (There are, however, exceptions such as the chain in previous paragraph.) Large monomer networks, e.g. with 9 mutual bonds, have to conform in shape to the surface of the catalyst, which additionally slows or prevents them from being templated.
We observe that for the N0=7 parent catalyst (Fig. S1A), the maximal number of bonds that a templated monomer network can have is nc=8. There are only 4 out of 22 different possible monomer networks having 8 bonds (Fig. S2). Of these four, one forms with overwhelming probability compared with the other three, and this particular monomer network can only fold into the octahedron ground state (Fig. S2). We note that candidate Fig. S1C can efficiently catalyze an octahedron, but has a trivially small catalytic network.

Dissociation Criterion in Simulations

In our DPD simulations, the dissociation of a monomer network from the catalyst happens when monomers have a certain number, nc, of bonds among themselves. As described in SI Text, Possible Catalysts for the Octahedron, when the N0=7 parent cluster catalyzes an octahedron, the value nc=8 is maximal and gives the highest chance of folding into an octahedron. Therefore, in our simulations that demonstrate the success of this particular catalysis process (Fig. 2), we fix the value nc=8. We confirmed that with other values, such as nc=6,7, the parent catalyst does catalyze the octahedron, just with reduced efficiency.
In a general catalytic system, it is, however, natural to consider having a distribution of different values of nc. On the one hand, the dissociation process could be driven by an external (day–night) cycling or by tuning of the intrinsic timescale of time-dependent particle interactions. On the other hand, the size of the monomer network that assembles on the catalyst in some fixed time interval is random, so even a fixed timescale for dissociation leads to random values of nc.
Therefore, in our simulations that demonstrate emergent catalytic reactions (Fig. 3), we replace the various time-dependent factors of dissociation of a monomer network by simply taking the number of bonds in the network to be random within a prescribed limited distribution. The precise procedure is as follows: Given any cluster, at every time step, we count how many monomers are attached to it, say, M, and then dissociation happens if the number M+3nc, where nc is a random number drawn from a prescribed distribution. The distribution is nonzero for nc=514, and peaked at ncmax=12. Our choice of demanding M+3 bonds among M monomers provides a balance between quicker templating (correlates with less bonds) and more successful folding into a rigid cluster (correlates with more bonds). The choice ncmax=12 sets a preference for templating clusters of size N=9, which are the smallest possible clusters necessary to lead to templating a N=11 cluster. No rigid cluster of N>11 can result from templating of the given particle species (SI Text, Maximal Cluster Size Consistent with the Interaction Matrix of N0 = 7 Parent Catalyst).

Maximal Cluster Size Consistent with the Interaction Matrix of N0 = 7 Parent Catalyst

Let us consider all of the possible products of catalysis in a system that starts from an N0=7 parent catalyst. The particle species are then determined by the interaction matrix of the N0=7 parent cluster. Since the valence of particles of species labeled 1–6 are fixed at one, any subsequent product of catalysis can contain at most a single particle of each of the flavors 1–6. Crucially, due to valence of two of the particle of the 7th species, products of continued catalysis can contain any number of particles of species 7.
In any given catalysis product, each of the particles of species 7 has to have at least 3 contacts (for rigidness), which must be with particles of species 1–6, since particles of the same species repel each other. Therefore, we need to determine the triangular facets, i.e., arrangements of three particles to which a fourth one can be rigidly attached, that can be formed out of species 1–6 particles.
Let us first consider building up from all rigid clusters that can be made out of N=6 particles. These are polytetrahedron and octahedron. Starting from the octahedron, we first note that it has 8 triangular facets. According to the interaction matrix, the only way to build it is to use 4 particles having species, e.g., 2, 3, 4, and 5, and additionally use 2 particles of species 7 (Fig. 1D). Since the latter two particles participate in each of the 8 triangular facets of the octahedron, we cannot attach any more of the species 7 particles. However, to this octahedron we can attach particles of the unused species (in this example, 1 and 6), remembering that valence limits us to just one particle of each. Checking all of the facets on which these can be attached, we find that the two attached particles create maximally 3 more spatial positions where additional species 7 particles could be added. In total, this would be 5 particles of species 7, forming a maximally sized rigid cluster at N=11.
In case of the polytetrahedron, there are multiple colorings, i.e., choices of particle species that form it, yet the above reasoning leads to the same result of maximally N=11 particle rigid cluster.
In fact, we applied this argument to all rigid clusters of N>6 particles and confirmed that N=11 is the largest rigid cluster that can form in the catalytic system generated by the N0=7 parent catalyst.

Templating Without Valence 2

Fig. S3 shows all possible templating processes in a system starting from an N0=7 parent cluster, with the restriction that each particle of the catalyst can bind at most one monomer at a time, i.e., valence of all particle species is 1. The limit on valence does not allow templating any structure with N>7. We consider all of the possible products of catalysis in terms of interaction matrix and the geometry of templated monomer network, just as in SI Text, Possible Catalysts for the Octahedron. Note that there are no catalytic cycles in this case.

Catalytic Relations with Complementary Species

In our system, due to the complementary nature of stickers on particles, each particle species has a complementary one: See Fig. 3, pairing the first species with the eighth, the second with the ninth, and so on. The sets of complementary species have the exact same mutual interactions. Physically, this means that every cluster has exactly one additional complementary version. Listing all of the rigid clusters that can result from catalysis in our system, the above implies that the list in Fig. 4 is doubled, with each cluster having a second version built of particles of complementary species.
We now compare the simplified description in Emergent Catalytic Reactions to the full set of clusters and the full graph of their templating connections. For convenience, we here label the 135 clusters by α, and their complementary versions by β. Since our templating is based on attachment of complementary monomers to the catalyst, the templated cluster always contains the complementary species compared with the catalyst. This means that every templating connection in Fig. 4 still exists but goes from the α catalyst to the β product. For every such templating connection, there is another one, from the catalyst version β to the product version α. Therefore, the total number of nodes and the total number of edges in the graph are doubled compared with the situation in the main text.
Most importantly, we consider the catalytic cycles in the graph. It is obvious that a simple catalytic cycle involving an even number of clusters from the main text has an even number of catalytic connections. If clusters switch from α to β in every reaction, the length and nature of such a cycle is preserved. In fact, there is a second cycle obtained by exchanging α and β for each cluster. Therefore, the number of even cycles also doubles with the number of clusters.
On the other hand, every simple catalytic cycle involving an odd number of clusters in the main text gives rise to a hypercycle of double the length. Namely, starting from an arbitrarily chosen initial cluster in the cycle, of version α, after odd number of reactions we obtain the initial cluster in the β version. The catalytic cycle can continue another round, involving all of the complementary versions, until we finish with the initial α cluster. The simplified odd cycle of the main text is in this way replaced by a hypercycle, which contains two copies of it executed in complementary versions of the clusters. The hypercycle therefore does not have a complementary version. In conclusion, the doubling of the number of clusters implies the replacement of odd cycles by hypercycles of double the length.
In terms of the number of catalytic cycles per cluster and in the overall connectivity of the graph of templating connections, we see that simplified description in the main text is similar to the detailed version, which includes complementary species. We therefore expect that all our results are fully applicable to the more detailed description.

Kinetic Rate Constants

We consider two models for the matrix, 𝐌, of first-order production rates, whose element Mij equals the kinetic rate constant of the reaction in which the j-th cluster acts as a catalyst for templating the i-th cluster (i,j=1...135). If there is no templating connection from j to i in the graph (Figs. 4 and 5), that Mij is always strictly zero. In the following, we describe models exclusively for the other (nonzero) entries of 𝐌.
In the first model, we consider all (nonzero) kinetic rate constants to be of the same order. We assign them random values from the uniform distribution between 0.9 and 1.1. The randomness represents the existence of many varying factors in templating, e.g., entropic, dynamic, and so on. Numerically, the limit of all kinetic rate constants being strictly equal leads to large degeneracies and is therefore undesirable. The amount of randomness we used induces small quantitative changes of the eigenvalues and the structure of the dominating eigenvector (Exponential Growth of Catalytic Cycles).
In the second model, each (nonzero) kinetic rate constant is a sum of two independent terms. The first term is the random term described in the previous paragraph. For the second term, we divide the Mij into blocks, for which i and j clusters have some fixed sizes. The second term has the same value for all rate constants within a block. To each block, we assign a random number drawn from an exponential distribution with exponent 2. The second term roughly corresponds to geometrical effects in templating, e.g. the fact that the average time needed for building a monomer network of certain size on the catalyst depends on the number of monomers involved. Note that the second term here is, on average, larger than the first term. We find that the structure of the dominating eigenvector depends on which matrix blocks got assigned the highest values of the second term. On the other hand, we find that the realness of the highest eigenvalue and the fact that it is much larger than the real parts of all other eigenvalues, do not depend on the random realizations of the matrix in this model. For the results presented in Exponential Growth of Catalytic Cycles, the two matrix blocks with the highest values represent the production of N=11 from N=9 catalysts and vice versa.


Z.Z. thanks Philippe Nghe for useful discussions. This research was supported by National Science Foundation Grant DMR-1435964; Harvard Materials Research Science and Engineering Center Grant DMR-1420570; and Division of Mathematical Sciences Grant DMS-1411694. M.P.B. is an investigator of the Simons Foundation.

Supporting Information

Supporting Information (PDF)


E Schrödinger What Is Life? (Cambridge Univ Press, Cambridge, UK, 1944).
FJ Dyson, A model for the origin of life. J Mol Evol 18, 344–350 (1982).
AI Oparin The Origin of Life (Moscow Worker, Moscow, 1924).
D Segre, D Lancet, Composing life. EMBO Rep 1, 217–222 (2000).
SA Kauffman, Autocatalytic sets of proteins. J Theor Biol 119, 1–24 (1986).
JD Farmer, SA Kauffman, NH Packard, Autocatalytic replication of polymers. Physica D 22, 50–67 (1986).
E Mossel, M Steel, Random biochemical networks: The probability of self-sustaining autocatalysis. J Theor Biol 233, 327–336 (2005).
Z Zeravcic, VN Manoharan, MP Brenner, Size limits of self-assembled colloidal structures made using specific interactions. Proc Natl Acad Sci USA 111, 15918–15923 (2014).
CA Mirkin, RL Letsinger, RC Mucic, JJ Storhoff, A DNA-based method for rationally assembling nanoparticles into macroscopic materials. Nature 382, 607–609 (1996).
AP Alivisatos, et al., Organization of ‘nanocrystal molecules’ using DNA. Nature 382, 609–611 (1996).
M Valignat, O Theodoly, J Crocker, Reversible self-assembly and directed assembly of DNA-linked micrometer-sized colloids. Proc Natl Acad Sci USA 102, 4225–4229 (2005).
PL Biancaniello, AJ Kim, JC Crocker, Colloidal interactions and self-assembly using DNA hybridization. Phys Rev Let 94, 58302 (2005).
JT McGinley, I Jenkins, T Sinno, JC Crocker, Assembling colloidal clusters using crystalline templates and reprogrammable DNA interactions. Soft Matter 9, 9119–9128 (2013).
FJ Martinez-Veracoechea, BM Mladek, AV Tkachenko, D Frenkel, Design rule for colloidal crystals of DNA-functionalized particles. Phys Rev Lett 107, 045902 (2011).
PE Theodorakis, C Dellago, G Kahl, A coarse-grained model for DNA-functionalized spherical colloids, revisited: Effective pair potential from parallel replica simulations. J Chem Phys 138, 025101–025113 (2013).
JD Halverson, AV Tkachenko, DNA-programmed mesoscopic architecture. Phys Rev E Stat Nonlin Soft Matter Phys 87, 062310 (2013).
LD Michele, E Eiser, Developments in understanding and controlling self assembly of DNA-functionalized colloids. Phys Chem Chem Phys 15, 3115–3129 (2013).
L Feng, LL Pontani, R Dreyfus, P Chaikin, J Brujic, Specificity, flexibility and valence of DNA bonds guide emulsion architecture. Soft Matter 9, 9816 (2013).
S Sahu, P Yin, JH Reif, A self-assembly model of time-dependent glue strength. Algorithmic Bioprocesses, eds A Condon, D Harel, JN Kok, A Salomaa, E Winfree (Springer, Berlin), pp. 185–204 (2009).
D Yao, et al., Integrating DNA-strand-displacement circuitry with self-assembly of spherical nucleic acids. J Am Chem Soc 137, 14107–14113 (2015).
N Arkus, VN Manoharan, MP Brenner, Minimal energy clusters of hard spheres with short range attractions. Phys Rev Lett 103, 118303 (2009).
N Arkus, VN Manoharan, MP Brenner, Deriving finite sphere packings. SIAM J Discrete Math 25, 1860–1901 (2011).
MC Holmes-Cerfon, Enumerating rigid sphere packings. SIAM Rev 58, 229–244 (2016).
S Hormoz, MP Brenner, Design principles for self-assembly with short-range interactions. Proc Natl Acad Sci USA 108, 5193–5198 (2011).
WM Jacobs, A Reinhardt, D Frenkel, Rational design of self-assembly pathways for complex multicomponent structures. Proc Natl Acad Sci USA 112, 6313–6318 (2015).
A Reinhardt, D Frenkel, Numerical evidence for nucleated self-assembly of DNA brick structures. Phys Rev Lett 112, 238103 (2014).
B Alberts, et al. Molecular Biology of the Cell (Garland Science, 5th Ed, New York, 2008).
G Swiegers Mechanical Catalysis: Methods of Enzymatic, Homogeneous, and Heterogeneous Catalysis (John Wiley & Sons, New York, 2008).
Z Zeravcic, MP Brenner, Self-replicating colloidal clusters. Proc Natl Acad Sci USA 111, 1748–1753 (2014).
P Nghe, et al., Prebiotic network evolution: Six key parameters. Mol Biosyst 11, 3206–3217 (2015).
AV Tkachenko, S Maslov, Spontaneous emergence of autocatalytic information-coding polymers. J Chem Phys 143, 045102 (2015).
S Angioletti-Uberti, P Varilly, BM Mognetti, D Frenkel, Mobile linkers on DNA-coated colloids: Valency without patches. Phys Rev Lett 113, 1–14 (2014).
JD Halverson, AV Tkachenko, Sequential programmable self-assembly: Role of cooperative interactions. J Chem Phys 144, 094903 (2016).
ME Leunissen, R Dreyfus, NC Seeman, DJ Pine, PM Chaikin, Towards self-replicating materials of DNA-functionalized colloids. Soft Matter 5, 2422 (2009).
T Wang, et al., Self-replication of information-bearing nanoscale patterns. Nature 478, 225–228 (2011).
AJ Turberfield, B Yurke, AP Mills, FC Simmel, JL Neumann, A DNA-fuelled molecular machine made of DNA. Nature 406, 605–608 (2000).

Information & Authors


Published in

Go to Proceedings of the National Academy of Sciences
Go to Proceedings of the National Academy of Sciences
Proceedings of the National Academy of Sciences
Vol. 114 | No. 17
April 25, 2017
PubMed: 28396424


Submission history

Published online: April 10, 2017
Published in issue: April 25, 2017


  1. catalytic cycles
  2. DNA-coated colloids
  3. metabolism
  4. templating


Z.Z. thanks Philippe Nghe for useful discussions. This research was supported by National Science Foundation Grant DMR-1435964; Harvard Materials Research Science and Engineering Center Grant DMR-1420570; and Division of Mathematical Sciences Grant DMS-1411694. M.P.B. is an investigator of the Simons Foundation.


This article is a PNAS Direct Submission. A.V.T. is a Guest Editor invited by the Editorial Board.



Zorana Zeravcic1 [email protected]
Soft Matter and Chemistry Department, École Supérieure de Physique et de Chimie Industrielles de la Ville de Paris, Paris Sciences et Lettres Research University, 75005 Paris, France;
Living Matter Laboratory, The Rockefeller University, New York, NY 10065;
Harvard John A. Paulson School of Engineering and Applied Sciences, Harvard University, Cambridge, MA 02138;
Kavli Institute of Bionano Science and Technology, Harvard University, Cambridge, MA 02138
Michael P. Brenner
Harvard John A. Paulson School of Engineering and Applied Sciences, Harvard University, Cambridge, MA 02138;
Kavli Institute of Bionano Science and Technology, Harvard University, Cambridge, MA 02138


To whom correspondence should be addressed. Email: [email protected].
Author contributions: Z.Z. and M.P.B. designed research, performed research, analyzed data, and wrote the paper.

Competing Interests

The authors declare no conflict of interest.

Metrics & Citations


Note: The article usage is presented with a three- to four-day delay and will update daily once available. Due to ths delay, usage data will not appear immediately following publication. Citation information is sourced from Crossref Cited-by service.

Citation statements



If you have the appropriate software installed, you can download article citation data to the citation manager of your choice. Simply select your manager software from the list below and click Download.

Cited by


    View Options

    View options

    PDF format

    Download this article as a PDF file


    Get Access

    Login options

    Check if you have access through your login credentials or your institution to get full access on this article.

    Personal login Institutional Login

    Recommend to a librarian

    Recommend PNAS to a Librarian

    Purchase options

    Purchase this article to get full access to it.

    Single Article Purchase

    Spontaneous emergence of catalytic cycles with colloidal spheres
    Proceedings of the National Academy of Sciences
    • Vol. 114
    • No. 17
    • pp. 4267-E3585







    Share article link

    Share on social media