# Rotational diffusion affects the dynamical self-assembly pathways of patchy particles

See allHide authors and affiliations

Edited by William Bialek, Princeton University, Princeton, NJ, and approved November 4, 2015 (received for review July 6, 2015)

## Significance

Recent experiments show that the rotational diffusion of proteins and patchy colloids does not always follow the Stokes–Einstein relation, especially in crowded environments or with the use of external fields. Because cellular cytoplasm is very crowded, this finding can have consequences for protein complex formation such as viruses. We study the kinetic network of simple models for proteins and patchy colloids and find that their dynamical self-assembly pathways change with varying the rotational diffusion constant. In particular, lowering the rotational diffusion avoids frustrated intermediate states. Such control of kinetic networks would also benefit the design of new self-assembled functional materials.

## Abstract

Predicting the self-assembly kinetics of particles with anisotropic interactions, such as colloidal patchy particles or proteins with multiple binding sites, is important for the design of novel high-tech materials, as well as for understanding biological systems, e.g., viruses or regulatory networks. Often stochastic in nature, such self-assembly processes are fundamentally governed by rotational and translational diffusion. Whereas the rotational diffusion constant of particles is usually considered to be coupled to the translational diffusion via the Stokes–Einstein relation, in the past decade it has become clear that they can be independently altered by molecular crowding agents or via external fields. Because virus capsids naturally assemble in crowded environments such as the cell cytoplasm but also in aqueous solution in vitro, it is important to investigate how varying the rotational diffusion with respect to transitional diffusion alters the kinetic pathways of self-assembly. Kinetic trapping in malformed or intermediate structures often impedes a direct simulation approach of a kinetic network by dramatically slowing down the relaxation to the designed ground state. However, using recently developed path-sampling techniques, we can sample and analyze the entire self-assembly kinetic network of simple patchy particle systems. For assembly of a designed cluster of patchy particles we find that changing the rotational diffusion does not change the equilibrium constants, but significantly affects the dynamical pathways, and enhances (suppresses) the overall relaxation process and the yield of the target structure, by avoiding (encountering) frustrated states. Besides insight, this finding provides a design principle for improved control of nanoparticle self-assembly.

In nature, self-assembled complex structures and networks often provide function. Prime examples are virus capsids, where capsomer proteins with specific interaction sites self-assemble into various structures, such as icosahedrons and dodecahedrons. Protein complexes can spontaneously form in the living cell, e.g., in signal transduction networks. Self-assembly of small designed building blocks can provide novel (bio)materials with desired properties. Such building blocks can consist of proteins, synthetic polypeptides, but also of colloidal particles. Particularly, the advent of novel synthesis routes for colloidal particles with a valence, so-called “patchy particles” opened up avenues for designing colloidal superstructures. Numerous experimental, theoretical, and numerical studies have enabled understanding of the phase behavior of these particles, predicting not only interesting building blocks for new functional materials, but also demonstrating new physics (1⇓⇓⇓–5). Design principles for colloidal superstructures can predict which structure is the most thermodynamically favorable state (6). However, the fact that kinetics often trumps thermodynamics can hamper such design of colloidal superstructures. Strong directional binding and slow dissociation (7) can kinetically trap patchy-particle systems in a malformed state, rendering it unable to reach the designed equilibrium (ground) state. Controlling the self-assembly of colloidal particles thus requires understanding how the system evolves toward equilibrium, which is dictated by the kinetic network spanning all the states in which the system can occur. Several studies have demonstrated that the assembly toward the final ground state is affected by changing the interaction between patchy particles (8, 9), which affects both the thermodynamics and kinetics of the system. In contrast, here we study how the pathways change upon changing the dynamics only. Colloidal dynamics is often stochastic, and well described by overdamped Langevin (Brownian) dynamics (10). The rotational and translational diffusion constants of anisotropic particles are under standard conditions coupled via the Stokes–Einstein relation. However, in environments with high molecular crowding or in external fields the Stokes–Einstein relation is not necessarily valid anymore (11⇓–13). Depending on the molecular crowder, the ratio between the rotational and translation diffusion constant can go up or down. Whereas crowding agents can in principle also change the binding equilibrium constants, this effect is expected to be weak if the specific interactions are sufficiently strong. Moreover, through the use of external fields, e.g., magnetic or electric, rotational motion could be controlled to a much further extent (14). In this work we investigate how varying the ratio of translational versus rotational diffusion constant influences the equilibrium kinetic network for small self-assembled clusters of colloidal patchy particles and how it could affect the design of self-assembling biological or artificial functional building blocks. Such particles provide a simple model for self-assembly of protein complexes such as in viruses or signal-transduction networks. Understanding and prediction of the colloidal self-assembly mechanisms requires the rates and pathways for all possible dissociation and association events in the kinetic network. However, on the timescale of the dynamics of the microscopic particles, binding and certainly dissociation processes are usually rare events due to high free-energy barriers caused by strong directional binding. As straightforward dynamical simulation is extremely inefficient, we used the Single Replica Transition Interface Sampling (SRTIS) algorithm to collect all possible (un)binding trajectory ensembles relevant to the patchy colloid assembly (15).

Surprisingly, even for the dimerization of a one-patch particle we already find an effect of the rotation on the formation dynamics. Next, we investigate a dimerization of two-patch particles, which exhibits an intermediate state and multiple pathways of formation. However, the effect of the rotation becomes truly important if metastable intermediates are possible, such as in tetrahedron formation of a three-patch particle. Here we find that varying rotational diffusion favors one pathway over the other, without changing the equilibrium constants. Finally, we investigate the entire nine-state kinetic network of a tetrahedron cluster, and show that a change in the rotational diffusion shifts the preferred self-assembly pathways significantly. Whereas for low rotational diffusion the overall rate of tetrahedron formation decreases, frustrated states are avoided, leading to significantly less kinetic trapping.

In principle, this effect results in a better yield of the designed ground state, and in less kinetic trapping. Whereas in our examples the intermediate states are not truly kinetic traps, we envision that for more complex target structures, e.g., virus capsids, which compete against disordered aggregates, the price paid for slowing down the dynamics is outweighed by a preference for the correct self-assembly route and leads to higher relative yields. Experiments on protein complex formation and colloidal assembly should be able to test this prediction, for instance by searching for nonmonotonous behavior upon dilution of the crowding environment. Including the interplay between rotational and translational diffusion in the self-assembly design of new supracolloidal structures opens up new opportunities for improved control of bottom-up synthesis of functional materials. Controlling the dynamics could also lead to a design principle for functional materials where the target structure is not necessarily the thermodynamic ground state. Besides crowding agents, such control could be exerted on the rotational motion of the particles via external fields, such as in ref. 14. Moreover, this work helps in understanding how rotational diffusion influences self-assembly processes in naturally occurring crowded environments such as the cell. As it has been established that a crowded protein environment can decrease the ratio of rotational diffusion over the translation diffusion up to more than a factor of 3 (13), our simulations provide (an additional) explanation why protein complex assembly does not suffer more from kinetic trapping as one would naively expect.

## Results and Discussion

### Dimer of One-Patch Particles.

Two hard spheres of radius *σ* with a single, relatively narrow, attractive patch (*Materials and Methods*), analogous to a binding site on a globular protein, can form a dimer. Although, due to the simplicity of the system, no energetically kinetic trap can exist, it is insightful to understand whether and how rotational diffusion can change the dynamical pathway taken during dissociation or association. We sample the path ensemble for this two-state binding process with SRTIS. The interfaces around stable states are defined by the energy of the system. We construct the free-energy landscape, as well as the (un)binding rate constants (*Materials and Methods*). Fig. 1 shows the free energy as a function of the distance between the particles, *B*) is the free-energy minimum at *ϕ*, indicating particles can freely rotate in the unbound state (*U*) where the interaction vanishes. Between these minima an (entropic) binding free-energy barrier is visible. Naturally, the free-energy landscape does not change with the rotational diffusion as it solely depends on the interactions between the particles. In contrast, varying the rotation diffusion ratio *B* and *U*. Fig. 1 visualizes the projection of these dynamical paths, i.e., the reactive path density *Materials and Methods* for a detailed definition) in the *U* to *B*. Note that the path densities are the average of all reactive diffusive pathways (several individual trajectories are shown in the *SI Appendix*, section III). For fast rotational diffusion the path density at large distance becomes almost uniform. The observed enhancement at *ϕ* at the boundary of the unbound state (Fig. 1). This distribution indicates the probability of dimer orientation at dissociation at *ϕ* for low rotational diffusion, while becoming symmetric for high rotational diffusion. The change in reactive pathways is also visible in the reactive current (Fig. 1), which gives the average velocity at which the particles move along the reactive pathways, and hence does not contain off-pathway excursions. Both the reactive path density and the reactive current reveal that for low

The rate constants are given by *I*, and *Materials and Methods*. The measured fluxes out of the stable states, the crossing probabilities and the rates are given in Table 1. Increasing rotational diffusion by 2 orders of magnitude enhances the absolute rates by a factor of 10. Nevertheless, the ratio between the rates does not change within the error bar, as this is determined by thermodynamics. Note that while the bond energy is *U* is still relatively stable with respect to *B*. Of course, this is also dependent on the width of the attractive patch and the simulation volume (particle concentration).

### Dimer of Two-Patch Particles.

We subsequently performed SRTIS of dimer formation for particles with two patches, where association occurs via a multiple step mechanism. Fig. 2 shows the free energy and the reactive path density as a function of the bond distances

### Tetrahedron Assembly.

To understand the effect of rotational diffusion on a more complex kinetic network, we study the formation of a tetrahedron consisting of particles with three patches. The patch vectors are set at an angle of 60*B* and unbound state *U*, there is now also an intermediate state *I* where the fourth particle “wrongly” binds with two bonds to the trimer (Fig. 3, *Inset*). To escape from this frustrated intermediate state, a bond has to be broken before the complete tetrahedron can form. We study how the kinetic network between these three states changes when the rotational diffusion is varied. The rate calculation is more intricate than for the simple dimer system (*SI Appendix*, section I).

We performed SRTIS and computed the rate matrix *SI Appendix*, section IV, and can be used to compute the time evolution of the population by

Next, we examine the self-assembly of four free particles into a tetrahedron. We define a set of nine states: {Unbound state (*U*), one dimer (*SI Appendix*, section IV). Applying transition path theory (TPT) to these rate matrices results in flux matrices, as well as committors for each state for the transition from *U* to *SI Appendix*, section I for details). Fig. 4 shows a graphical summary of this information for different values of *U* and *SI Appendix*, section V we show that the correctly formed tetrahedron pathway dominates for slow rotational diffusion by analyzing the rate matrices using such strong traps.

In Fig. 4 the committor *i* state *U*, is shown as function of

We note that the overall process always becomes slower when decreasing rotational diffusion, even when scaled with the dwell time in the unbound state (*SI Appendix*, section V). However, instead of keeping the translational diffusion fixed and varying the rotational diffusion constant, we can view our results in terms of a fixed rotational diffusion constant and a varying translational diffusion. Then the timescale for tetrahedron formation becomes shorter for lower rotational diffusion.

## Summary

Our path-sampling simulations demonstrate the proof-of-concept that altering the rotational diffusion can affect the kinetic network and hence the mechanism of self-assembly for simple patchy-particle systems. For the dimer system of one-patch and two-patch particles the mechanism of reactive pathways changes significantly with rotational diffusion, from a translational reactive pathway for slow rotational diffusion, to a rotational reactive pathway for fast rotational diffusion. The pathways do not need to follow the free-energy landscape, and avoid the saddle point when the dynamics are changed. In the assembly of a constraint tetrahedron of four particles decorated with three patches, faster rotational diffusion increases the likelihood of a route via the intermediate state. This effect is even stronger for the full tetrahedron assembly, where the frustrated states and kinetic traps are avoided when rotational diffusion is suppressed. We argue that these results hold beyond the tetramer, and can be seen as a generic effect governed by the different roles that rotation and translation play in the self-assembly. Whereas translational motion is important in particle binding, rotation is crucial for transitions between misaligned frustrated states. We support this argument by analyzing the population dynamics of a simple Markov model that mimics a much larger assembly process (*SI Appendix*, section VI). This analysis clearly shows that slower rotational diffusion results in less frustration and in orders of magnitude higher yields, and generalizes our finding into a genuine principle for controlling and understanding complex self-assembly kinetics.

## Materials and Methods

### Model.

We model the particles as patchy hard spheres with diameter *σ*. For a center-to-center distance,

where *k* controls the range of the interaction, *α* and the vector connecting the center of *α* and the center of the patch on the other particle, and *w* is a parameter that controls the patch width (Fig. 1). Note that the arccosine is needed for the angle calculation. As the patch widths in this work are all small, we approximated the arccosine to speed up the simulation significantly.

### Dynamic Monte Carlo.

In our simulations we used dynamic Monte Carlo (DMC) to evolve the system in time via translational and rotational MC moves. A translation move displaces a randomly chosen particle by a random amount in the interval [*N* translational and rotational MC moves, where *N* is the number of particles. The maximum translational displacement is set to

### SRTIS.

SRTIS requires a set *S* of chosen stable states and predefined interfaces around these stable states. The method samples unbiased dynamical trajectories that leave a particular state by the traditional path sampling shooting move, but in addition imposes a bias that enhances crossing of the interfaces via a Wang–Landau scheme. Special moves for switching between stable states enable efficient sampling of all allowed path ensembles in the entire kinetic network. From these ensembles we can subsequently calculate the free energy of the system, get mechanistic insight, and calculate the entire rate network (15, 20, 21).

### Simulation Settings.

#### One-patch dimer.

Two hard spheres of radius *σ* are each decorated with one patch of range *w*, the potential depth *B*), when the energy of the system, *U*), when the particles are separated more than *B* also when the particles are not properly aligned.

#### Two-patch dimer.

The potential parameters are identical to one patch case, except that the patch width is decreased (

#### Tetramer.

Model parameters for this system are

### Rates.

The rate of the system can be obtained via (22)

where *I* and *I* to the first interface of state *J*, which can be factorized into the crossing probability *SI Appendix*, section I.

### Free Energy from the Path Ensemble.

We collect the full path ensemble for dimerization via SRTIS and from this ensemble we can calculate the free-energy landscape (21):

where *C* is an arbitrary constant, *L* time frames **3** is the reweighted path ensemble with weights

where *j* of state *I* and the minus interface, respectively. The constants *i* and

### Reactive Path Density.

The reactive path density is defined as

The function *k*th time step, and zero otherwise;

### Reactive Path Current.

The reactive path current is (17)

where **q**. This current has the same properties as the reactive path density, except for the fact that dead ends are averaged out, which results in the mean direction a path will follow.

### Nested States and Flux Analysis.

Details on the rate calculation for nested states and the flux analysis using TPT are given in the *SI Appendix*, section I.

## Acknowledgments

This work is part of the research program of the Foundation for Fundamental Research on Matter (FOM), which is part of the Netherlands Organisation for Scientific Research (NWO).

## Footnotes

- ↵
^{1}To whom correspondence should be addressed. Email: p.g.bolhuis{at}uva.nl.

Author contributions: J.G., W.K.K., and P.G.B. designed research; A.C.N. and P.G.B. performed research; A.C.N. analyzed data; and A.C.N., W.K.K., and P.G.B. wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1513210112/-/DCSupplemental.

## References

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵.
- Kraft DJ, et al.

- ↵
- ↵
- ↵.
- Northrup SH,
- Erickson HP

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵.
- Noé F,
- Schütte C,
- Vanden-Eijnden E,
- Reich L,
- Weikl TR

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵

## Citation Manager Formats

## Article Classifications

- Physical Sciences
- Physics

- Biological Sciences
- Biophysics and Computational Biology