The role of collective motion in examples of coarsening and self-assembly

The role of collective motion in examples of coarsening and self-assembly

Stephen Whitelam    Edward H. Feng    Michael F. Hagan    Phillip L. Geissler Systems Biology Centre, University of Warwick, Coventry CV4 7AL, UK
Department of Chemistry, University of California at Berkeley, Berkeley, CA 94720, USA
Physical Biosciences and Materials Sciences Divisions, Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA
Department of Physics, Brandeis University, Waltham, MA, USA
July 5, 2019

The simplest prescription for building a patterned structure from its constituents is to add particles, one at a time, to an appropriate template. However, self-organizing 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 off-lattice 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: self-assembly and collective motion

Self-assembly refers to the generation of patterns or aggregates through the interaction of autonomous components Whitesides and Grzybowski (2002). The supramolecular self-assembly 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 Fraenkel-Conrat and Williams (1955). Biological components can also self-assemble when removed from their natural environment: protein complexes called chaperonins, for example, form large-scale sheet-like Paavola et al. (2006) and string-like Trent et al. (1997) structures in vitro. Self-assembly is widespread in non-biological contexts, from the myriad patterns formed by soap and oil in water to the ribbon-like 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.



Figure 1: The competitive nature of self-assembly: an illustration of the antagonism between the requirements for structural stability and kinetic accessibility. We consider a schematic system of ovoid nanoparticles (mimicking, for example, inorganic nanorods), equipped with pairwise end-to-end interactions of strength (darker shading of rod ends denotes greater interaction strength) and specificity, or inverse angular tolerance, (inversely proportional to shaded area). Stable structures of the required symmetry in general require strong, specific interactions (top left). However, such requirements tend to frustrate the kinetics of assembly (top right): too strong an interaction and structures fail to relax as they grow; too specific an interaction and productive binding events are rare. The competition between kinetics and thermodynamics dictates the region of viable assembly (bottom).

The success of self-assembly 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 self-assembly 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 inter-component 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, self-assembling components generally possess anisotropic interactions of maximal strength much greater than and range much less than a particle diameter. Strong, short-ranged interactions place stringent limits on the maximum integration time step able to preserve numerical stability. Under such conditions, simulations of large-scale assemblies are very time-consuming, 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 coarse-grained 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 single-particle 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 ‘virtual-move’ 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 single-particle 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 self-assembly 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 self-assembly of the homogeneous phase. We find a range of temperatures at which inter-particle 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 single-particle 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 three-dimensional off-lattice system of hard spheres with isotropic pairwise square-well interactions, a model of strongly-associating 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. Inter-particle forces are sufficiently strong that substantial correlated motion emerges on the fundamental timescale , but the geometry of inter-particle 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 ‘virtual-move’ Monte Carlo Algorithm



Figure 2: Illustration of the VMMC procedure applied to a collection of pairwise-interacting oval nanoparticles. (a) Starting in microstate we pick a seed particle (shaded) from which to ‘grow’ a pseudocluster, together with a virtual move map (denoted by the arrow). The pseudocluster is a set of linked particles that will experience a trial move. The virtual move map is both a device for computing neighboring potential energy gradients, and defines the move that the pseudocluster will execute. The pseudocluster growth procedure is an iterative linking scheme akin to the Swendsen-Wang (SW) algorithm Swendsen and Wang (1987): we propose links between the seed and all particles with which it interacts, and continue iteratively until links have been proposed between the pseudocluster and all particles with which its constituents interact. Our procedure differs from the SW algorithm in that here links are formed not on the basis of pairwise energies, but instead on the basis of pairwise energy gradients. These gradients are computed by executing virtual moves of shaded particles (example shown in lower panel), and recording neighboring pairwise energies before and after those moves. Conditioning the acceptance criterion upon reverse virtual moves (not shown) ensures that detailed balance is preserved. The linking criterion is similar to that of the geometric cluster algorithm of Liu and Luijten Liu and Luijten (2004), although here we propose links only between particles that interact in the initial configuration. We seek to effect dynamically realistic local and collective motion, rather than the collective and nonlocal equilibration-speeding motion of Ref. Liu and Luijten (2004). With the pseudocluster so defined (b) we displace it according to the virtual move map, resulting in microstate (c). We then evaluate the Monte Carlo acceptance criterion (see text) and accept or reject the move.

Here we summarize the virtual-move 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


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 link-forming 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


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


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 self-assembly: the attractive lattice gas

In this section we consider the quench-driven 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 solvent-mediated nanoparticle aggregation Rabani et al. (2003). We regard coarsening within the attractive lattice gas as a prototype of self-assembly 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 single-particle Monte Carlo algorithms, inducing transitions between microstates by moving a single particle to an unoccupied nearest-neighbor 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 nearest-neighbor sites. We express temperature in units of . Particles are dispersed randomly on the lattice with concentration and are evolved using virtual-move 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 single-particle 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 two-dimensional 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’ Zinke-Allmang et al. (1992)), . If exceeds , we expect collective modes of motion to significantly influence assembly.

Within the virtual-move algorithm the probability of whole-cluster motion is approximately


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




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 .



Figure 3: Estimate of the size of the largest cluster for which collective motion is at least as probable as the unbinding of its constituent monomers, as a function of particle binding energy . Estimate is derived from a simple argument (Equation (7)) pertaining to the 2 lattice gas, assuming constant coordination number and 3. For the two larger values of , is appreciable even for relatively modest values of .

In we have that , while . For sufficiently large we therefore expect single-particle binding events to dominate. To determine the cluster size at which unbinding and collective motion are equally likely, we equate Equations (4) and (5):


Clusters smaller than will likely move as a whole. To estimate at low temperatures, we ignore the term logarithmic in to obtain


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 .



Figure 4: Estimate of the regime of collective motion for the 2 attractive lattice gas for the given values of cluster size and diffusion exponent . We show the intersection of the ‘clustodal’ (locus of ) with surfaces of constant cluster size as a function of and . For given and , collective motion dominates the behavior of clusters of size (or smaller) in the region below the corresponding line. Hence, collective motion dominates this kinetic phase diagram for a broad range of temperatures and particle concentrations, with quantitative differences evident upon changing the exponent (governing cluster diffusivity) from to 1. We assume constant coordination number .

We next estimate the timescale upon which two clusters of size collide. The cluster diffusion constant is


where is the diffusion constant of a monomer. The timescale upon which two clusters of size encounter each other through diffusion is approximately




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


We compare the cluster-cluster encounter timescale with the timescale required for monomers to unbind from clusters of size and encounter other structures of size :


At low temperature the unbinding timescale is much larger than the timescale for diffusion of a single particle between clusters.


[width= 7.5cm]collective2 \includegraphics[width= 7.5cm]collective

Figure 5: Top row: Kinetic behavior of individual clusters within the attractive lattice gas at different . Starting from disordered configurations at concentration we evolve each system according to virtual-move Monte Carlo dynamics with . We ‘capture’ clusters of size or 30 as they develop, and subject these clusters repeatedly to the same collective-move Monte Carlo procedure, recording (but not making) accepted moves. We plot the probability of the concerted motion of a sub-cluster of size , , normalizing data by setting . At the highest temperature the probability of monomer unbinding, , is rapid relative to that of whole-cluster motion, , allowing assembled structures to relax as they grow. Note however that correlated motions contribute measurably to relaxation dynamics even at high temperature. At lower temperatures, collective motion of the whole cluster predominates. Bottom row: typical 30-member clusters obtained at the three temperatures considered.

[width= 8.5cm]lattice4

Figure 6: Configurations of the lattice gas at fixed time following a quench from a disordered state, evolved using Monte Carlo dynamics with cluster diffusivity and particle concentrations (top panel) and 0.4 (bottom panel). The equilibrium states (obtained using nonlocal moves, panel ‘equil.’) are in all cases phase-separated configurations. Top panel: At high temperature, collective ()- and single-particle () motion give rise to visually similar assembly behavior, while at low temperature collective motion drives the formation of gel-like kinetically trapped structures. Gelation is avoided by suppressing collective motion. Similar behavior is seen at higher concentrations (bottom panel), although collective motion at the lower temperatures is hindered sterically. Coarsening images captured after 2 million MC steps per particle (1.5 million steps for and ).

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 self-assembly 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, single-particle 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 single-particle unbinding for clusters of size less than


which can be very large at low temperature and low to moderate densities. Under these conditions, whole-cluster motion will dominate the system’s assembly dynamics on large length and timescales.



Figure 7: Growth of domain size with Monte Carlo time within the attractive lattice gas, for cluster diffusivity or at particle concentration on a lattice of size . We perform quenches to , , and . Quantitative differences between dynamical protocols are seen even at the highest temperature. Data correspond to the mean of 10 trajectories; error bars are shown sparsely for clarity.

The regimes of single-particle-dominated motion and cluster-dominated 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 virtual-move 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 single-particle unbinding is rapid on the timescale of whole-cluster 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 single-particle algorithms. The non-monotonicity of reflects the fact that moves of sub-clusters 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 single-particle dynamics . The equilibrium states (‘equil.’) are in all cases phase-separated 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 single-particle 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 large-scale 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, , single-particle 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 single-particle motion. The mean cluster size does not become large enough that whole-cluster diffusion is negated as a viable means of assembly. Interestingly, and counter-intuitively, 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 single-particle unbinding events; the two algorithms behave similarly in this regime.

At the two lowest temperatures, and , single-particle 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 fractal-like clusters. Such clusters coarsen chiefly by way of rare single-particle unbinding events, a mechanism that leads to smaller increases in domain length with time than do mutual cluster-cluster collisions.

Our results show that the attractive lattice gas displays a rich variety of coarsening or self-assembly 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 large-lengthscale 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, short-time cooperativity versus long-time 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 off-lattice 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 polymer-mediated 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 phase-separate or gelate.

We consider as a simple model of this system a three-dimensional 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 ‘virtual-move’ 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 large-scale 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 ‘freely-draining’ with ‘Stokesian’ cluster diffusion scalings. Freely-draining 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.


[width= 6cm]colloid_kinetics2

Figure 8: Kinetic data for the self-assembly of model colloids. We show number of clusters as a function of Monte Carlo time for systems A ) and B ), under ‘freely-draining’ (Brownian) and ‘Stokesian’ cluster diffusion scalings. System A assembles or coarsens chiefly through binding, unbinding and diffusion of monomers, and displays little dependence upon cluster diffusivity. These data superpose on data generated using a single-particle Monte Carlo algorithm (SPM). System B instead assembles through the concerted motion of aggregates (when collective motion is permitted), and displays a dependence upon cluster diffusivity. When denied collective motion, clusters coarsen via very slow monomer unbinding events. Data are averaged over 5 stochastic trajectories, each generated using 1058 particles. Error bars are displayed sparsely for clarity.

[width= 6cm]colloid_kinetics3

Figure 9: Kinetic data for the self-assembly of model colloids. We show the number of clusters of given set of sizes , , as a function of Monte Carlo time for system B under Brownian and Stokesian cluster diffusion scalings, and for single-particle moves (SPM). Here , denotes cluster sizes from 7 to 10, and denotes cluster sizes from 26 to 64. Dynamical pathways for the two realizations of collective motion differ even for relatively small aggregates. Data are averaged over 5 stochastic trajectories, each generated using 1058 particles. Error bars are displayed sparsely for clarity.

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 single-particle 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 slowly-coarsening 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 phase-separate or gel. The authors of Ref. Charbonneau and Reichman (2007) discuss possible reasons for the disparity between this experimental observation and the cluster-cluster 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.


[width= 4cm]ep2p6_vmmc \includegraphics[width= 4cm]ep2p6_spm \includegraphics[width= 4.5cm,angle=90]ep12_vmmc_bonds \includegraphics[width= 3.7cm]ep12_spm_bonds

Figure 10: Structures obtained from Monte Carlo simulation of 6500 particles equipped with attractive interactions of type A (top row, excluded-volume view) and B (bottom row, bond view), permitting (left column) or forbidding (right column) explicit collective (Stokesian) motion. The behavior of system A differs under the two dynamical protocols only at late times, where collective motion leads to the aggregation of crystalline clusters. For system B, collective motion induces gelation while single-particle motion results in small, slowly-ripening clusters. These behaviors are similar to those of the lattice gas (see Figure 6). Times of image capture (clockwise from top left), in units of MC steps, are 1.4, 2.9, 1.7 and 0.6.

[width= 2.3cm]ep2p6_vmmc_early \includegraphics[width= 2.3cm]ep2p6_vmmc_mid \includegraphics[width= 2.3cm]ep2p6_vmmc

Figure 11: Assembly pathway for system A under explicit collective moves. Assembly is driven principally by single-particle binding, unbinding and diffusion, with the effects of collective motion apparent only at late times when phase-separated clusters fuse. Times of image capture, in units of MC steps, are from left to right 0.2, 0.6 and 1.4.

V A model of viral capsid self-assembly

The ability of interacting proteins to spontaneously form icosahedral capsids in vivo and under some conditions in vitro is a striking example of biological self-assembly. 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’, 60-member 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 short-ranged, 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 virtual-move Monte Carlo a collection of 1000 subunits. These are initially randomly oriented and dispersed within a three-dimensional 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 inter-unit 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 virtual-move 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 non-monotonic 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 subunit-subunit binding events that are too rare to induce assembly on the timescales simulated. We show in Figure 13 example configurations obtained from a well-assembled and a badly-assembled system.



Figure 12: Final yield of capsid model as a function of potential specificity (left column) or strength (right column). Left panels: yield at subunit concentration for attraction strength 16 (top left) or 14 (bottom left); Right panels: yield at binding specificity for concentration (top right) or (bottom right). Circles denote Brownian dynamics results Hagan and Chandler (2006); squares denote results obtained using collective Monte Carlo dynamics; triangles denote results obtained using single-particle Monte Carlo dynamics at times equal to (dashed lines, right column) those at which collective-move data were sampled, or at later times (if available) when yield had appeared to saturate (dot-dashed lines, left column). Times of data capture for Monte Carlo simulations: left panels, 4.2 MC steps; right panels 9 MC steps. Exceptions are the late-time single-particle data captured at 9 MC steps (top left panel, triangles) and MC steps (bottom left panel, triangles).

Figure 13: Configurations of the capsid model generated using virtual-move Monte Carlo. Left: configuration illustrating high yield obtained at parameter set , , . Right: malformed shell of 76 particles obtained at parameter set , , .

Figure 14: Capsid kinetics at binding energy and concentration evolved using Monte Carlo dynamics with different degrees of collective motion permitted. We show yield and energy per particle as a function of time. Assembly under collective moves is more efficient than under single-particle moves. Restoring collective motion of dimers recovers a substantial fraction of the efficiency of fully collective motion. Data in left panel are the mean of 10 trajectories, data in right panel are the mean of 4 trajectories. Error bars are displayed sparsely for clarity.

We can examine the role of collective motion in this example of self-assembly 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 single-particle Monte Carlo dynamics (triangles) at Monte Carlo times equal to (dotted lines) the times at which collective-move Monte Carlo data were sampled, and for times at which yields appeared to saturate (dot-dashed lines). Saturated yields were not obtained (after 300 hours of simulation) for single-particle 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 higher-order clusters set equal to zero) and for single-particle motion (diffusion constant of dimers and higher-order 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 non-complementarity of model chaperonin-chaperonin 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. Single-particle moves, which suppress such cluster diffusion, give rise to small, isolated, well-formed 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 self-assembling 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 phase-separated. Suppressing collective motion at the same thermodynamic states instead leads to nucleation and growth of clusters. Collective motion plays a similar equilibration-frustrating role within a simple off-lattice model of associating colloids in three dimensions, driving gelation when interactions are strong. By contrast, correlated motions of anisotropically-interacting subunits within a model of viral capsid assembly lead to more efficient self-assembly under all conditions considered. Our results suggest that gelation of homogeneously-interacting particles might be regarded as coarsening in the presence of collective motion, and that the interplay of phase-separation and gelation depends both upon the thermodynamic state, and upon the rate of diffusion of self-assembled aggregates.

Control of collective motion in real systems is very difficult, but might be possible in special circumstances. Individual colloids in high-salt polymer solutions interact via only short-ranged depletion attractions, but buildup of charge on large aggregates can lead to many-body electrostatic repulsions. It is possible that fine-tuning 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. Many-body effects of a different nature can be induced by long-ranged hydrodynamic interactions in, for example, sedimenting colloidal suspensions Padding and Louis (2004). In the case of viral capsids, productive multimer-multimer 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 HHMI-NIBIB Interfaces Initiative grant to Brandeis University, and Brandeis University startup funds. PLG acknowledges funding from the US Department of Energy.

Viii Appendix A

The virtual-move 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 pseudocluster-formation 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).

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

  2. 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 ‘pre-link’ between and (not an automatic link, as described in the main text).

    • If the pre-link 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 pre-link 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 .

  3. 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


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


and so =1. Lastly, the product in the final line of Equation (15) runs over all fully-formed 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


The link-formation and link-failure 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.


  • Whitesides and Grzybowski (2002) G. Whitesides and B. Grzybowski, Science 295, 2418 (2002).
  • Fraenkel-Conrat and Williams (1955) H. Fraenkel-Conrat 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).
  • Zinke-Allmang et al. (1992) M. Zinke-Allmang, 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).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description