The role of collective motion in examples of coarsening and selfassembly
Abstract
The simplest prescription for building a patterned structure from its constituents is to add particles, one at a time, to an appropriate template. However, selforganizing molecular and colloidal systems in nature can evolve in much more hierarchical ways. Specifically, constituents (or clusters of constituents) may aggregate to form clusters (or clusters of clusters) that serve as building blocks for later stages of assembly. Here we evaluate the character and consequences of such collective motion in a set of prototypical assembly processes. We do so using computer simulations in which a system’s capacity for hierarchical dynamics can be controlled systematically. By explicitly allowing or suppressing collective motion, we quantify its effects. We find that coarsening within a two dimensional attractive lattice gas (and an analogous offlattice model in three dimensions) is naturally dominated by collective motion over a broad range of temperatures and densities. Under such circumstances, cluster mobility inhibits the development of uniform coexisting phases, especially when macroscopic segregation is strongly favored by thermodynamics. By contrast, the assembly of model viral capsids is not frustrated but is instead facilitated by collective moves, which promote the orderly binding of intermediates consisting of several monomers.
I Introduction: selfassembly and collective motion
Selfassembly refers to the generation of patterns or aggregates through the interaction of autonomous components Whitesides and Grzybowski (2002). The supramolecular selfassembly of thermodynamically stable structures plays a central role in biology, notably in the hybridization of nucleic acid strands, the organization of lipids to form cell membranes, and the assembly of proteins to form casings for viruses and other genetic material FraenkelConrat and Williams (1955). Biological components can also selfassemble when removed from their natural environment: protein complexes called chaperonins, for example, form largescale sheetlike Paavola et al. (2006) and stringlike Trent et al. (1997) structures in vitro. Selfassembly is widespread in nonbiological contexts, from the myriad patterns formed by soap and oil in water to the ribbonlike structures that assemble from cobalt nanoparticles in solution Puntes et al. (2001). The formation of such structures inspires the design of materials with novel properties, the chief goal of nanotechnology.
The success of selfassembly is determined by an interplay of thermodynamics and dynamics. Components equipped with strong and highly directional interactions may stabilize the thermodynamically preferred structure, but are not guaranteed to assemble spontaneously into such a structure. Overly strong attractions may prevent structural relaxation (by impairing unbinding Licata and Tkachenko (2007); Jack et al. (2007) or ‘reversibility’ Rapaport (2008)), resulting in malformed aggregates. Components that must bind via precise alignment with a neighbor may take a prohibitively long time to do so, resulting in slow structural growth. This competition between the requirements for structural stability and kinetic accessibility in general restricts viable assembly to small regions of parameter space, an idea illustrated in Figure 1.
Identifying and controlling regimes of viable assembly pushes the envelope of current experimental capabilities. In support of such work, computer simulation provides a powerful means of understanding selfassembly and its potential for creating new materials. Simulation permits exhaustive trials of model systems at little expense Hagan and Chandler (2006); Horsch et al. (2005); Zhang et al. (2003); Wilber et al. (2007); Zhang et al. (2005). It can reveal both the nature of intercomponent forces that lead to thermodynamically stable structures, and the dynamics through which such components associate Licata and Tkachenko (2007); Hagan and Chandler (2006); Rapaport (2008). When assessing the assembly properties of a given model it is desirable to evolve that model in order to approximate the dynamics that the corresponding physical system would execute. Molecular dynamics Frenkel and Smit (2002); Rapaport (2004) algorithms evolve components according to Newton’s laws of motion, and so are a natural choice for simulating particle systems. However, selfassembling components generally possess anisotropic interactions of maximal strength much greater than and range much less than a particle diameter. Strong, shortranged interactions place stringent limits on the maximum integration time step able to preserve numerical stability. Under such conditions, simulations of largescale assemblies are very timeconsuming, forcing the simulator to choose between focusing on dynamics on smaller scales or starting from forcibly equilibrated samples.
One way to circumvent this problem is to use a coarsegrained dynamical procedure to move particles according to potential energy gradients without explicitly integrating equations of motion. The Monte Carlo technique provides a flexible framework in which to do so Frenkel and Smit (2002); Binder (1997); Kikuchi et al. (1991); Tiana et al. (2007); Berthier and Kob (2007); Berthier (2007). However, conventional Monte Carlo techniques involve sequential moves of individual particles and so neglect the correlated motion of particles on timescales less than the fundamental discrete time step (the time corresponding to a typical discrete particle displacement). For some systems this neglect of collective motion appears to be unimportant: in Ref. Berthier (2007), for example, molecular dynamics results were reproduced using a singleparticle Monte Carlo protocol. For many systems, however, chiefly those whose constituents possess interactions whose strengths vary strongly with angle or distance, neglecting motion correlated on timescales less than leads to unphysical relaxation, particularly at long times.
To address such problems, ‘cluster’ algorithms have been used extensively to effect correlated or collective motion Swendsen and Wang (1987); Wolff (1989); Wu et al. (1992); Orkoulas and Panagiotopoulos (1999); Amar et al. (1988); Liu and Luijten (2004); Mak (2005); Babu et al. (2008); Krauth (2006). In general, such algorithms identify collections of particles to be moved in concert by recursively ‘linking’ particles according to a set of criteria, such as the pairwise energy or degree of proximity of neighboring particles. One such algorithm, the ‘virtualmove’ Monte Carlo (VMMC) procedure of Ref. Whitelam and Geissler (2007), is designed to effect correlated displacements and rotations according to potential energy gradients or forces experienced under ‘virtual’ moves of neighboring particles. We describe this idea in the following section and in the Appendix. This procedure reduces to singleparticle motion when particles experience small energy changes on timescale , but effects collective motion on arbitrarily large lengthscales when particles experience large changes in energy on this timescale.
In this paper we use VMMC to study the qualitative effect of collective motion on selfassembly by explicitly allowing or suppressing collective moves of particles. In Section III we consider the 2 attractive lattice gas, which coarsens upon a temperature quench via the selfassembly of the homogeneous phase. We find a range of temperatures at which interparticle forces are large enough to encourage assembly but not so large that motion is strongly correlated on timescales less than . In this case assembly is driven by singleparticle binding and unbinding events, and little qualitative effect is observed upon accounting explicitly for correlated motion. However, at low temperatures we find that motion on a timescale is strongly correlated according to potential energy gradients. Neglecting or allowing explicit collective motion under these conditions selects drastically different fates for the system: strongly collective motion impairs assembly via the formation of kinetic traps associated with the binding of large clusters. In the language of coarsening, these kinetic traps represent the arrest of phase separation by gelation. Suppressing collective motion suppresses gelation. We present a simple argument designed to estimate the importance of collective motion in the space of temperature and particle concentration.
In Section IV we examine a threedimensional offlattice system of hard spheres with isotropic pairwise squarewell interactions, a model of stronglyassociating colloids. Here, as for the lattice gas, collective motion tends only to impair assembly (or promote gelation) by inducing awkward binding events between large clusters. However, in Section V, we study a model of viral capsid assembly that displays qualitatively different behavior. Interparticle forces are sufficiently strong that substantial correlated motion emerges on the fundamental timescale , but the geometry of interparticle association is such that these collective motions improve assembly by inducing productive collisions between small intermediates larger than monomers. Suppressing this correlated motion slows assembly, but does not strongly impact the final capsid yield at most thermodynamic states. Collective motion therefore appears to play a qualitatively different role within different model systems: it drives the formation of kinetic traps and facilitates orderly growth.
Ii A ‘virtualmove’ Monte Carlo Algorithm
Here we summarize the virtualmove Monte Carlo algorithm. We consider a dimensional collection of particles equipped with pairwise interactions. This algorithm is a dynamic procedure designed to identify, on the basis of potential energy gradients explored on a fundamental timescale , the extent to which the motion of one particle is correlated with that of its neighbors. If so correlated, this motion is effected with a frequency designed to approximate a physical dynamics. An illustration of the procedure is given in Figure 2. We start in microstate . We define a ‘pseudocluster’ (a set of particles to be moved in concert) by choosing as its first member a ‘seed’ particle . We link the seed to a neighbor with probability , which in general depends on a ‘virtual’ move of that defines a notional new microstate, . Particles linked to members of the pseudocluster join the pseudocluster. We proceed iteratively, until no more members are added to the pseudocluster. We accept the move with probability
(1)  
Here is a factor we impose in order to modulate the diffusivity of pseudoclusters according to size; is the energy of the system in microstate ; is the probability of not linking particles and ; and denotes a particular realization of formed and failed links. The linkforming procedure is aborted in situ if the pseudocluster size exceeds a specified cutoff , the smallest integer larger than . Here is a random variable drawn uniformly from the interval . The subsequent rejection of the move is enforced by the factor . This rejection procedure ensures that all particles experience proposed moves with approximately equal frequency. The products over variables quantify the probabilities of not forming proposed links, internal to and external to (for forward and reverse moves), while the products over variables quantify the probabilities of linking together members of the pseudocluster. We choose to link particles and with a probability
(2)  
that depends on a virtual move (translation or rotation) of relative to . Here is the pairwise energy of the bond in microstate , and is the bond energy following the virtual move of . The factor is unity if and interact in microstate , and zero otherwise; the factor terminates link formation if . Linking particles in this fashion ensures that neighbors exert mutual forces proportional to the gradient of their pairwise energies, with motion ‘linked’ or correlated if is large. Equation (2) implies that the acceptance rate (1) reduces to
(3)  
The label identifies particle pairs that start () in a noninteracting configuration and end () with positive energy of interaction (overlapping), or which start () in an overlapping configuration and end () in a noninteracting one. The final line of Equation (3) accounts explicitly for links and failed links internal to the pseudocluster (note that the factors internal to the pseudocluster were omitted in error in Equations (11) and (13) of Ref. Whitelam and Geissler (2007)). In Appendix A we outline a procedure in which this acceptance rate is simplified by ‘symmetrizing’ link formation at the level of link generation.
Iii A prototype of selfassembly: the attractive lattice gas
In this section we consider the quenchdriven coarsening of the 2 attractive lattice gas with conserved particle number. Lattice gas models, whose thermodynamics can be related to the Ising model Chandler (1987); Huang (1987), are used to caricature a diverse range of physical systems, from binary metallic alloys to solventmediated nanoparticle aggregation Rabani et al. (2003). We regard coarsening within the attractive lattice gas as a prototype of selfassembly in which attractive interactions drive the organization of a homogeneous phase. Our aim is to identify the range of temperatures and particle concentrations where collective motion (motion correlated on a timescale less than the discrete time step ) strongly influences assembly.
Many authors have studied the dynamics of the lattice gas using singleparticle Monte Carlo algorithms, inducing transitions between microstates by moving a single particle to an unoccupied nearestneighbor site. The assumption underlying these studies is that sequential moves of single particles represent a good approximation of the dynamics that the corresponding physical system might execute. This assumption is likely accurate for systems that relax via transport of mass from interfaces of high to low curvature by diffusion through the intervening medium Bray (2002). Scaling arguments and simulations based on this physical picture predict the typical domain size to grow as Bray (2002); Amar et al. (1988).
However, in many settings, the transport of mass between domains by the evaporation and diffusion of constituent monomers is not the only possible mode of relaxation. Indeed, we frequently encounter the concerted motion of domains of one phase in another: witness the relaxation of polymers in solution or the the diffusion of nanocrystal aggregates on graphite Ge and Brus (2001). To model such behavior within the lattice gas, we must explicitly account for motion correlated on timescales of order .
We consider a collection of particles on a simple square lattice of sites. Two particles may not occupy the same site, and interact with binding energy when occupying nearestneighbor sites. We express temperature in units of . Particles are dispersed randomly on the lattice with concentration and are evolved using virtualmove Monte Carlo translations. We enforce a scaling of the diffusion constant of for clusters of size . We focus on the difference between the case , corresponding to singleparticle moves, and , corresponding to diffusion akin to Brownian dynamics.
The simplicity of the model allows us to estimate the values of temperature and particle concentration for which we expect motion correlated on a timescale to be important. We present this argument for the dimensional hypercubic lattice gas. We assume for simplicity that growing clusters are compact (as a twodimensional illustration, consider the left half of Figure 6, top row). We wish to estimate the timescale on which two clusters of size encounter each other as a result of their collective diffusion, , and compare this estimate with the timescale on which such clusters exchange mass via evaporation (Ostwald ‘ripening’ ZinkeAllmang et al. (1992)), . If exceeds , we expect collective modes of motion to significantly influence assembly.
Within the virtualmove algorithm the probability of wholecluster motion is approximately
(4) 
where is the probability of linking two particles following a virtual move, and is the change in energy resulting from separating those particles. We impose a factor of , which ensures that particles suffer attempted moves with approximately equal frequencies, and a factor of , to account for our chosen diffusion constant. The number reflects the efficiency with which the recursive algorithm forms links within a cluster. This number is approximately for the version of the algorithm described in the main text, and is zero for the version of the algorithm described in the appendix. This distinction is unimportant at low temperature, and we shall for simplicity take .
The probability that a single particle breaks away from its host cluster is
(5) 
where
(6) 
is a geometric factor quantifying the likelihood that a chosen particle lies on the surface of a cluster (we expect this approximation to be reasonable for ), and is the typical coordination number of a particle on the surface of a cluster of size .
In we have that , while . For sufficiently large we therefore expect singleparticle binding events to dominate. To determine the cluster size at which unbinding and collective motion are equally likely, we equate Equations (4) and (5):
(7) 
Clusters smaller than will likely move as a whole. To estimate at low temperatures, we ignore the term logarithmic in to obtain
(8) 
In Figure 3 we plot obtained from numerical solution of Equation (7) for , , and three values of . We observe that is large even for relatively modest values of when .
We next estimate the timescale upon which two clusters of size collide. The cluster diffusion constant is
(9) 
where is the diffusion constant of a monomer. The timescale upon which two clusters of size encounter each other through diffusion is approximately
(10) 
where
(11) 
is a measure of the distance separating clusters of size . Here is the monomer radius, and is an exponent measuring the increase in distance between structures due to clustering (our simulations indicate that in ). For brevity we write . Then
(12) 
We compare the clustercluster encounter timescale with the timescale required for monomers to unbind from clusters of size and encounter other structures of size :
(13)  
At low temperature the unbinding timescale is much larger than the timescale for diffusion of a single particle between clusters.
For a cluster of size , Equations (12) and (13) quantify the respective timescales for mass transport by collective motion, and for the unbinding and diffusion of monomers. We view the ratio as a measure of the propensity for collective motion at a given scale . The selfassembly of the homogeneous phase proceeds in stages via the appearance of structures of size . At each stage, we expect collective motion to be important if is greater than unity. At high temperature, singleparticle motion dominates: when . At low temperature, instead, collective motion dominates assembly up to large values of : . For the choices and , we expect collective motion to be a more effective means of mass transport than singleparticle unbinding for clusters of size less than
(14) 
which can be very large at low temperature and low to moderate densities. Under these conditions, wholecluster motion will dominate the system’s assembly dynamics on large length and timescales.
The regimes of singleparticledominated motion and clusterdominated motion are separated by the locus of points , defined for a given by the equation . We call this locus the ‘clustodal’, by analogy with ‘binodal’ (locus of phase equilibrium), and ‘spinodal’ (locus of the onset of spontaneous decomposition).
In Figure 4 we show a plot of the propensity for collective motion, , in the plane, based on our simple argument. We consider two values of (quantifying the cluster diffusion rate). The intersections of the clustodal with surfaces of constant indicate where in state space collective motion dominates assembly dynamics on a scale . Thus for a thermodynamic state , say, our estimate suggests that assembly dynamics (for ) for clusters of size 30 (but not size 300) is dominated by motion correlated on times less than the discrete time step .
In simulations of the lattice gas we indeed find that collective motion is important for a broad range of temperatures and particle concentrations. In Figure 5 we offer a microscopic perspective on the meaning of the clustodal by analyzing the (averaged) motion of individual lattice gas clusters of sizes and 30 at three temperatures. Starting from randomly dispersed monomers at concentration , we capture clusters of size as they assemble. We use the virtualmove algorithm (enforcing a diffusion constant scaling ) to calculate the likelihood of the correlated motion of a subset of size of a given cluster, normalized by the likelihood of monomer unbinding from that cluster. At high temperature singleparticle unbinding is rapid on the timescale of wholecluster motion, while at low temperature the converse is true. However, even at the highest temperature we find that motions correlated on timescales less then contribute appreciably to the ‘spectrum’ of dynamic relaxation. Such motion is ignored within singleparticle algorithms. The nonmonotonicity of reflects the fact that moves of subclusters of size involve the breaking of many bonds, and so are amongst the least favorable processes in an energetic sense; the asymmetry of results from the fact that cluster diffusion constants decrease with increasing . Each data set in Figure 5 was obtained by applying moves to a cluster of size (recording any accepted move but not making that move), and averaging over such clusters.
In Figure 6 we show lattice gas configurations at fixed time following a quench from high temperature to a low final temperature, for four values of this final temperature and for two particle concentrations (here and in subsequent sections we used VMD Humphrey et al. (1996) to render simulation configurations). We employ both collective dynamics and singleparticle dynamics . The equilibrium states (‘equil.’) are in all cases phaseseparated configurations, which we deduce by using nonlocal Monte Carlo moves Chen and Siepmann (2001). At the highest temperatures, coarsening dynamics are visually similar. In this regime the system lies ‘outside’ the clustodal for all but the smallest values of , and singleparticle evaporation and diffusion is the dominant mode of relaxation. At lower temperatures and the lower of the two particle concentrations, striking differences emerge: for sufficiently low temperatures collective modes of motion () induce kinetic frustration through the mutual collisions of clusters that bind awkwardly and fail to relax before encountering similar structures. If collective motion is suppressed then instead isolated, compact structures are formed. At the higher particle concentration, steric effects partially frustrate largescale collective motion, and coarsening patterns are qualitatively more similar (although differences can be clearly seen).
We quantify the influence of collective motion upon assembly dynamics by measuring the characteristic lengthscale of domains as a function of time at , for and . We calculate the domain lengthscale from the first moment of the structure factor Amar et al. (1988); we display in Figure 7 results for four temperatures. At the highest temperature, , singleparticle and cluster algorithms show the same qualitative behavior, with both and increasing monotonically with time. Both dynamics produce compact clusters, as may be seen in Figure 6. However, even at this high temperature we observe quantitative differences between the two algorithms, with collective motion giving rise to larger domains at any given time than does singleparticle motion. The mean cluster size does not become large enough that wholecluster diffusion is negated as a viable means of assembly. Interestingly, and counterintuitively, the growth in lenghscale under the two algorithms is more similar at late times at a slightly lower temperature of . Here we observe the same qualitative assembly behavior under the two algorithms as at , but now the mean cluster size is larger at late times because of the stronger thermodynamic driving force. Even when collective dynamics is permitted, collisions involving these larger clusters are sufficiently rare that increases in are driven chiefly by singleparticle unbinding events; the two algorithms behave similarly in this regime.
At the two lowest temperatures, and , singleparticle dynamics displays the same qualitative behavior as it does at and , but cluster moves show different behavior. While still increases monotonically with time, the slope has at late times a constant value at , and decreases at . The corresponding pictures for in Figure 6 suggest that this decrease in slope results from the steric hinderance associated with percolating fractallike clusters. Such clusters coarsen chiefly by way of rare singleparticle unbinding events, a mechanism that leads to smaller increases in domain length with time than do mutual clustercluster collisions.
Our results show that the attractive lattice gas displays a rich variety of coarsening or selfassembly behaviors, and that these behaviors depend sensitively upon the dynamical protocol used to evolve the system. Figure 6 demonstrates that when one accounts for collective motion , nucleation and growth mechanisms () yield to gelation as temperature is reduced ( and ). Gelation at low temperature is driven by explicit collective motion, or motion correlated on the fundamental timescale : when we forbid such motion , gelation is avoided in favor of very slow nucleation and growth. It is interesting that collective motion appears to be crucial to the formation of largelengthscale gels within the attractive lattice gas at low particle concentration, but that such explicit correlated motion is not required to observe glassy behavior within a model of (dense) silica Berthier (2007). These results imply a difference in the nature of the dynamical cooperativity associated with glasses and attractive gels, respectively, shorttime cooperativity versus longtime cooperativity.
Iv The attractive lattice gas generalized to continuous space: model associating colloids
The effects of collective motion seen in the lattice gas are also apparent in a simple offlattice model of hard spheres equipped with isotropic interactions. Such a model is a caricature of colloidal particles in solution with small polymers; the colloids associate by virtue of a polymermediated depletion attraction. In Ref. Lu et al. (2006), attractive colloidal particles of this nature were observed in experiment to assemble into clusters whose geometry depended upon the range and strength of the depletion attraction. Strikingly, at low colloid concentration these clusters appeared to be stable, defying the expectation that components with strong pairwise attractions should phaseseparate or gelate.
We consider as a simple model of this system a threedimensional collection of hard spheres of diameter equipped with an attractive square well of range and strength . Spheres are placed randomly within the simulation box and occupy 4% of its volume. We impose periodic boundary conditions. We consider two of the parameter sets described in Ref. Lu et al. (2006): a potential of moderate strength and range, A ), and a potential of considerable strength and short range, B ). We carried out Monte Carlo simulations of systems of and 6500 particles, drawing particle displacement magnitudes (in units of ) uniformly from the interval , and rotation angles uniformly from a distribution with maximum . We used ‘virtualmove’ translations (employing the algorithm described in the appendix), and used a static linking scheme to effect rotations about the center of mass of a chosen pseudocluster. These choices of displacement magnitude and rotation angle imply a basic timescale sufficiently large that to a good approximation largescale cooperative motion cannot occur from uncorrelated moves of single particles; such moves effect only local structural relaxation and the binding, unbinding and diffusion of monomers. Our aim is to compare such dynamics with pathways accessible to explicitly correlated motion. When we consider cluster moves we also compare ‘freelydraining’ with ‘Stokesian’ cluster diffusion scalings. Freelydraining motion implies , where is the number of monomers comprising the cluster and is the cluster moment of inertia about the rotation axis. For Stokesian scalings we take , with a measure of the cluster radius of gyration perpendicular to the translation vector or axis of rotation.
As in the lattice gas, collective motion plays a dominant role when potential energy gradients encountered on a timescale are large, and is less important than singleparticle unbinding events for small energy gradients. In Figures 8 and 9 we show kinetic measures for systems of 1058 particles. System A displays no statistically significant difference in the evolution of the number of clusters with time if collective motion is allowed or suppressed. System B, however, experiences dramatically different fates under collective and individual particle motion, with the former inducing gelation and the latter resulting in a slowlycoarsening collection of isolated clusters.
In Figures 10 and 11 we display snapshots of these systems for collections of 6500 particles (in these images system A was evolved using a maximum translation of 0.3). System undergoes phase separation into crystalline clusters. The effect of collective motion is apparent only at late times when clusters fuse; such fusing is in qualitative agreement with Brownian dynamics simulations of a similar system Charbonneau and Reichman (2007). System B forms under collective motion stringy, frustrated aggregates that merge and form a gel, with the timescale for gelation different for the two different cluster diffusivities modeled. In experiment, colloids with a strong depletion attraction (of nature similar to systems A and B considered here) form instead isolated clusters Lu et al. (2006), and do not phaseseparate or gel. The authors of Ref. Charbonneau and Reichman (2007) discuss possible reasons for the disparity between this experimental observation and the clustercluster aggregation seen in simulations in Ref. Charbonneau and Reichman (2007) (similar to those observed here): one such suggestion is that, in experiment, an accumulation of charge on large bodies might induce a repulsion that stabilizes a phase of isolated clusters.
V A model of viral capsid selfassembly
The ability of interacting proteins to spontaneously form icosahedral capsids in vivo and under some conditions in vitro is a striking example of biological selfassembly. Protein subunits assemble under a variety of conditions, avoiding both kinetic and thermodynamic traps. Understanding the mechanisms that render assembly so robust is an essential step towards designing synthetic analogs of viral capsids. In addition, a comprehensive understanding of the viral capsid formation process will spur the design of antiviral drugs and new drug delivery systems; the latter could possess the ability to assemble and disassemble around their cargo without requiring explicit external control.
In Ref. Hagan and Chandler (2006) a class of simple models was introduced in order to study the mechanism by which interacting protein subunits might form ‘capsids’, 60member closed shells having icosahedral symmetry. We focus here on the ‘B3’ model of that reference. Units interact via a pairwise potential that models an excluded volume and a shortranged, angularly specific attractive interaction designed to stabilize capsids.
Here we examine the role of collective motion within the assembly dynamics of this model. We evolved via virtualmove Monte Carlo a collection of 1000 subunits. These are initially randomly oriented and dispersed within a threedimensional simulation box; we take periodic boundaries in each dimension. We scaled cluster diffusivities in order to approximate Brownian dynamics. As found in Reference Hagan and Chandler (2006), we observe a regime of interunit potential strength and specificity within which assembly is robust. We show in Figure 12 the capsid yield (fraction of units residing in complete capsids) at fixed observation time as a function of potential strength or specificity . Yield data obtained via virtualmove Monte Carlo agree with Brownian dynamics results to within statistical error (considerable variations in yield are observed within each algorithm at a given thermodynamic state). Particle displacement magnitudes are drawn from a uniform distribution with maximum equal to a length unit ; particle interaction range is . Rotations are scaled accordingly. Yields are nonmonotonic functions of these parameters, for the reasons outlined in Figure 1: overly strong or insufficiently specific interactions promote malformed intermediates that fail to assemble into complete structures; overly weak or specific interactions result in productive subunitsubunit binding events that are too rare to induce assembly on the timescales simulated. We show in Figure 13 example configurations obtained from a wellassembled and a badlyassembled system.
We can examine the role of collective motion in this example of selfassembly by explicitly restricting or forbidding collective modes of relaxation within Monte Carlo dynamics. We find that in many cases final yields are not strongly affected by doing so, but that assembly dynamics are impaired. In Figure 12 we show yields obtained using singleparticle Monte Carlo dynamics (triangles) at Monte Carlo times equal to (dotted lines) the times at which collectivemove Monte Carlo data were sampled, and for times at which yields appeared to saturate (dotdashed lines). Saturated yields were not obtained (after 300 hours of simulation) for singleparticle moves in the right panels of Figure 12.
In Figure 14 we present measures of assembly kinetics obtained at two parameter sets for fully collective motion (we move clusters according to the algorithm described in Section II, with the diffusivity of mers chosen to be in order to approximate Brownian dynamics), for motion allowing explicit moves of monomers and dimers only (as for fully collective motion, but with the diffusion constant of trimers and higherorder clusters set equal to zero) and for singleparticle motion (diffusion constant of dimers and higherorder structures set equal to zero). To generate these data we drew particle displacement magnitudes from a uniform distribution with maximum equal to . While at these thermodynamic states the final yields do not depend strongly upon the availability of correlated moves, the kinetics of assembly is more rapid when collective motion is allowed. Interestingly, restoring explicit moves of dimers alone is sufficient to recover much of the ease of assembly afforded by a ‘full’ spectum of correlated motions. These results support the observations made in Ref. Hagan and Chandler (2006). There it was found that B3 capsids grow in part through events consisting of collisions between intermediates larger than monomers, implying that collective motion plays an important role in the model’s assembly dynamics (see also Ref. Zhang and Schwartz (2006)). Further, the most frequent intermediate binding events for the B3 model involved dimers. Indeed, we find here that explicit dimer motion renders assembly much more facile than if moves are uncorrelated. We note also that the tendency of subunits to form closed shells suppresses in large part the kinetic trapping seen in the other models studied in this paper, where large aggregates bind in an awkward fashion and frustrate equilibration.
These capsid assembly results may be contrasted with the assembly properties of a model of interacting protein complexes called chaperonins, studied in Ref. Whitelam and Geissler (2007). There the noncomplementarity of model chaperoninchaperonin interactions coupled with collective modes of motion induce a high degree of kinetic frustration in some regions of the phase diagram: large clusters bind awkwardly, which slows or prevents equilibration. Singleparticle moves, which suppress such cluster diffusion, give rise to small, isolated, wellformed structures. We conclude that particle interaction geometry and the extent of organized aggregates play a decisive role in shaping the effects of collective motion in selfassembling systems.
Vi Conclusions
We have examined the role of collective motion in Monte Carlo simulations of three model systems. We find that collective motion is responsible for gelation at low particle concentrations within the two dimensional lattice gas, at thermodynamic states for which the equilibrium configurations are phaseseparated. Suppressing collective motion at the same thermodynamic states instead leads to nucleation and growth of clusters. Collective motion plays a similar equilibrationfrustrating role within a simple offlattice model of associating colloids in three dimensions, driving gelation when interactions are strong. By contrast, correlated motions of anisotropicallyinteracting subunits within a model of viral capsid assembly lead to more efficient selfassembly under all conditions considered. Our results suggest that gelation of homogeneouslyinteracting particles might be regarded as coarsening in the presence of collective motion, and that the interplay of phaseseparation and gelation depends both upon the thermodynamic state, and upon the rate of diffusion of selfassembled aggregates.
Control of collective motion in real systems is very difficult, but might be possible in special circumstances. Individual colloids in highsalt polymer solutions interact via only shortranged depletion attractions, but buildup of charge on large aggregates can lead to manybody electrostatic repulsions. It is possible that finetuning of solution conditions could be used to select the smallest aggregate lengthscale for which repulsions become significant, thereby controlling the extent to which multimers collide and bind. Manybody effects of a different nature can be induced by longranged hydrodynamic interactions in, for example, sedimenting colloidal suspensions Padding and Louis (2004). In the case of viral capsids, productive multimermultimer bindings occur when aggregates’ exposed contacts meet; changing the number of contacts per subunit (by mutation, for instance), thereby changing the sticky surface presented by a capsid of a given size, might be used to affect the importance of multimer binding to assembly pathways.
Vii Acknowledgements
We thank Ludovic Berthier, Patrick Charbonneau and Thomas Ouldridge for correspondence. SW was supported initially by the US Department of Energy and subsequently by the BioSim European Union Network of Excellence, and acknowledges a Royal Society conference grant that made possible a collaborative visit. Computing facilities were provided in part by the Centre for Scientific Computing at the University of Warwick with support from the Science Research Investment Fund. EHF thanks the Miller Institute of Basic Research in Science for financial support. MFH acknowledges funding from the HHMINIBIB Interfaces Initiative grant to Brandeis University, and Brandeis University startup funds. PLG acknowledges funding from the US Department of Energy.
Viii Appendix A
The virtualmove scheme discussed in Ref Whitelam and Geissler (2007) and the main text accounts for the asymmetry of link formation at the level of the acceptance rate. Here we present a modification of this scheme in which link formation is forbidden if the corresponding link would not form under the reverse move. We find that the acceptance rate for collective motion is made simpler. Consider first a generic dynamic pseudoclusterformation procedure, moving from state to state , in which a link is formed between particles and with probability . This probability is computed by making a virtual move of . We subsequently specialize to the particular choice of given in Equation (2).

Start in state . Choose a seed particle, say , and a move map. Add to the pseudocluster, the list of particles to be moved.

Choose a neighbor of not in the pseudocluster, and with which has not in the current move proposed a link. Call this neighbor . With probability (computed by moving under its virtual map) form a ‘prelink’ between and (not an automatic link, as described in the main text).

If the prelink does not form, we consider that a link has failed to form (outright failure). Choose another neighbor of , say , and return to stage 2, with the replacement .

If the prelink forms, compute the reverse link balance factor .

With probability , form a full link between and . Add to the pseudocluster. Go to particle and proceed from step 2, with the replacement and .

If a full link fails to form, record the link as frustrated. Do not add to the pseudocluster. Choose another neighbor of , say , and return to stage 2, with the replacement .



Proceed until no more links remain to be tested, and evaluate the acceptance probability for the move.
The acceptance probability for the revised algorithm follows from a modification of Equation (1) in the main text. We choose to balance the total rates for forward and reverse moves involving a given realization of internal pseudocluster links, and a realization of internal failed links ‘blind’ to the nature of those failed links (whether frustrated or outright failed). The acceptance rate for the dynamic linking procedure described here, for a generic choice of , is
(15)  
The first two lines of Equation (15) are as Equation (1) of the main text. Variables denote outright failed links between and its environment. When we take as in Equation (2) of the main text, such variables cancel the Boltzmann bond weights for all but a specialized class of moves (see main text). Variables denote frustrated links between and its environment. Frustrated links cannot form between the pseudocluster and its environment for both forward and reverse moves, and should such links form during the forward move then that move must be rejected. Variables denote unformed links internal to the pseudocluster, whether frustrated or outright rejected. We have that
(16)  
and so =1. Lastly, the product in the final line of Equation (15) runs over all fullyformed links. By construction of the linking procedure each quotient in this product is unity: we ensure that links formed during the forward move can also form during the reverse move.
If we take as in Equation (2) of the main text, than the acceptance rate for the collective move is
(17)  
The linkformation and linkfailure factors internal to the pseudocluster have canceled. The factor in the third line of Equation (17) is unity if no frustrated links join the pseudocluster to its environment, and zero otherwise. This factor is required because a frustrated link indicates that a move of one particle relative to another has effectively been rejected. If we then do not form the link, in order to see if both particles are incorporated into the pseudocluster, we must reject any move in which both particles do not end up in the pseudocluster. The advantage of this scheme relative to that presented in the main text is that here links internal to the pseudocluster that form under the forward move, but do so with zero probability under the reverse move, do not automatically result in the rejection of that move.
References
 Whitesides and Grzybowski (2002) G. Whitesides and B. Grzybowski, Science 295, 2418 (2002).
 FraenkelConrat and Williams (1955) H. FraenkelConrat and R. Williams, Proceedings of the National Academy of Sciences of the United States of America 41, 690 (1955).
 Paavola et al. (2006) C. Paavola, S. Chan, Y. Li, K. Mazzarella, R. McMillan, and J. Trent, Nanotechnology 17, 1171 (2006).
 Trent et al. (1997) J. Trent, H. Kagawa, T. Yaoi, E. Olle, and N. Zaluzec, Proceedings of the National Academy of Sciences 94, 5383 (1997).
 Puntes et al. (2001) V. Puntes, K. Krishnan, and A. Alivisatos, Science 291, 2115 (2001).
 Licata and Tkachenko (2007) N. Licata and A. Tkachenko, Physical Review E 76, 41405 (2007).
 Jack et al. (2007) R. Jack, M. Hagan, and D. Chandler, Physical Review E 76, 21119 (2007).
 Rapaport (2008) D. Rapaport, Arxiv preprint arXiv:0803.0115 (2008).
 Hagan and Chandler (2006) M. Hagan and D. Chandler, Biophysical Journal 91, 42 (2006).
 Horsch et al. (2005) M. Horsch, Z. Zhang, and S. Glotzer, Physical Review Letters 95, 56105 (2005).
 Zhang et al. (2003) Z. Zhang, M. Horsch, M. Lamm, and S. Glotzer, Nano Letters 3, 1341 (2003).
 Wilber et al. (2007) A. Wilber, J. Doye, A. Louis, E. Noya, M. Miller, and P. Wong, The Journal of Chemical Physics 127, 085106 (2007).
 Zhang et al. (2005) Z. Zhang, A. Keys, T. Chen, and S. Glotzer, Langmuir 21, 11547 (2005).
 Frenkel and Smit (2002) D. Frenkel and B. Smit, Understanding Molecular Simulation: From Algorithms to Applications (Academic Press, 2002).
 Rapaport (2004) D. Rapaport, The Art of Molecular Dynamics Simulation (Cambridge University Press, 2004).
 Binder (1997) K. Binder, Reports on Progress in Physics 60, 487 (1997).
 Kikuchi et al. (1991) K. Kikuchi, M. Yoshida, T. Maekawa, and H. Watanabe, Chemical Physics Letters 185, 335 (1991).
 Tiana et al. (2007) G. Tiana, L. Sutto, and R. Broglia, Physica A: Statistical Mechanics and its Applications 380, 241 (2007).
 Berthier and Kob (2007) L. Berthier and W. Kob, Journal of Physics, Condensed Matter 19, 205130 (2007).
 Berthier (2007) L. Berthier, Physical Review E 76, 11507 (2007).
 Swendsen and Wang (1987) R. Swendsen and J. Wang, Physical Review Letters 58, 86 (1987).
 Wolff (1989) U. Wolff, Physical Review Letters 62, 361 (1989).
 Wu et al. (1992) D. Wu, D. Chandler, and B. Smit, The Journal of Physical Chemistry 96, 4077 (1992).
 Orkoulas and Panagiotopoulos (1999) G. Orkoulas and A. Panagiotopoulos, The Journal of Chemical Physics 110, 1581 (1999).
 Amar et al. (1988) J. Amar, F. Sullivan, and R. Mountain, Physical Review B 37, 196 (1988).
 Liu and Luijten (2004) J. Liu and E. Luijten, Physical Review Letters 92, 35504 (2004).
 Mak (2005) C. Mak, The Journal of Chemical Physics 122, 214110 (2005).
 Babu et al. (2008) S. Babu, J. Gimel, and T. Nicolai, Arxiv preprint arXiv:0801.4447 (2008).
 Krauth (2006) W. Krauth, Statistical mechanics: algorithms and computations (Oxford University Press, Oxford, 2006).
 Whitelam and Geissler (2007) S. Whitelam and P. Geissler, The Journal of Chemical Physics 127, 154101 (2007).
 Chandler (1987) D. Chandler, Introduction to modern statistical mechanics (Oxford University Press New York, 1987).
 Huang (1987) K. Huang, Statistical Mechanics, J (Wiley, New York, 1987).
 Rabani et al. (2003) E. Rabani, D. Reichman, P. Geissler, L. Brus, et al., Nature 426, 271 (2003).
 Bray (2002) A. Bray, Advances in Physics 51, 481 (2002).
 Ge and Brus (2001) G. Ge and L. Brus, Nano Lett 1 (2001).
 ZinkeAllmang et al. (1992) M. ZinkeAllmang, L. Feldman, and M. Grabow, Surface Science Reports(The Netherlands) 16, 377 (1992).
 Humphrey et al. (1996) W. Humphrey, A. Dalke, and K. Schulten, Journal of Molecular Graphics 14, 33 (1996).
 Chen and Siepmann (2001) B. Chen and J. Siepmann, J. Phys. Chem. B 105, 11275 (2001).
 Lu et al. (2006) P. Lu, J. Conrad, H. Wyss, A. Schofield, and D. Weitz, Physical Review Letters 96, 28306 (2006).
 Charbonneau and Reichman (2007) P. Charbonneau and D. Reichman, Physical Review E 75, 11507 (2007).
 Zhang and Schwartz (2006) T. Q. Zhang and R. Schwartz, Biophys. J. 90, 57 (2006).
 Padding and Louis (2004) J. Padding and A. Louis, Physical Review Letters 93, 220601 (2004).