Self-replicating colloidal clusters Zorana Zeravcica,b,1 and Michael P. Brennera,b a

School of Engineering and Applied Sciences and bKavli Institute for Bionano Science and Technology, Harvard University, Cambridge, MA 02138

We construct schemes for self-replicating clusters of spherical particles, validated with computer simulations in a finite-temperature heat bath. Each particle has stickers uniformly distributed over its surface, and the rules for self-replication are encoded into the specificity and strength of interactions. Geometrical constraints imply that a compact cluster can copy itself only with help of a catalyst, a smaller cluster that increases the surface area to form a template. Replication efficiency requires optimizing interaction energies to destabilize all kinetic traps along the reaction pathway, as well as initiating a trigger event that specifies when the new cluster disassociates from its parent. Although there is a reasonably wide parameter range for self-replication, there is a subtle balance between the speed of the reaction, and the error rate. As a proof of principle, we construct interactions that self-replicate an octahedron, requiring a two-particle dimer for a catalyst. The resulting selfreplication scheme is a hypercycle, and computer simulations confirm the exponential growth of both octahedron and catalyst replicas. self-assembly

| catalytic cycle

T

he ability to invent materials that replicate themselves could lead to a paradigm shift in materials discovery. The exponential amplification of biological molecules, followed by mutation and selection, has allowed the development of powerful protocols for evolving proteins with improved catalytic properties (1, 2). However, whereas modern materials science excels at synthesis, there has been much less success at building artificial self-replicating materials. If there were methods for self-replicating complex materials, then selection–amplification cycles would surely discover solutions with greatly enhanced properties. The logical framework required for objects to replicate themselves was given by von Neumann. He gave an explicit construction of a self-replicating 2D lattice of coupled cellular automata (3), in which a finite area of the lattice replicates itself on an adjacent region. Over the years, von Neumann’s schemes have been refined (4–8), but an efficiently self-replicating artificial system has never been physically realized. To date, artificial selfreplicating systems have focused on linear chain-like structures, where the replicate is templated on the original, closely analogous to DNA, whose replication machinery was unknown at the time of von Neumann’s writings. For example, Seeman and coworkers (9) have recently described and experimentally demonstrated the creation of a self-replicating DNA-based material, using tile motifs with specific binding at their edges and faces. In their scheme, the replication of each generation is achieved by manipulating chemistry and temperature, with complementary motifs manually separated from the initial sequence. The newly formed generation is an accurate copy about one-third of time. Although this is a step in the desired direction, this setup does not achieve exponential growth, and moreover this type of protocol is restricted to linear chain-like molecules. In this paper, we present explicit examples of physical interactions between spheres in a finite-temperature bath that lead to self-replication of clusters, in a manner that is geometrically distinct from replicating reactions of linear chains. Although our examples were inspired by recent experiments with colloidal spheres or nanoparticles coated with DNA, we believe that our construction is sufficiently general that it might prove implementable with a broader class of materials, such as protein complexes

www.pnas.org/cgi/doi/10.1073/pnas.1313601111

with designed interactions (10). We make no attempt to construct a self-replicating material of minimal complexity, but instead endeavor to show that efficient schemes exist, where the errors caused by thermal fluctuations into misfolded states do not inhibit the replication reaction. In our examples, each sphere has stickers distributed uniformly over its surface, causing short-ranged, specific interactions with other spheres. The rules for self-replication are encoded into the specificity and strength of the interactions. Bond strengths must be chosen not only to respect the geometry of the desired self-replication reaction, but also, critically, to avoid kinetic traps. We find that self-replicating a compact cluster requires a specifically designed catalyst, a smaller cluster that allows the parent to template itself from a finite-concentration monomer bath. Additionally, efficient self-replication requires specifying an allosteric trigger event, where the daughter structure separates from the parent, and the bonds within the daughter stabilize. The self-replication reactions we outline replicate both the parent cluster and the catalyst, so that a single catalyst and cluster are sufficient for exponential growth of both. Computer simulations of interacting particles in a thermal bath verify that the design interactions lead to efficient self-replication. For definiteness, we focus on a detailed scheme for the specific case of a self-replicating octahedron, as shown in Fig. 1A. The first step is the design of the initial seed octahedron (red particles). Octahedron is one of the two rigid (ground-state) structures, having 12 contacts, that can be formed out of 6 colloidal particles, the other being a polytetrahedron. When self-assembling from identical particles, the yield of octahedra is ∼4% and that of the polytetrahedra is ∼96%, the difference stemming from rotational entropy (11, 12). If, however, particles are not identical but have stickers that are chosen to satisfy certain interaction rules (13), octahedron becomes the only ground-state structure that can be self-assembled, and therefore the yield is improved Significance One of the hallmarks of living systems is self-replication. Mimicking nature’s ability to self-replicate would not only give more insight into biological mechanisms of self-replication but also could potentially revolutionize material science and nanotechnology. Over the past 60 y, much research, both theoretical and experimental, has been focused on understanding and realizing self-replicating systems. However, artificial systems that efficiently self-replicate remained elusive. In this paper, we construct schemes for self-replication of small clusters of isotropic particles. By manipulating the energy landscape of the process, we show how exponential replication can be achieved. As a proof of principle, we show exponential selfreplication of an octahedral cluster using finite-temperature computer simulations. Author contributions: Z.Z. and M.P.B. designed research, performed research, analyzed data, and wrote the paper. The authors declare no conflict of interest. This article is a PNAS Direct Submission. 1

To whom correspondence should be addressed. E-mail: [email protected].

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10. 1073/pnas.1313601111/-/DCSupplemental.

PNAS Early Edition | 1 of 6

PHYSICS

Edited by Paul M. Chaikin, New York University, New York, NY, and approved December 23, 2013 (received for review July 18, 2013)

Ia

B Oh

tra

p

A

tic

ki

A

B

C

C

A

ne

B* C*

B*

I

B

II

C* C

C 1x

2x

m is

IIa

fo

l

di

Ground state

ng

III

1x

2x

2x

2x

Local minimum

Fig. 1. Self-replication scheme for an octahedron. (A) Self-sustained reaction scheme for self-replication of an octahedron Oh (red particles). A minimal catalyst needed for the replication is a dimer (blue particles). Only one appropriately labeled monomer (gray particles) per particle in Oh/catalyst can be attached (stage Ia). Attached monomers can now interact between themselves and form substructures templated by the Oh and the catalyst (stage II). Once a certain number of bonds between these monomers forms, melting occurs, separating starting Oh and catalyst from attached monomers (stage IIa). New structures formed produce a new catalyst and a nonrigid cluster that folds into a new octahedron (stage III). Along this cycle, there are two undesirable branches: (i) after stage Ia, the Oh–monomer complex can end up in a kinetic trap, and (ii) after the melting occurs (stage II), the new structure can misfold into a local minimum. (B) Interaction matrices for the octahedron and the catalyst. Bonds within these clusters are irreversible. (C) All 9-contact floppy transient structures made by removing one bond from the 10-contact floppy transient structure. Numbers next to these structures indicate their stoichiometric factors. Nine out of 10 structures fold into an octahedron. The one framed can also fold into the lowest-energy local minimum.

(see SI Text and Figs. S1 and S2 for more details). This choice of interaction rules (Fig. 1B) has three types of particles, labeled blue (A), green (B), and red (C), where like particles do not stick to each other, but do stick to the other particle types. Generally, a scheme for self-replication requires the seed octahedron to act as a template for the formation of another octahedron. Each particle in the seed can bind to another (complementary) particle from the solution of monomers. Attached particles can interact between themselves, resulting in templated substructures. However, these interactions alone are still insufficient for producing an octahedron, because the geometrical backbone of the seed is too rigid for the attached particles to flex into an octahedral shape. This can instead be achieved with the presence of a catalyst, and we find that it is sufficient for the catalyst to consist of a dimer (blue particles in Fig. 1A). The dimer has two particle types labeled purple ðC* Þ and orange ðB* Þ, with the interaction matrix shown in Fig. 1B. As shown in Fig. 1A, the catalyst can bind two (appropriate) monomers from the solution (stage Ia) and together with four particles bound to the seed octahedron form a six-particle monomer structure (stage II), which gives rise to the octahedron replica when “melting” is triggered (stage IIa and III). This cycle is self-sustained because the remaining two monomers bound to the seed octahedron become a catalyst dimer after melting (stage IIa and III). The melting process refers to breaking of the single bond that connects a monomer to the seed octahedron or the catalyst, and is triggered for all monomers in stage II when the total number of bonds between them reaches a predetermined value nc . Once melting occurs, we assume that the bonds formed between the monomers become irreversible. Setting the value of nc is crucial because it strongly influences the balance between replication efficiency and error rate. By the very nature of templating from a 3D cluster, the melting produces a floppy transient structure of monomers (having nc bonds), which ideally should fold into a geometrical copy of the octahedron (see stage IIa in Fig. 1A). In our construction 2 of 6 | www.pnas.org/cgi/doi/10.1073/pnas.1313601111

(Fig. 1A), the maximum geometrically possible value of nmax is c 10. For nc = 10, the transient structure can only fold into the octahedron. For nc = 9, Fig. 1C shows that 9 out of 10 possible outcomes of folding is the octahedron. Although this gives a rough estimate of the error rate to be 10%, our simulations indicate it is smaller, implying that entropic factors are significant. In general, with decreasing nc , the melting happens more frequently but with a higher error rate. An optimal value of nc depends on the energy landscape of the cluster in question, and we will discuss this in more detail below. The geometric scheme for self-replication analyzed so far (Fig. 1A) can be turned into a precise reaction scheme (Fig. 2A) with 10 different particle types. The particles A, B, and C in the seed octahedron can bind to their complementary particle type A′, B′, and C′ from the solution of monomers if the particles were coated with DNA, a monomer attaching to a seed particle would have the appropriate complementary strand. Similarly, the catalyst particles C* and B* bind to complementary monomers C*′ and B*′ from the solution. The complete interaction matrix for all particle types is shown in Fig. 2B. The interaction rules imply that an A; B; C-type octahedron seed with a C* ; B* -type dimer creates a new octahedron made of particle types A′; B′; C′ and a new dimer of type C*′ ; B*′ . The newly created octahedron and dimer can now bind monomers and take part in a melting process that will create their complements, i.e., an A; B; C-type octahedron and a C* ; B* -type dimer, closing the twostep replication cycle. As Fig. 2A illustrates, the entire self-replication process of two complementary types of octahedra is a hypercycle (14–17), consisting of two linked catalytic cycles. We test the reaction scheme using dissipative particle dynamics (DPD) techniques (18, 19). Our simulation contains colloidal spheres of diameter D, with an interaction range of 1:05D. The colloids are immersed into a DPD solvent of smaller particles. We simulated systems with 96, 256, and 512 colloidal particles, out of which 6 comprise the parent octahedron and 2 the catalyst. Simulations are run at a fixed temperature with Zeravcic and Brenner

A

B IIIa

Ia

II

I

A

A C

B*

C A

C

B

A’ B’

B’

C’

A

A’

C’

B

C

C*

IV

III A’

C’

B’

B

B

A B C B*’ C*’ A’ B’ C’ B* C*

C*’

A’

B*’

B’

B B*’

C’ A’

C*’

C

C’

B’

C’ A’

C*

B*

A B C

IIa

IVa

A B’

B C

A

C* B*

B*’

C*’

B*’ C*’ A’ B’ C’ B* C*

a volume fraction of colloids ϕcoll = 0:01, and a larger volume fraction of solvent ϕsol ≈ 0:2. More details are given in SI Text. Even with a geometrically compatible scheme for self-replication, a successful implementation still requires an extensive search for both bond strengths and the characteristics of the disassociation reaction nc . With bond optimization, and the value

Oh Clusters

A

20.0 15.0

nc = 9, Fig. 3 presents typical snapshots of a simulation showing successful self-replication, where different phases of the reaction cycle can be recognized. Fig. 3A shows exponential growth of the number of octahedra, as expected for a hypercycle. Such robust results rely critically on avoiding traps along the self-replication reaction pathway, requiring careful choice of

B1

B2

B4

B5

B7

B8

10.0 7.0 5.0 3.0 2.0 1.5 1.0 0

B3

B6

2

4

6

8 10 12 14

time s

Fig. 3. DPD simulation of octahedron self-replication using 512 colloidal particles. (A) The total number of octahedra (both ABC and A′B′C′) as function of time. The black line ∼ ex is a guide for the eye. (B1–B8) Typical snapshots of a simulation (see SI Text for the movie). Colloids are colored if they comprise an Oh or a catalyst (color scheme like in Fig. 2A). Monomers are in general not shown unless they are attached to an Oh/catalyst (transparent spheres). These snapshots reveal phases of the reaction cycle (Fig. 2A). For example, B3 shows a melting process, after which in B4 a new octahedron is formed together with a new catalyst. In B8, there are 12 octahedral clusters apart from the initial seed. One of the replicas (marked with a red arrow) is a misfolded octahedron.

Zeravcic and Brenner

PNAS Early Edition | 3 of 6

PHYSICS

Fig. 2. A detailed self-replication scheme for an octahedron. (A) The reaction scheme. (B) Full interaction matrix between particle types present in the simulation. Interactions between particle types A=A′ and B,C,B* ,C* =B′,C′,B*′ ,C*′ (colored light orange) are weaker than the other interactions (dark orange). The ratio of these bond strengths is optimized (SI Text). Once particles form an octahedron or a catalyst, their bonds become irreversible.

4 of 6 | www.pnas.org/cgi/doi/10.1073/pnas.1313601111

kkinetic trap = const 1

0.02 0.08

0

Log (kmisfolding)

both bond strengths and the melting criterion. The first kinetic trap that can impede the self-replication of the octahedron occurs during the initial attachment of the monomers (after stage Ia in Fig. 1A); the attached monomers can form two triangles attached on opposite faces of the octahedron. This configuration is relatively frequent and stable due to the octahedral symmetry. The bond optimization must break this symmetry of bonds within a triangle, disfavoring this configuration. Kinetic traps in a later stage of the hypercycle (after stage II in Fig. 1A) can also occur due to the choice of nc . The maximum value of nc is determined by the geometry of the template. For the octahedron with a dimer enzyme, nmax = 10. There is also a minimum value of nc = ngs c c , where the floppy transient structure obtained by melting can fold only into the desired cluster. The value of ngs c depends on the geometry of the lowest energy local minima where kinetic traps can occur. For the octahedron, the lowest energy local minimum (Fig. 1C) has two bonds less (10) than the ground state (12), and three bonds need to be broken to transform the ground state into that local minimum. This implies ngs c = ð12 − 3Þ + 1 = 10. The mentioned lowest-energy local minimum can actually still act as a template for replicating an octahedron, albeit with a smaller replication rate because, e.g., the red–green–red–green face and blue–green–blue–green face (Fig. 1C) cannot template a replica. The property that the local minimum can only template a copy of an octahedron is exceptional, and holds because the octahedron is a small and highly symmetric cluster. How large is the parameter range available for self-replication? One way to assess this is to consider a situation where there is a fixed concentration of the different monomers, and the catalyst dimer, and write the kinetic differential equations that govern transitions between different states of the hypercycle as labeled in Fig. 1A (see SI Text for explicit equations). Several of the parameters (e.g., the attachment of monomers onto the parent octahedron) are set by diffusion, so as a choice of normalization we set these equal to unity. Three key parameters can be tuned with the bond strengths and nc : (i) the rate kkinetic trap , at which the stage I octahedron transitions to the initial kinetic trap; (ii) the rate kmelting , at which the melting criterion is satisfied from the stage II octahedron–catalyst complex; and (iii) the rate kmisfolding , at which the melted structure misfolds into something other than an octahedron. Keeping in mind that there is a lower bound on the bond strengths set by the temperature of the experiment, we can choose bond strengths so that the initial kinetic trap is effectively suppressed. By solving the differential equations (SI Text) and extracting the exponential growth rate of octahedra, we arrive at the replication phase diagram in Fig. 4, which shows that exponential self-replication occurs whenever misfolding is sufficiently low, independently of the value of melting rate constant. Generally, with a more stringent melting criterion (i.e., larger nc ), kmisfolding decreases because there are less pathways available for the floppy structure to misfold. Quantitatively, the misfolding rate is roughly the ratio of numbers of ways the floppy structure can fold into a local minimum versus into the desired cluster. Fig. 1C shows that, for the octahedron, if nc = 10, kmisfolding = 0, if nc = 9, kmisfolding ∼ 1=10, if nc = 8, kmisfolding ∼ 1=6, and so on, where the rates are measured in the units of a typical rate constant set by diffusion. This means that small kmisfolding (i.e., exponential replication) can be generically achieved by simply choosing a sufficiently large nc . The rate kmelting likewise decreases with increasing nc , and hence is not independent of kmisfolding . Fig. 4 also shows that, once kmisfolding allows exponential replication, kmelting just sets the exponential growth rate of the replication reaction. These conclusions are robust for a wide range of kkinetic trap , confirming that the melting and misfolding rates basically control the efficiency of self-replication. As a general principle, the second role of bond optimization is to increase the probability of formation of transient nc -contact structures that can fold into the desired cluster. In the octahedron

0.04

0

1

2 0.02

3

3

2

1

0

1

Log(kmelting) Fig. 4. Contour map of the replication growth rate as a function of melting and misfolding rates. Kinetic trap rate is kept constant. The blue contours correspond to negative growth and pink to positive growth rates.

example, using the value nc = 9, this effect is less prominent. Our bond optimization is naively not expected to suppress the undesired nine-contact transient structure (marked red in Fig. 1C), but still the replication rate is much improved due to suppression of the kinetic trap. For a given cluster and choice of nc , bond optimization can be guided by suppression of kinetic traps and boosting of final folding, although possibly only a balance between the two will give the optimal replication and error rates. It is worth noting that the existence of the trade-off between misfolding (error) rate and melting (replication) rate is well established in protein synthesis and other biosynthetic activities (like chromosome replication and transcription) (20, 21). In our system as well, this balance is intrinsic, reflecting the fact that replication proceeds autonomously. Following the principles explained above, we created schemes for self-replicating three other clusters: (i) two different schemes for the polytetrahedron ðN = 6Þ (Fig. 5 A and B) requiring a dimer catalyst; (ii) pentagonal bipyramid ðN = 7Þ (Fig. 5C) requiring a trimer catalyst; (iii) chiral (N = 7) cluster (Fig. 5D) requiring a dimer catalyst. It is noteworthy that the N = 7 chiral = 14. cluster can replicate, preserving its chirality when nc = nmax c If nc ≤ 13, the replication process creates copies with both chiralities. These examples suggest that the ability of a cluster to self-replicate might not depend so much on its geometry. The number of different particle types required for realizing any of these schemes depends on the design of actual clusters. Like the octahedron, to ensure that a desired cluster is the only ground state, we need to use a number of different particle types obeying certain interaction rules. For a given cluster, the solution to this problem is not necessarily unique (13). Even without going into details, we can give upper bounds for the number of different particle types, M, needed for self-replicating schemes in Fig. 5: (i) the polytetrahedron (schemes a and b) can be designed with at most 5 different particle types, implying the total number of different particle types is M = 2 · 5 + 2 · 2 = 14, where we account for the cluster, its catalyst, and their complements; (ii) the pentagonal bipyramid (scheme C) can be designed with at most 6 particle types, implying M = 2 · 6 + 2 · 3 = 18 because the catalyst is a trimer; (iii) the chiral cluster (scheme D) can be designed with at most 7 particle types, implying M = 2 · 7 + 2 · 2 = 18. Next, we discuss possible experimental realizations of our protocol. The replication schemas described here were inspired by colloidal particles coated with DNA that have already been realized and used for self-assembly (22–26). Although our description assumes particles with isotropic interactions, the use of Zeravcic and Brenner

Ia

A

Ia

B II

I

I

II

IIa

IIa

III

III

C

D Ia Ia

I

II I

II

IIa

IIa

III

Fig. 5. Self-replication schemes for various clusters. (A and B) Two different self-sustained reaction schemes for self-replication of a N = 6 polytetrahedral cluster (red particles). A minimal catalyst needed for the replication is a dimer (blue particles). (C) Self-sustained reaction scheme for self-replication of a N = 7 pentagonal bipyramid (red particles). A minimal catalyst needed for the replication is a trimer (blue particles). (D) Self-sustained reaction scheme for selfreplication of a chiral N = 7 cluster (red particles). A minimal catalyst needed for the replication is a dimer (blue particles).

particles with bond directionality is also possible (27–29) and would lower the number of different particles types required to implement our schemes. The essential ingredients could just as well be encoded with other materials, e.g., proteins (10). The most difficult part of the experimental realization appears to be the melting process, where we require that melting is locally triggered when attached monomers form nc bonds. In general, as nc is being reached, each formed bond decreases the energy of the floppy attached cluster. Once attached monomers have enough contacts between themselves, this floppy cluster signals both the seed cluster and the catalyst to release all of their extra bonds, which requires an input of energy (SI Text). With present technology for building blocks, it is not straightforward to trigger local melting. This would require the development of methods for achieving colloidal interactions with positive and negative cooperative binding and other allosteric mechanisms. Another, currently accessible, way for exploring our schemes is to modulate the particle interactions with light and temperature cycling; this is currently the standard experimental method for achieving replication with linear structures (9). To conclude, we emphasize that the vast majority of replication schemes studied so far, such as those based on linear molecules (9), are inspired by DNA and represent “symbolic 1. Giver L, Gershenson A, Freskgard P-O, Arnold FH (1998) Directed evolution of a thermostable esterase. Proc Natl Acad Sci USA 95(22):12809–12813. 2. Arnold FH (2001) Combinatorial and computational challenges for biocatalyst design. Nature 409(6817):253–257. 3. Von Neumann J, Burks AW (1966) Theory of Self-Reproducing Automata (Univ of Illinois Press, Urbana, IL). 4. Nobili R, Pesavento U (1994) John von Neumann’s Automata Revisited. Artificial Worlds and Urban Studies. eds Bessusi E, Cecchini A (DAEST Publication, Convegni 1, Venezia). 5. Sipper M (1998) Fifty years of research on self-replication: An overview. Artif Life 4(3): 237–257. 6. Chou H-H, Reggia JA (1998) Problem solving during artificial selection of self-replicating loops. Physica D 115:293–312. 7. Penrose LS (1959) Self-reproducing machines. Sci Am 200:105–114. 8. Penrose LS (1958) Mechanics of self-reproduction. Ann Hum Genet 23(1):59–72.

Zeravcic and Brenner

replication,” which is irrespective of size and shape. However, the replication mechanism of other biological objects, such as Golgi apparatus, endoplasmatic reticulum, and micelles (30–33), is firmly tied to their geometrical shape. Our schemas are closer to the latter mechanism of geometry-based replication, which could be a significant advantage in building bulk materials. Finally, we note that, although we have demonstrated the possibility of coding a replication reaction into the interactions between particles in compact clusters, an important open question is whether such schemes can be designed for arbitrarily large compact structures. Materials and Methods A detailed description of our simulation together with two simulation movies (Movies S1 and S2) and a movie that demonstrates the replication scheme shown in Fig. 1A (Movie S3) is included in SI Text. ACKNOWLEDGMENTS. We thank Arvind Murugan, Stanislas Leibler, and Vinothan N. Manoharan for helpful discussions. This research was funded by the George F. Carrier Fellowship, the National Science Foundation through the Harvard Materials Research Science and Engineering Center (DMR0820484), the Division of Mathematical Sciences (DMS-0907985), and by Grant RFP-12-04 from the Foundational Questions in Evolutionary Biology Fund. M.P.B. is an investigator of the Simons Foundation.

9. Wang T, et al. (2011) Self-replication of information-bearing nanoscale patterns. Nature 478(7368):225–228. 10. King NP, et al. (2012) Computational design of self-assembling protein nanomaterials with atomic level accuracy. Science 336(6085):1171–1174. 11. Meng G, Arkus N, Brenner MP, Manoharan VN (2010) The free-energy landscape of clusters of attractive hard spheres. Science 327(5965):560–563. 12. Arkus N, Manoharan VN, Brenner MP (2009) Minimal energy clusters of hard spheres with short range attractions. Phys Rev Lett 103(11):118303. 13. Hormoz S, Brenner MP (2011) Design principles for self-assembly with short-range interactions. Proc Natl Acad Sci USA 108(13):5193–5198. 14. Eigen M, Schuster P (1977) The Hypercyde: A principle of natural self-organization. Naturwissenschaften 64:541–565. 15. Eigen M, Schuster P (1978) The hypercycle part B. Naturwissenschaften 65:7–41. 16. Eigen M, Schuster P (1978) The hypercycle part C. Naturwissenschaften 65:341–369. 17. Müller-Herold U (1983) What is a hypercycle? J Theor Biol 102:569–584.

PNAS Early Edition | 5 of 6

PHYSICS

III

18. Hoogerbrugge PJ, Koelman JMVA (1992) Simulating microscopic hydrodynamic phenomena with dissipative particle dynamics. Europhys Lett 19(3):155–160. 19. Groot R, Warren P (1997) Dissipative particle dynamics: Bridging the gap between atomistic and mesoscopic simulation. J Chem Phys 107(11):4423–4435. 20. Johansson M, Lovmar M, Ehrenberg M (2008) Rate and accuracy of bacterial protein synthesis revisited. Curr Opin Microbiol 11(2):141–147. 21. Murugan A, Huse DA, Leibler S (2012) Speed, dissipation, and error in kinetic proofreading. Proc Natl Acad Sci USA 109(30):12034–12039. 22. Mirkin CA, Letsinger RL, Mucic RC, Storhoff JJ (1996) A DNA-based method for rationally assembling nanoparticles into macroscopic materials. Nature 382(6592):607–609. 23. Valignat MP, Theodoly O, Crocker JC, Russel WB, Chaikin PM (2005) Reversible selfassembly and directed assembly of DNA-linked micrometer-sized colloids. Proc Natl Acad Sci USA 102(12):4225–4229. 24. Biancaniello PL, Kim AJ, Crocker JC (2005) Colloidal interactions and self-assembly using DNA hybridization. Phys Rev Lett 94(5):058302. 25. Park SY, et al. (2008) DNA-programmable nanoparticle crystallization. Nature 451(7178): 553–556.

6 of 6 | www.pnas.org/cgi/doi/10.1073/pnas.1313601111

26. Macfarlane RJ, et al. (2011) Nanoparticle superlattice engineering with DNA. Science 334(6053):204–208. 27. Suzuki K, Hosokawa K, Maeda M (2009) Controlling the number and positions of oligonucleotides on gold nanoparticle surfaces. J Am Chem Soc 131(22): 7518–7519. 28. Kim J-W, Kim J-H, Deaton R (2011) DNA-linked nanoparticle building blocks for programmable matter. Angew Chem Int Ed Engl 50(39):9185–9190. 29. Wang Y, et al. (2012) Colloids with valence and specific directional bonding. Nature 491(7422):51–55. 30. Alberts B, et al. (2002) Molecular Biology of the Cell (Garland Science, New York). 31. Munro S (2002) More than one way to replicate the Golgi apparatus. Nat Cell Biol 4(10):E223–E224. 32. Bachmann PA, Luisi PL, Lang J (1992) Autocatalytic self-replicating micelles as models for prebiotic structures. Nature 357:57–59. 33. Fellermann H, Solé RV (2007) Minimal model of self-replicating nanocells: A physically embodied information-free scenario. Philos Trans R Soc Lond B Biol Sci 362(1486): 1803–1811.

Zeravcic and Brenner

Supporting Information Zeravcic and Brenner 10.1073/pnas.1313601111

 ωij =

SI Text Simulation Details Dissipative Particle Dynamics. We use dissipative particle dynamics (DPD) techniques (1, 2) to test the reaction scheme in Fig. 1A of the main text. Similar to refs. 3 and 4, our system consists of two types of particles: (i) solvent particles, modeled as standard DPD beads with soft repulsion, dissipative and random interaction; (ii) colloidal particles, which are larger DPD beads, that have the conservative force between two colloids replaced with 48 − 96 Lennard–Jones interaction. The dynamics of each particle in our simulation is governed by Newton’s equations of motion as follows:

dri = vi ; dt

mi

dvi = f i; dt

[S1]

where ri is the position vector of particle i, vi is its velocity, and mi is its mass. All of the particles in our simulation have equal mass m = 1 (unit of mass in our simulation). The force acting on particle i is composed out of three parts as follows: fi =

X j≠i

R C FD ij + Fij + Fij ;

[S2]

where FD is the dissipative force, FR is the random force, and FC is the conservative force. The dissipative force is as follows:    D ^rij · vij ^rij ; FD ij = − γω rij

[S3]

where γ is the viscosity coefficient, vij = vi − vj is the relative velocity of particles i and j, rij = jri − rj j is the distance between the centers of particles i and j, ^rij = rij =rij is a unit vector, and ωD is a distance-dependent weight function. The random force is as follows:  pffiffiffiffiffi   FRij = 1= Δt σωR rij θij^rij ; [S4] where σ is the noise strength, ωR is a distance-dependent weight function, Δt is the simulation time step, and θij is a random variable. Instead of taking θij to be a variable with a Gaussian distribution and unit variance, pffiffiffi as is standard in DPD simulations, in this work we use θij = 3ð2 · ζ − 1Þ, where ζ is a uniformly distributed random number ζ ∈ Uð0; 1Þ (5). This choice makes the simulation very efficient and the results are basically indistinguishable from those calculated using Gaussian numbers. To ensure momentum conservation, in DPD algorithms θij = θji . One of the two weight functions ωD and ωR can be chosen arbitrarily, therefore fixing the other weight function (1, 2, 6). To ensure that the system has Gibbsian equilibrium, the viscosity and noise have to be related by a fluctuation dissipation theorem; this leads to the following relations:     2 ωD rij = ωR rij ; σ 2 = 2γkB T=m;

[S5]

where kB is the Boltzmann constant and T is the temperature. In our simulations, we use γ = 10. As is practiced in DPD simulations ωRij = ωij , where ωij is a simple (soft) weighting function that vanishes at some interaction range rcut as follows: Zeravcic and Brenner www.pnas.org/cgi/content/short/1313601111

1 − rij =rcut 0



if rij ≤ rcut otherwise:

[S6]

Because we simulate two types of particles, the interaction range will depend on which pair is interacting. The range of interaction CC = 1:5D for two colloidal particles, where D is rcut equals: (i) rcut SS the diameter of colloids; (ii) rcut = 0:5D for two solvent particles; CS and (iii) rcut = 1:0D for the interaction between a colloidal and a solvent particle. For the same reason, the conservative forces will differ depending on which particles are interacting. The colloid–colloid conservative force is modeled by interaction of 48 − 96 Lennard–Jones spheres with a unit diameter D = 1 and a short interaction range of r CC = 1:05D as follows: 8 ! > > < 96e 1 − 1 ^rij if rij ≤ r CC rij96 rij48 rij FCcc [S7] ij = > > : 0 otherwise; where e is the bond strength (energy scale in our simulations). The conservative force between two solvent (DPD) particles and between a solvent particle and a colloid is a soft repulsion: FCij = aij ωCij^rij ;

[S8]

CS = 0:75D, and aij is the repulsion pawhere ωCij = ωij but with rcut rameter between particles i and j. Following ref. 2, we use aij = 25kB T for both solvent–solvent and solvent–colloid interactions and ρsol = Nsol =V = 3 for the solvent number density, where the volume of the simulation box V follows from the colloid volume fraction ϕcoll as V = ðNcoll πD3 Þ=ð6  ϕcoll Þ. Our simulation box has periodic boundaries in all three dimensions, and we use the linked-cell algorithm to speed up the calculations (7). To integrate equations of motion, we use the standard velocityVerlet algorithm (8). One of the advantages of DPD simulations is the use of large time steps in the integration. However, compared with soft DPD forces, the Lennard–Jones forces require significantly shorter time step for accurate integration. Therefore, to exploit the advantageous efficiency of DPD in our simulation, we use the multiple time step algorithm (7) with Δt = 0:035 for DPD forces and ΔtLJ = 0:0005 for Lennard–Jones forces, where pffiffiffiffiffiffiffiffiffi time is measured in the units of ðD=2Þ m=e.

Calibration of Simulation. The values for interaction ranges we use

are a result of calibration of our simulation based on the experimental results of ref. 9. In that work, the authors experimentally study spontaneous self-assembly of identical colloidal particles into small clusters and measure their yields. With the use of calibrated interaction ranges (values quoted above), our simulation reproduces the experimentally observed yields for all clusters of N = 6, 7, and 8 identical particles. To explain the calibration procedure, we start by considering our simulation results for N = 6 identical particles. As mentioned in the main text, 6 colloidal particles can form 2 rigid structures, having 12 contacts: polytetrahedron and octahedron (Fig. S1B, Insets). The relative yields of the two ground states in equilibrium are experimentally found to be ∼ 96% for the polytetrahedron and consequently ∼ 4% for the octahedron (dashed horizontal lines in Fig. S1B). These values match the statistical mechanics calculations of the ratios of their partition functions (9) (Fig. S1B, solid horizontal lines). The striking difference 1 of 6

between the two yields originates in the different rotational entropy of these clusters (9). In Fig. S1, we show (A) absolute and (B) relative yields as a function of temperature T (in units of bond strength e) for N = 6 particles, obtained from our calibrated simulations. Each data point is an ensemble average of 1,000 different initial condition realizations, run at a fixed temperature T for trun = 12; 000 time steps Δt. For these simulations, we use ϕcoll = 1=30 and Ncoll = 6. At the end of a simulation run, we identify whether a formed structure is one of the rigid clusters: Using the positions of colloids, we form an adjacency matrix A^ (a Ncoll × Ncoll matrix, with an element Aij = 1 if particles i and j are in contact and Aij = 0 otherwise), calculate its eigenvalues, and compare them with the referent values for polytetrahedron and octahedron. The eigenvalues of an adjacency matrix uniquely determine the cluster that is formed by the particles. Our simulations in general exhibit several regimes as a function of temperature T. At the highest temperatures (at T=e & 0:16 in Fig. S1A), the system is in equilibrium, but the bonds between the colloids are short-lived, leading to small absolute yields of the ground states. With a given ensemble size, this also leads to statistical noise in relative yields (Fig. S1B). The equilibrium regime extends down to T=e ∼ 0:1, where the noise in the relative yields is small. Below T=e ∼ 0:1, the relaxation time of clusters becomes comparable to trun , and the results are strongly influenced by kinetic effects. It is important to note that the colloid–colloid conservative interaction range r CC must not be too large, because that would allow bonds between particles comprising a cluster that are geometrically impossible when particles are hard spheres. This interaction range is experimentally measured to be r CC = 1:05D, and we use this value in all our simulations. Finally, we note that, in our simulations, the system completely melts into a gas for T=e > 0:25. Although one might naively expect for this to occur for T=e & 1, the stability of bonds is also influenced by the vibrational frequency in the Lennard–Jones potential well. The Self-Replication Simulation Designing the Octahedron. To ensure that the octahedron is the only ground state that can be assembled out of six particles, we need to introduce specific interactions between particles. Following ref. 10, we choose three types of particles (two particles per type) such that: (i) within the same type, particles interact unfavorably, and (ii) with other types, particles interact favorably (see the right interaction matrix in Fig. S2). To demonstrate that this particular selection of interactions improves the yield of the octahedron, we performed simulations with the same parameters as in the previous section and with imposed interaction rules, where the favorable interactions are the Lennard–Jones described above, and the unfavorable ones are modeled with the repulsive part of the Lennard–Jones potential. In Fig. S2, we show the yield for the designed octahedron: the improvement compared with the identical particle case is clearly visible. Simulation Parameters. We simulated systems with Ncoll = 96, 256, and 512 colloidal particles. In all of our simulations, the volume fraction of colloids is ϕcoll = 0:01 and that of a solvent is ϕsol ≈ 0:2. Six out of starting Ncoll colloids comprise the parent octahedron made out of A − B − C − A − B − C particles, and two, the parent catalyst made out of C* − B* particles. The rest of the Ncoll − 8 colloids are the monomer bath. In principle, the bath should consist of monomers of all different types (A, B, C, A′, B′, etc.), interacting according to the prescribed interaction matrix (Fig. 2B in the main text). We simplify and speed up the simulation by setting all monomers to be of a single new type “0,” which does not interact with itself. This does not change the physics of the system at the simulation temperature because Zeravcic and Brenner www.pnas.org/cgi/content/short/1313601111

monomer bonds are easily broken, unless they are bonds connecting to the particles comprising the parent clusters (see section below for more details). Furthermore, the first time a type “0” monomer comes into contact with a particle in the parent cluster (type A for example), it immediately assumes the complementary type (A′), and retains this type until the end of the simulation. All of our self-replication simulations are run at a fixed temperature T=e = 0:14 (see below for more details) and for trun = 107 time steps or until the monomer bath is depleted. Bond Strength Optimization and Simulation Temperature. The basic and general property of our approach to self-replication is that the self-replication happens at some fixed temperature T, with an appropriate optimization of bond strengths. To achieve efficient self-replication using our schemes, physically it is enough: (i) for the bonds between monomers to be at equilibrium, (ii) for the bonds between parents and attached monomers to be stronger but still at equilibrium at T, (iii) for the newly formed structures to have irreversible bonds, and (iv) to destabilize kinetic traps by using weaker bonds that are near their melting point at T. Along the self-replication reaction pathway, the first kinetic trap occurs during the initial attachment of monomers to the seed octahedron (after stage Ia in Fig. 1A of the main text). Because of the high symmetry of the octahedron cluster, this kinetic trap is frequent and stable. To disfavor it, bonds between attached monomers have to be optimized such that the symmetry of the kinetic trap is broken. In the kinetic trap, the attached monomers form two triangles on opposite faces of the octahedron. This trap can be suppressed by making one of these three particle types bind more weakly to the other two types. We choose the bonds A′ − B′ and A′ − C′ (A − B and A − C) to be weaker compared with the C′ − B′ ðC − BÞ bond. This particular breaking of the symmetry between the colors should be consistent with the particle types in the catalyst: By choosing the catalyst to be of C* − B* (C*′ − B*′ ) type, two out of four monomers that comprise the substructure that will become the octahedron replica are of type A′ (A) (see stage II in Fig. 2A of the main text). Our bond optimization makes it easier for the two A′ (A) to appear on the same side of the parent octahedron and therefore participate in the interaction with the catalyst. All of our self-replication simulations are run at a fixed temperature T=e = 0:14, where e is the bond strength between attached monomers A′ − B′ and A′ − C′ (A − B and A − C). As explained in the previous section, at this temperature, bonds of strength e are in equilibrium but are short-lived. Therefore, the bond strengths between C′ − B′ ðC − BÞ are set to be 1:5e, making them more difficult to break at the simulation temperature. To ensure that the melting does not happen before nc bonds are formed, we make the bonds between the attached monomers and the parent octahedron/catalyst stronger, and set to 2e. These bonds are still in equilibrium at the simulation temperature. Finally, the bonds within a formed octahedron/catalyst, which are irreversible at the simulation temperature, are set to 5e. We have obtained similar outcomes for self-replication by varying bonds strengths in accordance to our design rules described at the beginning of this section. Energetics of the Melting Process. A practical implementation of the melting rule requires the input of energy, because the free energy of the system increases due to the melting event. The input free energy required is Finput = Einput − TΔSinput . Here, the energy for breaking bonds is Einput = ð6 + 2Þ · eM−P , with eM−P , the bond strength between attached monomers and parents, and ΔSinput < 0, the entropy difference from melting. In our model, Einput is of order 16e, where e is the weaker bond between attached monomers; this is a lower bound on the change in free energy. The total free energy difference could be provided by either consumption of fuel (GTP, ATP) during the melting event or by, 2 of 6

e.g., cyclically changing the temperature of the overall system. The energy required to facilitate the transition depends on the precise design of the melting process. There are many possible implementations; for example, the melting could be cooperative. Movies of Our Simulations. We include here two movies of the selfreplication simulations described above: (i) as a clear visualization of the replication process, we include Movie S1 with Ncoll = 16, which shows one-half of the replication cycle shown in Fig. 2A of the main text, and (ii) Movie S2 with Ncoll = 512, which shows multiple replication cycles. The data and snapshots in Fig. 3 of the main text are taken from Movie S2. In both movies, the monomers that comprise the bath are not visible unless they attach to the already formed octahedra and catalysts. Apart from the simulation movies, we include a demonstration movie of the replication process (Fig. 1A of the main text) using magnetic sticks and balls (Movie S3).

Kinetic Model In this section, we study conditions for efficient self-replication within a generic and simple kinetic model. For simplicity, we start with the assumption that the catalyst and monomer concentrations are kept fixed. Focusing on the example of octahedron self-replication reaction (Fig. 1A of the main text), we introduce the following reactants: A, octahedron (stage I in Fig. 1A of the main text); D, octahedron with attached monomers (stage Ia); F, structure D with attached catalyst (stage II); G, floppy structure obtained by melting (stage IIa); A~, misfolded floppy structure ~ octahedron with attached (between stage II and IIa); and D, monomers in a kinetic trap (between stage Ia and II). The set of reactions describing the self-replication cycle then follows:

6 Gk!   A~:

[S14]

Note that the reactions for catalysts and monomers are omitted, because their concentrations are assumed constant. The last two reactions describe two generic sources of “error” in the replication process: stuck in a kinetic trap, and folding into a wrong structure, respectively. The above reactions lead to the following rate equations: c_ A = − k1 cA + k3 cF + k4 cG

[S15]

c_D = k1 cA − k2 cD − k5 cD

[S16]

c_ F = k2 cD − k3 cF

[S17]

c_G = k3 cF − k4 cG − k6 cG ;

[S18]

3 Fk!   A+G

[S11]

4   A Gk!

[S12]

5 ~   D Dk!

[S13]

where cI is the concentration of reactant I. The rate constants k1, k2, and k4 are all related to diffusion of particles, and it holds that k1 = k2 ≈ k4 . For simplicity and as a choice of normalization, we set k1 = k2 = k4 ≡1. The rate constant k3 is determined by the melting criterion nc ðkmelting ≡ k3 Þ, whereas the constant k5 corresponds to the rate of getting into the initial kinetic trap ðkkinetic trap ≡ k5 Þ. The rate constant k6 corresponds to misfolding of the melted structure into a local minimum ðkmisfolding ≡ k6 Þ, and this constant is related to the melting criterion nc as well. By solving the differential equations [S15] to [S18] and extracting the exponential growth rate of octahedra, we arrive at the replication phase diagram show in Fig. 3 of the main text. This kinetic model was simplified by assuming constant concentrations of catalyst and monomers. The monomer bath is expected to be large, and large changes in monomer concentration should not be physically relevant. However, the behavior of catalyst concentration is important. Because the above model already leads to exponential growth of cluster, one expects even more efficient replication when the catalyst is replicated as well (i.e., the full replication scheme in Fig. 1A of the main text). When we extend the above analysis by including the catalyst and monomer concentrations, we find that typically both the octahedra and catalyst quickly switch from exponential growth to faster than exponential.

1. Hoogerbrugge PJ, Koelman JMVA (1992) Simulating microscopic hydrodynamic phenomena with dissipative particle dynamics. Europhys Lett 19(3):155–160. 2. Groot R, Warren P (1997) Dissipative particle dynamics: Bridging the gap between atomistic and mesoscopic simulation. J Chem Phys 107(11):4423–4435. 3. Dzwinel W, Yuen DA (2000) A two-level, discrete-particle approach for simulating ordered colloidal structures. J Colloid Interface Sci 225(1):179–190. 4. Whittle M, Travis KP (2010) Dynamic simulations of colloids by core-modified dissipative particle dynamics. J Chem Phys 132(12):124906. 5. Nikunen P, Karttunen M, Vattulainen I (2003) How would you integrate the equations of motion in dissipative particle dynamics simulations? Comput Phys Commun 153:407–423.

6. Español P, Warren P (1995) Statistical mechanics of dissipative particle dynamics. Europhys Lett 30(4):191–196. 7. Allen MP, Tildesley DJ (1987) Computer Simulation of Liquids (Oxford Univ Press, Oxford). 8. Swope WC, Andersen HC, Berens PH, Wilson KR (1982) A computer simulation method for the calculation of equilibrium constants for the formation of physical clusters of molecules: Application to small water clusters. J Chem Phys 76(1):637–649. 9. Meng G, Arkus N, Brenner MP, Manoharan VN (2010) The free-energy landscape of clusters of attractive hard spheres. Science 327(5965):560–563. 10. Hormoz S, Brenner MP (2011) Design principles for self-assembly with short-range interactions. Proc Natl Acad Sci USA 108(13):5193–5198.

1   D Ak!

[S9]

2 Dk!   F

[S10]

Zeravcic and Brenner www.pnas.org/cgi/content/short/1313601111

3 of 6

Fig. S1. (A) Absolute yield as a function of rescaled temperature T =e, measured from simulations as explained in the text above. (B) Relative yield of polytetrahedron (black disks) and octahedron (red disks) as a function of rescaled temperature T =e. The horizontal lines are referent values obtained from the experiments and partition function calculations.

Fig. S2. Absolute yield as a function of rescaled temperature T =e for an octahedron made out of three different particle types. Absolute yield when particles are identical is included as well.

Zeravcic and Brenner www.pnas.org/cgi/content/short/1313601111

4 of 6

Movie S1. Visualization of the self-replication simulation with Ncoll = 16, which shows one-half of the replication cycle shown in Fig. 2A of the main text. The monomers that comprise the bath are not visible unless they attach to the already formed octahedra and catalysts. The particle color scheme is as in Fig. 2A of the main text.

Movie S1

Zeravcic and Brenner www.pnas.org/cgi/content/short/1313601111

5 of 6

Movie S2. Visualization of the self-replication simulation with Ncoll = 512, which shows multiple replication cycles. The data and snapshots in Fig. 3 of the main text are taken from this movie. The monomers that comprise the bath are not visible unless they attach to the already formed octahedra and catalysts. The particle color scheme is as in Fig. 2A of the main text.

Movie S2

Movie S3. A demonstration movie of the replication process shown in Fig. 1A of the main text, using magnetic sticks and balls. Movie S3

Zeravcic and Brenner www.pnas.org/cgi/content/short/1313601111

6 of 6