1887
banner image
No data available.
Please log in to see this content.
You have no subscription access to this content.
No metrics data to plot.
The attempt to load metrics for this article has failed.
The attempt to plot a graph for these metrics has failed.
The full text of this article is not currently available.

F

Studying protein assembly with reversible Brownian dynamics of patchy particles
Rent:
Rent this article for
Access full text Article

Abstract

Assembly of protein complexes like virus shells, the centriole, the nuclear pore complex, or the actin cytoskeleton is strongly determined by their spatial structure. Moreover, it is becoming increasingly clear that the reversible nature of protein assembly is also an essential element for their biological function. Here we introduce a computational approach for the Brownian dynamics of patchy particles with anisotropic assemblies and fully reversible reactions. Different particles stochastically associate and dissociate with microscopic reaction rates depending on their relative spatial positions. The translational and rotational diffusive properties of all protein complexes are evaluated on-the-fly. Because we focus on reversible assembly, we introduce a scheme which ensures detailed balance for patchy particles. We then show how the macroscopic rates follow from the microscopic ones. As an instructive example, we study the assembly of a pentameric ring structure, for which we find excellent agreement between simulation results and a macroscopic kinetic description without any adjustable parameters. This demonstrates that our approach correctly accounts for both the diffusive and reactive processes involved in protein assembly.

Key Topics

Proteins
75.0
Bond formation
33.0
Biochemical reactions
32.0
Enzyme kinetics
32.0
Dissociation
31.0
0
75.0

Studying protein assembly with reversible Brownian dynamics of patchy particles

    Heinrich C. R. Klein1; Ulrich S. Schwarz1,2a)
    1: 1Institute for Theoretical Physics, Heidelberg University, 69120 Heidelberg, Germany; 2: 2BioQuant, Heidelberg University, 69120 Heidelberg, Germany
a)

a)Electronic mail: ulrich.schwarz@bioquant.uni-heidelberg.de

 

Assembly of protein complexes like virus shells, the centriole, the nuclear pore complex, or the actin cytoskeleton is strongly determined by their spatial structure. Moreover, it is becoming increasingly clear that the reversible nature of protein assembly is also an essential element for their biological function. Here we introduce a computational approach for the Brownian dynamics of patchy particles with anisotropic assemblies and fully reversible reactions. Different particles stochastically associate and dissociate with microscopic reaction rates depending on their relative spatial positions. The translational and rotational diffusive properties of all protein complexes are evaluated on-the-fly. Because we focus on reversible assembly, we introduce a scheme which ensures detailed balance for patchy particles. We then show how the macroscopic rates follow from the microscopic ones. As an instructive example, we study the assembly of a pentameric ring structure, for which we find excellent agreement between simulation results and a macroscopic kinetic description without any adjustable parameters. This demonstrates that our approach correctly accounts for both the diffusive and reactive processes involved in protein assembly.

Keywords
    Keywords
    • biochemistry,
    • biodiffusion,
    • Brownian motion,
    • molecular biophysics,
    • proteins,
    • stochastic processes
    Keywords
      Keywords
         

        Assembly of biomolecules into supramolecular complexes is at the heart of many biological processes and the dynamic interplay of the different components leads to biological functionality. The most important type of self-assembling biomolecules are proteins. Although they often have a globular shape, protein assemblies can be strongly anisotropic due to the localized binding interactions between the proteins. Like biological systems in general, protein assemblies are operative on many different scales. Their size ranges from the nanometer-scale (for example, for two-component complexes like Barnase and Barstar1) through tens of nanometers (for example, for viruses, which typically consist of small multiples of 60 components2,3) up to the micrometer scale (for example, the mitotic spindle4–6). Even in steady state most biological complexes remain highly dynamic, with association events being balanced by dissociation events. Prominent examples are nuclear pore complexes7–9 or focal adhesions.10–12 Another important example for the dynamic nature of biological complexes are actin filaments,13–15 which are called due to their continuous exchange dynamics.16 Similar to the broad range of length scales, association rate constants observed in biological systems span a wide spectrum ranging from <103 M−1 s−1 up to >109 M−1 s−1.17 The highly dynamic nature of biological assemblies as well as the large range of involved spatial and temporal scales renders this physical problem challenging but fascinating.

        Interestingly, recent advances in the fabrication of functionalized colloids with directional interactions () make it possible to design elementary building blocks of micrometer sizes which can be used to self-assemble complexes with new functionality. Experimental techniques ranging from DNA-mediated self-assembly18–21 to entropic depletion interactions22,23 have been rapidly advancing during the last decade, thus providing a plethora of possibilities to fabricate colloidal particles whose shape and interactions can be controlled in detail. Moreover, external stimuli such as temperature, light, or pH can be used to control the inter-particle interactions during the assembly process.24,25 These techniques allow for a state- or time-dependent switching of the interactions and can be used to steer the assembly process. Controlling the particle interactions during the assembly process can prevent kinetic trapping resulting in a higher yield of the desired structure, as has recently been shown in a computational study for virus assembly.26 To fully exploit the potential of these techniques, a detailed understanding of the dynamics of the assembly process (distribution of intermediates, relevant time scales) is of crucial importance.

        Understanding the mechanisms governing chemical reactions has a long tradition in theoretical physics and chemistry. The most powerful analytical technique in this context is the Fokker-Planck or Smoluchowski equation, which has been used early to study bimolecular association based on diffusive motion. This approach was pioneered by Smoluchowski who first calculated the maximum diffusion-limited reaction rate for a fully reactive spherical particle.27 An important generalization of this result was derived by Collins and Kimball28 who studied the effect of finite reactivity by introducing a radiation boundary condition relating the concentration at contact to the reactive flux. Another essential extension is the case of particles with anisotropic reactivity ().29–34 However, most of these generalizations focus on the calculation of association rates and only a few aim at describing reversible reactions.35,36

        To study the full dynamics of protein assembly one needs to consider both association and dissociation processes. Moreover, even for globular proteins the intermediates formed during the assembly process are often of highly non-spherical geometry, thus rendering an analytical treatment of assembly in the framework of the Fokker-Planck equation very difficult. In this case computer simulations provide a valuable alternative. While detailed molecular dynamics (MD) simulations have been successfully used to investigate the behavior of a single molecule in great detail, studying the dynamics of large protein complexes on biologically relevant time and length scales is prohibited by the high computational costs of these simulations. Therefore, coarse-grained models are required for this case. Brownian dynamics (BD) simulations, the numerical counterpart to the Fokker-Planck equation, have been extensively used to study bimolecular association kinetics with realistic protein shapes.1,37–45 These studies have been focused mainly on a detailed calculation of association rates in bimolecular reactions based on realistic protein shapes and binding interactions. To study reaction dynamics in a large system consisting of many proteins, various simulation frameworks have been developed.46 They range from space- and time-continuous BD simulations as, for example, in the Smoldyn framework47,48 through event-driven tools like Greens functions reaction dynamics (GFRD),49–51 which is based on the analytical solution of the Fokker-Planck equation, up to a space-discretized version of Gillespie's chemical master equation approach,52 as for example used in MesoRD.53–55 While these simulation frameworks have been successfully used to study reaction kinetics on a large scale, they cannot be used to study the details of assembly processes as all of these frameworks lack a detailed description of protein anisotropy, including their shape or directional interactions. To study the assembly of virus capsids as a paradigm of a biological self-assembly process for which these elements are essential, coarse-grained MD simulations,56–59 and Monte Carlo (MC) techniques60–63 have been used. Similarly, the interaction of colloidal particles has been investigated on various scales with different simulation techniques,21 including MD simulations,64 MC simulations,65,66 and BD studies.67,68

        Here we introduce a simulation framework which allows us to study the spatial and stochastic aspects of protein assembly to large detail and nevertheless is computationally relatively cheap. In our approach, we combine BD of realistic protein shapes with stochastic reactivity for both, association and dissociation. Proteins are described as assemblies of spheres equipped with reactive patches. Their motion in real and orientation space is based on the overdamped Langevin equation with an anisotropic diffusion tensor. Inspired by the notion of an encounter complex,1,17,28,36,43,69,70 we decompose the reaction process into a diffusion and a reaction part. Upon the diffusive formation of an encounter, two particles can stochastically react with a microscopic reaction rate and thus form a bond. Similarly, an existing bond can be disrupted with a microscopic dissociation rate. We first show how these microscopic reaction rates can be related to macroscopic, experimentally measurable rates. We then verify that our algorithm correctly reproduces the macroscopically expected equilibrium behavior for bimolecular reactions. Finally, we investigate the assembly of a pentameric ring structure and compare our simulation results to a macroscopic rate equation approach. Here we again find excellent agreement between the stochastic simulations and the macroscopic description if the macroscopic rates are calculated without any free parameter in the correct way that includes both the diffusive and reactive contributions. These results show the importance of including spatial and orientational effects as well as realistic diffusion properties to understand the dynamics governing protein assembly.

         
        A.  Encounter complex

        Our simulation approach is based on the concept of the encounter complex.1,17,28,36,43,69,70 In this concept, a bimolecular reaction is decomposed into two steps: the undirected diffusive motion of the two binding partners and until they stochastically reach the encounter state · and the reaction from the encounter state to the bound complex . In terms of the free energy landscape the encounter complex represents a barrier separating the flat diffusive energy landscape from the reaction funnel. Thus, it is a transient state which can be characterized by the onset of highly correlated motion between the binding partners.43 A barrier in the free energy landscape might arise, for example, due to a necessary rearrangement of the binding site or due to a reorganization of the water layer. In some cases, there might not exist an explicit barrier, but even then the encounter complex is a helpful theoretical concept. Considering the encounter complex as the “watershed” between diffusion and reaction processes, the bimolecular reaction can be split according to the following reaction scheme:

        A+Bkk+C,
        (1) 
        A+BkD,bkDABkdkaC.
        (2) 
        Here and represent the overall kinetic rates of the reaction. In Eq. (2), the binding and unbinding process is split into a diffusive part represented by the rates and and a reaction part represented by and , respectively. Separating the binding process in free diffusive motion without steering forces and localized binding requires short-ranged interactions. Thus, our approach for the simulation of protein association corresponds to a regime of high salt concentration screening long-ranged electrostatic forces. Indeed this is a reasonable assumption for physiological salt conditions. Moreover, our approach can be used to study the assembly of μm-sized colloids functionalized with reactive patches where the size of the reactive region is small compared to the particle size. As we consider a situation in which the interaction between the proteins is short-ranged and only affects the reaction probabilities in the encounter state, individual clusters undergo free diffusive motion without an additional drift term. Moreover, here we do not consider hydrodynamic interactions which also might influence the binding kinetics.

        B.  Patchy particles

        Our simulation system can be decomposed into four elementary structures: proteins, patches, bonds, clusters. Clusters are rigid objects that can consist of a single or of multiple proteins held together by bonds. Two clusters can react with each other by bond formation between adhesive patches. A cluster consisting of multiple proteins can decay into smaller clusters by bond dissociation. Proteins are described by sets of non-overlapping, hard spheres approximating their shape. The proteins are equipped with reactive patches which reflect the localized binding sites. A patch is defined starting from a sphere of radius whose center c⃗ p (relative to the protein) can be chosen independent of the center of the protein. In addition, each patch is described by an orientation vector o⃗ . As a simple example for a non-globular protein, in Fig. 1(a) we schematically show a dumbbell-shaped protein consisting of two hard spheres in contact. The protein is equipped with one reactive patch of radius whose center coincides with the center of one of the steric spheres and is thus shifted by c⃗ p from the protein center located at the contact point of the steric spheres. The orientation vector o⃗  points along the long axis of the dumbbell. Using the above definitions, we formulate the following set of equations to define the encounter between two clusters mediated by a pair of patches, compare Fig. 1(b),

        rpRp,1+Rp,2,
        (3) 
        acos(r⃗ o⃗ 1|r⃗ o⃗ 1|)θp,1,acos(r⃗ o⃗ 2|r⃗ o⃗ 2|)θp,2.
        (4) 
        Equation (3) means that the distance between the centers of the two patches has to be sufficiently small for an encounter to occur and can be implemented easily in a simulation. Equation (4) is more complex. It involves the center-to-center vector r⃗  and not only reflects that the patches are anisotropic, but also assumes that an encounter only occurs if the two partners have a favorable orientation relatively to each other which is defined by the two parameters θ and θ. In general an anisotropic protein can be described by multiple spheres and the center of the patches can be chosen independent of the steric spheres, allowing for a variable description of proteins and their localized interactions.

        FIG. 1.

        Model definition. (a) A dumbbell-shaped protein modeled by two hard spheres. The protein is covered with a patch of radius whose center coincides with the center of one of the steric spheres and is thus shifted by c⃗ p relative to the center of the protein (contact point of the spheres). The patch has an opening angle of θ around the orientation vector o⃗  pointing along the long axis of the dumbbell. (b) The encounter between two of the proteins of (a). Here r⃗  is the center-to-center(ctc)-vector between the proteins and r⃗ p is the center-to-center vector between the patches. (c) Illustration of the reorientation processes upon reaction. In a first step, the ctc vectors of the proteins are aligned with the desired ctc-vectors. In a next step, the projection of the torsion vectors t⃗ 1 and t⃗ 2 on a plane perpendicular to r⃗ ctc are aligned. Finally, the clusters are shifted so that the desired relative distance between the proteins is achieved.

         

        C.  Particle motion

        In our approach, particle motion is described by the six-dimensional, overdamped Langevin equation (Brownian motion) describing translational and rotational diffusion of an arbitrarily shaped but rigid object:71,72

        tXt=gtwith:gt=0,gtgt=2kBTMδ(tt).
        (5) 
        Here is a six-dimensional position vector describing the position and orientations at time t and is Gaussian white noise. is the mobility matrix which is calculated on-the-fly for any cluster shape71–73 following the method of de la Torre.74 Because intermediates continuously grow and shrink due to association and dissociation, the diffusive properties of the involved clusters are continuously changing during a simulation.

        Protein particles are modeled as hard spheres and their collision requires special attention. In a BD framework, particle velocities are not defined and thus a ballistic reflection scheme is not appropriate. Here we treat steric collisions by a rejection method similar to MC simulations. If a steric overlap between two clusters is created during a move step, this move step is rejected. One method to deal with this situation would be to repeatedly choose a new position of the two clusters, starting from their original location and orientation, until the move step is accepted. This method leads to incorrect simulation results as is demonstrated for the example of two hard spheres in a periodic box in Fig. 2. Here, a histogram of the distance distribution between the two spheres, scaled by the volume of a spherical shell, is shown in a range in which boundary effects are not observed. By choosing a new position until successful acceptance of the move step (green histogram), an effective repulsion between the two spheres is introduced, leading to a depletion zone at small particle distances instead of the expected uniform distribution. If, instead of choosing a new position, we set the two clusters involved in a steric collision back to their original position (and orientation) before the move step, we recover the expected uniform equilibrium distribution (shown in red in Fig. 2). Thus, this method ensures that the correct equilibrium distribution is realized.

        FIG. 2.

        Rejection methods. The distance distribution of two hard spheres of radius = 1 nm in a periodic box of volume Vbox=(8nm)3 is shown for two different rejection methods. In the first rejection method (green), new positions of the particles are chosen (starting from their positions before collision) until no overlap is observed. In the second rejection method (red) the spheres are set back to their position before the collision. The simulations were performed at a time-step resolution of Δ = 0.1 ns.

         

        D.  Local rules

        In our approach, clusters are rigid assemblies consisting of one or multiple proteins. Assuming that the assembly geometry is determined by unique local interactions, a set of local rules is necessary to describe the new relative orientation and position of two reacting clusters.75 These local rules are encoded in the bond that is formed between the clusters. Each bond carries the information of the desired center-to-center (ctc) vector of the two proteins directly involved in the binding process as well as their relative orientation (torsion vectors). Upon reaction the two clusters instantaneously flip into the desired relative configuration. This rule follows from the assumption that the short-ranged forces involved in the final binding step lead to a fast rearrangement on the time scale of our BD simulation. The reorientation process is schematically shown in Fig. 1(c). In a first step, the clusters are shifted and rotated so that the predefined ctc-vector matches the real ctc-vector. This procedure enforces the correct position of the two merging clusters relative to each other. In a next step, the relative orientation of the clusters is corrected by aligning the projection of the two torsion vectors on a plane perpendicular to the ctc-vector. The necessary rotation and translation of the clusters is distributed between them according to their diffusive weights. This means that for a small and a large partner, essentially only the small partner is moved, as one expects for physical reasons; for two similarly sized partners, both are moved to a similar extent. Since all interactions are local (e.g., local patches or constraints on orientation/torsion), the encounter complex corresponds to a region in configuration space surrounding the desired relative position and orientation. If this region is sufficiently small, the instantaneous flipping into the correct configuration is comparable to the resolution of the BD simulation. Reorientation of the clusters during the binding process can lead to a steric overlap either with another cluster or between the merging clusters. If this is the case, the reaction is rejected and the old positions and orientations of the clusters before the move step are resumed. The steric collision of two merging clusters reflects that during the assembly process incompatible fragments can prevent the formation of the desired structure. For example, if the desired structure encoded by the local bonds is a ring consisting of five proteins, two ring fragments each containing three proteins cannot bind with each other.

        E.  Microscopic and macroscopic rates

        In order to implement full reversibility, both reaction directions are treated as stochastic reactions with two corresponding microscopic rates. If the reactive patches of two clusters form an encounter, specified by the constraints in Eqs. (3) and (4), they can react with a bond-specific rate ka. Similarly, an existing bond can dissociate with a bond-specific rate kd. Here we assume that the waiting times for association or dissociation of a bond between two clusters in encounter are Poisson-distributed, P(t,ki)=kiexp(kit), = {, }. In this case, the probability that no association or dissociation has occurred after a timestep Δ is given by S(Δt,ki)=1Δt0P(t,ki)dt and hence the probabilities for bond dissociation or bond formation are given by50

        PassocPdissoc==1S(Δt,ka)kaΔt,1S(Δt,kd)kdΔt.
        (6) 
        The approximations used in Eq. (6) are valid if kdΔt1 and kaΔt1, which is the case throughout this work.

        We now show how the microscopic reaction rates ka and kd can be related to macroscopic reaction rates. In a macroscopic framework, the reaction scheme in Eq. (2) can be interpreted as a system of ordinary differential equations describing the changes in the concentrations ,,, and of the different species involved in the reaction.69 In this framework, the encounter is understood as a single, intermediate state connecting the bound complex and the unbound and particles. The encounter can either react to the final complex with the first order rate or decay into two separated particles with the first order rate . An encounter can be formed by the first order decay of cluster with a rate or by the diffusive association of an and particle described by the second order rate . Here has the dimension m3/s, if we consider particle concentrations in 3 dimensions. Using a steady state approximation of the encounter complex, the following relations are obtained for the overall forward and backward rates and the equilibrium constant:69

        k+=kDkaka+kD,b,k=kD,bkdka+kD,b
        (7) 
        Keq=k+k=kDkD,bkakd.
        (8) 
        In the case of a diffusion-limited reaction .69 In this case, the overall forward and reverse rates simplify into the purely diffusion-limited forward rate and the backward rate becomes /. In the case that the reaction is reaction-limited (), the overall forward rate becomes / and the overall backward rate is .

        In the framework of the Fokker-Planck equation, two spherical particles are considered to be in an encounter if they are in contact.27,28 Finite reactivity for particles in encounter was introduced by Collins and Kimball28 in the form of a radiation boundary condition in which the rate κ (termed κ in Ref. 70 and in Ref. 28) relates the concentration at contact to the reactive flux. By comparing the escape probabilities from the encounter in the Fokker-Planck framework and the macroscopic rate equation framework, Shoup and Szabo showed that κ can be related to macroscopic rates by identifying κ = / and Keq=κa/kd.70 Agmon and Szabo derived the same relation for the equilibrium constant (Keq=κa/kd) by considering an isolated pair of reactive spheres using a back-reaction boundary condition.36

        Here we show how the microscopic reaction rate used to describe reactions in our simulation approach can be related to the reaction scheme in Eq. (2) and to the equilibrium constant (Eq. (8)). Inspired by the work of Shoup and Szabo,70 we calculate the equilibrium constant for the formation of the encounter. In our approach the encounter is defined as a region in configuration space (Eqs. (3) and (4)) around the desired relative position in the bound complex instead of a boundary as commonly used in the Fokker-Planck picture. The configuration space of two rigid unbound particles in a periodic box of volume can be described by their relative position vector r⃗ , their center-of-mass vector R⃗  and 2 × 3 angular coordinates ω⃗ 1 and ω⃗ 2. Assuming free diffusive motion without reactions all configurations in the two particle configuration space are equally probable and they only interact by steric repulsion. By using the thermodynamic extremum principle for the Gibbs free energy, we can relate the equilibrium constant for the formation of an encounter to the partition sums of the free molecules and and to the partition sum of the encounter complex :44,76

        Kenceq=kDkD,b=(zAB/V)(zA/V)(zB/V),
        (9) 
        zA=zB=Vd3xd3ω=V×4π×2π,
        (10) 
        zAB=Vd3Rd3ri=12(d3ωi),
        (11) 
        Kenceq=d3ri=12(18π2d3ωi)=:V.
        (12) 
        For the single molecules all positions and orientations (neglecting the excluded volume) are accessible, resulting in a factor of from the integration over the translational degrees of freedom and a factor of 8π2 from the orientational degrees of freedom of the rigid molecules. To determine the partition sum of the encounter complex , we first integrate out the center-of-mass coordinate R⃗  of the two molecules resulting in a factor of . The boundaries of the remaining integrals over the relative position coordinate r⃗  and the orientation coordinates ω⃗ 1 and ω⃗ 2 are determined by the patch definition and the particle geometries. They follow from Eqs. (3) and (4) and have to be calculated for every specific case, as will be discussed in more detail below. The resulting reactive volume is the central concept to relate microscopic and macroscopic rates. From Eq. (12), we see that it equals the equilibrium constant Kenceq.

        If we now allow for reactions, not every configuration in the encounter region is equally populated as some particles will already react before they explore the inner region of . However, if on average the encounter region is well explored before an association, we can approximate the real distribution within by a uniform distribution. This approximation is equivalent to the underlying assumption in the work of Pogson 77 and Klann 78 who introduced a reactive volume (and an intrinsic association rate78) based on the assumption that the particles are randomly distributed in the box at all times.

        Alsallaq and Zhou used thermodynamic arguments to determine the equilibrium constant for the bound complex.44 However, in contrast to this work we consider all configurations within the generalized reactive volume as encounter configurations. By introducing finite reactivity with the microscopic reaction rate ka with which a bound complex can be formed from the encounter region (see Fig. 3(a)), our approach can account for reactions which are reaction- or diffusion-limited, as will be shown below in Sec. III. To study assembly, this is of great importance as intermediates emerging during the assembly process can have very different diffusion properties and a reaction that was reaction-limited for small clusters can become diffusion-limited for larger clusters.

        FIG. 3.

        Illustration of the reactive volume . The configuration space is defined by the relative position r⃗  of the particles and the 2 × 3 dimensional orientation vector Ω⃗ . (a) In the case of one reactive patch corresponds to the region in configuration space in which two particles are considered to be in encounter (Eqs. (3) and (4)). (b) For particles with multiple reactive patches different regions Vi,j in configuration space correspond to an encounter mediated by one particular patch combination ,. For identical microscopic rates and non-overlapping Vi,j the total reactive volume between the two particles is given by the sum over the encounter volumes associated with each patch combination.

         

        Thus, by identifying the equilibrium constant for the formation of the encounter with the reactive volume (Eq. (12)) and relegating bond association and dissociation to the reaction rates, we can identify the reaction rates ka and kd in the microscopic framework with the corresponding reaction rates and in the macroscopic framework (Eq. (8)) and the equilibrium constant is given by

        Keq=kDkD,bkakd=Vka/kd.
        (13) 
        In accordance with the work by Shoup and Szabo,70 the rate κ = is given by the product of the equilibrium constant for the formation of the encounter complex and the reaction rate for the formation of the final complex. This rate can be identified with the reaction rate used by Morelli and ten Wolde50 ( in Ref. 50).

        While the reactive volume = / is sufficient to predict the correct equilibrium constant (Eq. (13)), the kinetics of the reactions depend on the absolute values of the diffusive rates and . Thus in order to determine the macroscopic rates and defined in Eq. (7), we need to evaluate . This can be done within our simulation approach based on an algorithm by Zhou.39 In this algorithm, can be calculated from the survival probability () of two particles starting in encounter by39,43

        kD=limtκaS(t)1S(t).
        (14) 
        Given the diffusive on-rate , the diffusive off-rate can be simply calculated by = / (Eq. (9)).

        Any simulation algorithm with underlying reversible dynamics needs to satisfy detailed balance. To establish detailed balance we follow the approach by Morelli and ten Wolde50 who inferred from a detailed balance consideration that the relative position distribution of two spherical particles prior to a reaction has to be identical (after renormalization) to the position distribution of the two particles after dissociation. For two hard spheres covered with a reactive shell, it has been argued that detailed balance can be introduced by placing the two hard spheres according to a radial uniform distribution within the reactive shell followed by one additional move step.79 Here we generalize this idea as follows: Upon the dissociation of two clusters, a new configuration is chosen uniformly within the encounter region followed by a diffusive move step of the two clusters. To establish the uniform distribution in we exploit the time invariance symmetry of the diffusion propagator. This symmetry ensures that in a confined volume of the configuration space, a uniform distribution is established as the steady state distribution. As describes a configuration around the desired relative configuration of the proteins, we can generate a uniform distribution without prior knowledge of the exact shape of the reactive volume by performing various “pseudo-diffusion” steps of the two clusters starting from their predefined relative configuration. These “pseudo-diffusion” steps are confined to the region in the relative accessible configuration space of the two partners. We call this procedure “pseudo-diffusion” as the two dissociating clusters are propagated diffusively within , however, the timestep used to establish the uniform distribution has no physical meaning. This procedure for generating a uniform distribution within the reactive area does not depend on the particle or patch geometry. Thus it is especially suited to study assembly dynamics as parts of the reactive volume might become sterically blocked during the assembly process. These steric effects are automatically taken into account with our method. The effect of steric blocking and its consequences for assembly will be discussed in Sec. III. Moreover, this method would also remain valid when including long range electrostatic interactions. In this case, the distribution within the encounter complex in equilibrium would not be uniform. However, the “pseudo-diffusive” motion step would lead to the expected distribution within .

        In general, the encounter volume has to be understood as the volume of all two particle configurations resulting in an encounter between two clusters (see Fig. 3(b)). This means that for clusters with multiple patches the total encounter volume V=Vpi,pj where Vpi,pj is the encounter volume between a pair of patches. In the case of a dissociation of clusters with multiple patches a uniform configuration in the subvolume Vpi,pj of the patches that formed the bond is realized. If the different patches have a different microscopic reactivity, the subvolumes have to be treated separately and the equilibrium constant for the complex formation is given by the sum of the different equilibrium constants for different bonds.

        F.  Intrinsic association

        Proteins with multiple reactive patches can assemble in structures containing loops. In this case an already connected cluster can contain open bonds with two patches being in permanent encounter due to the geometry of the cluster. These patches can bind with the internal association rate kintraa, where we again assume Poisson-distributed waiting times for the bond formation. The bond formation in this case is fundamentally different from the above discussed process as no diffusion is involved. Since clusters are rigid objects, internal bond formation does not change the shape of the cluster. Thus it can be considered as an internal process stabilizing an already existing structure. Given a dissociation rate of for a specific bond, we can relate the internal association rate to the free energy of bond formation. For a closed loop containing bonds, there are different configurations with one open bond in the loop and the fraction of the open and closed loop in equilibrium is given by

        popenpclosed=nkdkintraa=ZopenZclosed=neβEkintraa=kdeβE.
        (15) 
        Thus, the internal association regulates the fraction of open and closed loop configurations. In a more realistic scenario, the dissociation of one bond in a loop will enable an enhanced internal movement. This will be reflected in an increase in the entropy which would shift the ratio of open and closed ring structures towards the open state (smaller kintraa). In principle it should be possible to calculate the change in energy and entropy upon dissociation of a bond using detailed MD simulations. Such simulations could be used to specify the internal association rate.

        G.  Outline of the algorithm

        In this part, we will describe the work flow of our algorithm. Once a starting configuration is seeded the following steps are repeated iteratively. First, every particle is translated and rotated according to its diffusive properties in the framework of free Brownian motion within a timestep Δ. With the new positions and orientations all clusters with steric overlaps as well as all clusters participating in a reaction are tagged. A reaction between two clusters occurs with the bond specific probability Δ if they are in encounter. In a next step, all steric overlaps (including box overlaps if the box is non-periodic) are corrected by setting the involved particles back to their original position before the move step. This might lead to further steric overlaps as a higher order effect. These overlaps are then also corrected. In this step, clusters tagged for a reaction click into their predefined relative positions and orientations (see Fig. 1(c)). If a reaction leads to a steric overlap either between the two merging clusters or with other clusters, the reaction is rejected and the particles reassume their configuration before the move step. Collisions with other clusters after the reorientations step is a higher order effect in the concentration of the particles. In dilute systems, especially for reactive configuration volumes around the desired relative configuration, this happens rarely. In a next step open, internal bonds can be closed with probability Paccintra=kintraaΔt. Finally, each existing bond can dissociate with the bond specific probability Δ. If this leads to two unconnected cluster fragments, the diffusive properties of the fragments are calculated. Subsequently, a uniform distribution within the encounter volume is realized by “pseudo-diffusive” motion preserving the encounter followed by one unconstrained diffusion step of the two clusters. All simulations have been performed at room temperature = 293 K and with the viscosity of water η = mPa s.

         
        A.  The reactive volume

        As discussed in Sec. II E, we can relate the microscopic rates to the macroscopic equilibrium constant using the concept of a generalized reactive volume in configuration space (Eq. (13)). Here we introduce a model system in which it is possible to analytically calculate and show that the results from our MC scheme (Eq. (17)) agree well with the analytical results. In this model system, a protein is described by one hard sphere of radius R equipped with one polar patch as depicted in Fig. 3(a). The patch is centered at the origin of the protein. It has a radius of and an opening angle of θ ∈ [0, π] around the z-axis of the proteins. With our definition of the encounter complex for two such particles (Eqs. (3) and (4)), we can analytically calculate by decomposing the radial and orientational part, similar to the potential model proposed by Kern and Frenkel,80

        VVradVori===Vrad×Vori43π((Rp)3R3),(14π2π0dϕθp0sin(θ)dθ)2=14(1cos(θp))2.
        (16) 
        The square in the orientation part arises from the fact that we have to consider two independently rotating particles and only if the orientation of both particles is correct, an encounter is reached (see Eq. (4)). For θ = π the orientation constraint vanishes and the spherical result is recovered.

        As analytically calculating the reactive volume is only feasible for some special geometries, it can be numerically pre-calculated for the elementary assembly blocks using a MC integration scheme. In this scheme, we sample different configurations by randomly positioning two clusters with random orientations in a periodic box of volume Vbox. We repeat this procedure times and can estimate by counting the fraction of trials n=Nencounter/Ntotal that result in an encounter configuration. Then we simply have

        VMC=nVbox.
        (17) 
        In Fig. 4, the configuration space volume for two identical, spherical proteins equipped with a polar patch ( = 1 nm and = 1.1 nm) is shown as a function of the opening angle θ. Comparing the analytical result (black line) and MC simulation (Eq. (17)) we see that our simple MC scheme correctly estimates . For proteins equipped with multiple, non-overlapping patches which can bind with each other via different patch combinations, we have to distinguish between the encounter volume Vpi,pj of two specific patches and and the total encounter volume of two proteins. The total configuration volume is given by the union of all volumes Vpi,pj for the encounter of two different patches , that can bind to each other. If the encounter volumes of different patch combinations are disjoint, the total configuration space volume for the encounter of two clusters is given by summing over all Vpi,pj. For multiple, identical patches that can bind with each other, the total reactive volume is given by multiplying the configuration volume Vpi,pj of one patch combination with a combinatorial prefactor accounting for the multiplicity of the different patch configurations if the corresponding volumes are not overlapping in configuration space.

        FIG. 4.

        for two identical hard spheres of radius = 1 nm equipped with one polar patch of radius = 1.1 nm as a function of the patch angle θ. The analytical result from Eq. (16) is shown with a dashed line, while the MC estimates (according to Eq. (17)) are represented by red points with error bars.

         

        B.  Bimolecular reactions

        To verify our algorithm and to show the importance of detailed balance, we first analyze a system with a small number of proteins only. In detail we consider a situation in which one particle of type and particles of type are placed in a periodic box of volume Vbox. In this case, we can relate the concentration of the complex to the fraction of time in which particle is bound to particle ( = /( + )), while the concentration of particle can be related to the fraction of time in which is unbound ( = /( + )). The unbound particle is surrounded by a concentration of = / particles of type . Here V=VboxVexcl is the freely accessible volume. For proteins of identical radius ( = = ) the excluded volume can be approximated by Vexc=4πNB(2R)3/3. Thus we can relate the probability that is bound to the equilibrium constant of the reaction by36,50

        Keq=cCcBcA=pboundNBV(1pbound)pbound(NB)=KeqNBKeqNB+V.
        (18) 
        Equation (18) follows from elementary thermodynamic arguments and is valid irrespective of the details of the binding interactions. In Subsections III B 1 and III B 2, we will show that the correct equilibrium bound probability is reproduced with our algorithm for spherically reactive and patchy particles when appropriately taking the generalized encounter volume into account. We also show that detailed balance is essential to reproduce from Eq. (18).

        1.  Spherically reactive particles

        Here we compare our BD algorithm with the mean field result (Eq. (18)) for the case of spherically reactive and particles ( = = 1 nm, = = 1.1 nm, θ = θ = π, = 11.09 nm3) and show that our algorithm fulfills detailed balance in this case. Given an association rate of = 10 ns−1 we chose so that Keq=Vkakd=V. In Fig. 5(a), the probability of particle being bound to particle is shown as a function of the number of particles and compared to the mean field result (black line). Placing particles uniformly in the reactive shell by “pseudo-diffusive” motion and allowing for one unconstrained move step of the two proteins involved in the dissociation according to the algorithm described above (circles) leads to very good agreement between the mean field theory and the simulation results. Only for a very large timestep Δ = 0.01 ns leading to a large reaction probability of 0.1, a small deviation between the simulation results and the expected mean field theory is observed. On the contrary, violating detailed balance by placing particles in contact after the dissociation (stars) leads to an overestimation of pbound and the simulation results cannot be reconciled with the mean field theory. Although reducing the timestep from Δ = 0.01 ns to Δ = 0.001ns leads to a modest improvement for the simulation results with contact dissociation, no further improvement is observed when decreasing the timestep further. This shows that detailed balance is essential in the dissociation step for the agreement between simulation results and macroscopic theory. These findings agree with those obtained by Morelli and ten Wolde who performed a similar simulation to test their algorithm.50

        FIG. 5.

        Equilibrium bound probabilities and radial distribution for spherically reactive particles. (a) The probability of one particle being bound to a particle is shown as a function of the number of particles according to Eq. (18) (black line) for = . The following parameters have been used for the simulations: = = 1 nm, = = 1.1 nm, θ = θ = π leading to ≈ = 11.09 nm3 (Eq. (16)). For a given = 10ns−1, the dissociation rate is chosen so that Vkakd=Vbox=8000 nm3 (Eq. (13)). Circles depict the simulation results obtained with our algorithm for different time resolutions, while stars correspond to simulations in which detailed balance is violated by placing particles in contact after dissociation. (b) The normalized radial distribution Δ before reaction (circles) and after dissociation including the additional unconstrained move step (stars) is shown for two different association rates = 10 ns−1 (red) and = 0.1 ns−1(green) with a time resolution of Δ = 0.01 ns. The inset shows the comparison of the uniform radial distribution established by “pseudo-diffusive” motion constraint to (red points) and the expected uniform distribution in a spherical shell (black line).

         

        In Fig. 5(b), we verify that our algorithm indeed obeys detailed balance by comparing the radial distribution of particles before a reaction (circles) and after dissociation (stars) for two different association rates. Irrespective of the choice of we observe the same radial distribution before reaction and after dissociation and thus our algorithm satisfies detailed balance for this quasi-one-dimensional case. In the inset of Fig. 5(b), we show a histogram of the distance between two dissociating particles after the “pseudo-diffusive” motion constraint to without the additional, unconstrained move step of the two particles. Here we see that we are indeed able to generate a uniform distribution within by exploiting the symmetry of the diffusion propagator as the histogram of the simulated positions agrees very well with the expected uniform distribution.

        2.  Patchy particles

        After verifying that our algorithm reproduces the expected results for spherically reactive particles, we will now investigate the effect of anisotropic reactivity and show that our generalized definition of the reactive volume enables us to compare our simulation results with the mean field theory, which is the same as for spherically reactive particles. Here we consider the same setup as in the previous part with one particle of type and particles of type in a periodic box ( = = 1 nm and = = 1.1 nm). This time, however, both particles are equipped with a polar patch (θ = π/4) allowing only for an association of two particles if the patches are sufficiently aligned with the ctc-vector of the proteins (Eq. (4)). In this case, the reactive volume reduces to = 0.23 nm3 according to Eq. (16). Given the reactive volume and the association rate = 10 ns−1, we again choose in such a way that Keq=V. The comparison between the simulation and mean field theory is shown in Fig. 6. By using the appropriate scaling with the generalized volume for the patchy particles, we observe perfect agreement of the mean field results from Eq. (18) (black line) with the simulation results when placing the particles according to the proposed detailed balance algorithm (circles). When placing the particles in contact (alignment of patch orientation with ctc-vector and radial contact), we drastically overestimate the bound probability (stars). Compared to the spherical case (Fig. 5), this overestimation of pbound is much stronger. Thus, in the case of localized reactivity it is even more important to carefully consider the relative position and orientation of dissociating proteins. Moreover, our results confirm the relation between microscopic rates and the macroscopic equilibrium constant given in Eq. (13).

        FIG. 6.

        Equilibrium bound probabilities for patchy particles. The probability of one particle being bound to a particle is shown as a function of the number of particles according to Eq. (18) (black line) for Keq=Vbox. The mean field theory is independent of the specific particle details. The following parameters have been used for the simulations: = = 1 nm, = = 1.1 nm, θ = θ = π/4 leading to ≈ = 0.23 nm3 (Eq. (16)). For a given = 10 ns−1 the dissociation rate is chosen so that Vkakd=Vbox=8000 nm3. Circles show the results of our algorithm for various timesteps, while stars correspond to simulations in which detailed balance is violated by placing particles in contact and with perfect alignment of the patches after dissociation.

         

        FIG. 5.

        Equilibrium bound probabilities and radial distribution for spherically reactive particles. (a) The probability of one particle being bound to a particle is shown as a function of the number of particles according to Eq. (18) (black line) for = . The following parameters have been used for the simulations: = = 1 nm, = = 1.1 nm, θ = θ = π leading to ≈ = 11.09 nm3 (Eq. (16)). For a given = 10ns−1, the dissociation rate is chosen so that Vkakd=Vbox=8000 nm3 (Eq. (13)). Circles depict the simulation results obtained with our algorithm for different time resolutions, while stars correspond to simulations in which detailed balance is violated by placing particles in contact after dissociation. (b) The normalized radial distribution Δ before reaction (circles) and after dissociation including the additional unconstrained move step (stars) is shown for two different association rates = 10 ns−1 (red) and = 0.1 ns−1(green) with a time resolution of Δ = 0.01 ns. The inset shows the comparison of the uniform radial distribution established by “pseudo-diffusive” motion constraint to (red points) and the expected uniform distribution in a spherical shell (black line).

         

        In Fig. 7, we again show the distribution of the relative position and orientation of two particles equipped with a polar patch before a reaction and after dissociation. Here ϕ corresponds to the angle between the orientation vectors of the patches and is the distance between the centers of the two particles. Due to the symmetry of the particles around the z-axis these two parameters provide a good description of the relevant configuration space. The normalized probability distribution (, ϕ)/(4π2 × 2πsin(ϕ)) is shown for two microscopic association rates = 0.1 ns−1 (Fig. 7 upper part) and = 10 ns−1 (Fig. 7 lower part) at a time resolution of Δ = 0.01 ns−1. In the left part of Fig. 7, the distribution before reaction is shown, while in the right part the corresponding distribution after dissociation is depicted. As expected the distribution has its maximum for particles in contact and anti-parallel aligned patch orientations ( = 2.0 nm and ϕ = π). Although the distributions fluctuate more for a smaller association constant (Figs. 7(a) and 7(b)), as can be seen by looking at the contour lines, very good agreement between the distribution before association and after dissociation is observed for both association parameters. This shows that our algorithm maintains detailed balance in this case of non-spherical reactivity. Changing the association rate by two orders in magnitude does not affect the distributions suggesting that our algorithm works well for a large spectrum of microscopic parameters. This first example of non-spherically reactive particles shows that it is essential to appropriately scale the macroscopic parameters using the concept of a generalized reactive volume in configuration space. Moreover, careful consideration of the relative position and orientation of the dissociating particle is necessary to obtain quantitative agreement between BD simulations and mean field theory or experimental results.

        FIG. 7.

        Histograms of the particle positions before reaction and after dissociation as a function of the relative distance and the angle between the orientations of the two patches ϕ. This histogram has been weighted by the generalized bin volume (4π2 × 2πsin (ϕ)ϕ). Several contour lines are shown for better comparability. In (a) and (c), the distribution of particle positions before a reaction is shown for = 0.1 ns−1 and = 10 ns−1, respectively. In (b) and (d), the corresponding distributions of particle positions after dissociation is shown. The distributions have been recorded with a timestep resolution of Δ = 0.01 ns.

         

        C.  Beyond bimolecular reaction: Assembly of a pentameric ring

        In this section, we analyze the assembly of a pentameric ring structure consisting of 5 proteins (Fig. 8(b)) as an example of an assembly process involving more than two partners. As elementary building blocks we consider spherical proteins of radius = 1 nm, each equipped with two reactive patches (Fig. 8(a)). The patches are centered at the origin of the protein and both have a radius of = 1.1 nm and an opening angle of θ = π/5. The geometry of the ring is reflected in the different orientation vectors o⃗  associated with each patch, and the relative position and orientation of the proteins are encoded in the bond structure.

        FIG. 8.

        Simulation of the assembly of a pentameric ring structure. (a) Elementary spherical building block of radius = 1 nm for the pentameric ring equipped with two patches. The angle between the centers of the patches is determined by the desired ring structure. Each patch has an opening angle of θ = π/5 and a radius of = 1.1 nm. (b) Fully assembled ring structure. (c) Simulation snapshot of a typical box containing 500 proteins.

         

        Here we develop a rate equation approach based on a set of ordinary differential equations (ODEs) to describe the changes in concentration of the different fragments sizes. We demonstrate how the macroscopic reaction rates can be calculated from the microscopic reaction parameters and diffusive properties of the intermediates. The parameter free comparison of our simulation results with the rate equation predictions shows excellent agreement. This demonstrates that we are indeed able to predict macroscopic reaction rates from our simulations taking into account changes in diffusive properties of the assembly intermediates as well as steric effects arising during the assembly process.

        1.  Macroscopic rates

        In Eq. (7), the reaction rates and are given in the case of a bimolecular reaction. We now specify the macroscopic rates ki,j+ and ki,j for the reaction between two ring fragments with and proteins. These rates are given by

        ki,j+=ki,jDkaki,jD,b+kak˜i,j+(112δi,j),
        (19) 
        ki,j=ki,jD,bkdki,jD,b+kak˜i,j(2δi,j).
        (20) 
        Here ki,j+ describes the association of two ring fragments of size and . As the diffusive properties and the accessible reaction volume change with the size of the fragments ki,j+ is different for different combination of fragments, albeit the microscopic reaction rate remains constant. For the reaction of identical particles the macroscopic rate is reduced by a prefactor of 1/2 compared to the simulation rate . This reflects the fact that for reactions of identical components ( + ) the observed macroscopic rate is by a factor of 1/2 smaller than the rate used in the simulation due to the number of collisions being proportional to ( − 1)/2!.47,52 ki,j describes the macroscopic dissociation rate of a cluster of size + into two ring fragments of size and , respectively. To determine the macroscopic reaction rates given in Eqs. (19) and (20) we need to calculate the diffusive on- and off-rates ki,jD and ki,jD,b=ki,jD/Vi,j which change for different combination of fragments. Here Vi,j is the configuration volume associated with the reaction of two fragments of size and .

        The reactive volume for two different fragment configurations and can be calculated by MC sampling as has been described in Eq. (17). The resulting volumes Vi,j for two fragments of size and are shown in Table I. For the reaction of two monomers the reactive volume V1,1 can also be evaluated by using Eq. (16) for the encounter of one specific pair of patches. Taking into account the multiplicity of bond interactions 22 results in a total encounter volume of V11=0.4 nm3. While fragments not combining into a full ring ( + < 5) have the same configuration volume as the monomeric fragments, the configuration volume for fragments that can combine into a full ring is reduced. This reflects the fact that, independent of the particle geometry, the encounter volume and thus the number of configurations in this volume remains unchanged, unless parts of the reactive volume become sterically blocked, which is the case for fragments that can assemble into the full ring structure. Moreover, for these fragments some configurations with a simultaneous overlap of different patch combinations exist. As the encounter volumes for different patch combinations are not disjoint in this case, the total encounter volume is smaller than the sum of the subvolumes. These effects are naturally included in our simulation algorithm and for the simulations only knowledge of the reduced volume for the basic building blocks V1,1 is required. However, in order to compare our simulation results to a macroscopic rate equation approach, we need to determine the changes in the encounter volumes for different assembly intermediates.

        Table I.

        Generalized reactive volumes Vi,j in nm3 for the encounter of different ring fragments calculated by MC sampling according to Eq. (17) and analytically according to Eq. (16) for V1,1, respectively.

         

        Moreover, in order to capture the reaction kinetics with the above defined rates, we need to determine the diffusive on-rates ki,jD. In order to calculate the diffusive on-rate ki,jD, we follow a scheme which was originally proposed by Zhou.39 Here, two ring fragments of size and are initially placed in a random configuration in encounter (e.g., a random configuration within Vi,j is chosen). Starting from this configuration the two fragments are propagated diffusively within our simulation framework. If the two fragments form an encounter, they can react with the probability Δ. In this case the run is terminated. By recording the survival probability , () (the fraction of runs which survive until time ) ki,jD can be estimated according to Eq. (14) with κi,ja=kaVi,j. In order to calculate , ( → ∞) we use a microscopic reaction rate of = 1 ns−1 and an adaptive timestep scheme with a minimum timestep of Δtmin=0.01ns which ensures that all steric collision and encounters are detected on the resolution of the minimum timestep.39 The diffusive off-rate for the decay of a fragment into two smaller fragments of size and is given by ki,jD,b=ki,jD/Vi,j. The absolute values for the different ring fragments are given in Table II. As fragments with ( + ) > 5 cannot combine with each other due to steric collisions, the diffusive on-rate in this case is ki,jD=0.

        Table II.

        Diffusive on-rates for different ring fragments in nm3/ns. The rates have been calculated from the survival probabilities of fragments starting in a encounter according to the algorithm proposed by Zhou39 with a minimum timestep resolution of Δ = 0.01ns.

         

        2.  Macroscopic rate equation approach

        After having defined the macroscopic reaction rates for the association and dissociation of different ring fragments (Eqs. (19) and (20)) we will now introduce a system of ODEs to describe the changes in concentration of a ring fragment of size . This description assumes that the system is homogeneous and well mixed at all times. Stochastic fluctuations arising from the finite number of proteins are neglected. While for ring fragments containing less than 5 proteins only one state exists, two states exist for the full ring: an open ring with 4 bonds and a closed ring with 5 bonds. The closed ring can be considered as an internal state as the last bond formation in our framework does not involve any diffusion and it is connected to the open state by the internal association rate kintraa (see Eq. (15)). The concentration of the open state is denoted by , while the concentration of the closed state is denoted by c5. The changes in concentration for a ring consisting of proteins can be described by the following set of ODEs:

        ċ i=l=1Niki,l+(1+δi,l)cicll=1lili1kl,ilcikintraaδi,NcN+l=1lili1kl,il+clcil+l=1Niki,l(1+δi,l)cl+i+Nδi,NkdcN
        (21) 
        =l=1Nik˜i,l+cicll=1i1k˜l,ilcikintraaδi,NcN+12l=1i1k˜l,il+clcil+2l=1Nik˜i,lcl+i+Nδi,NkdcN
        (22) 
        ċ N=NkdcN+kintraacN.
        (23) 
        In our case = 5. Equation (21) shows the different processes which lead to a change in the concentration of a fragment containing proteins. The first three terms lead to a decrease in concentration. The first term describes the change in concentration due to the reaction of a fragment of size with another fragment of size . If = , a twofold decrease in concentration is observed. The second term describes the change in concentration due to the decay of a fragment of size into two fragments of size and . To prevent the double counting of dissociation events the constraint for the summation is used here. In the third term, the internal bond formation from the open to the closed state is described. Similarly in the last three terms all contributions which lead to an increase in concentration due to association and dissociation processes are described. In Eq. (23), the dynamics of the closed pentameric ring state is separately described. In the derivation of Eqs. (22) from (21) the symmetry of the matrices ki,l±=kl,i± has been exploited. The above described set of ODEs obeys mass conservation (i(ci+δi,5c5)i=const.). Thus, one of the concentrations could be eliminated as it is dependent on the other concentrations.

        3.  Comparison between rate equations and BD simulations

        After having defined the macroscopic rate equation approach for the changes in fragment concentrations in Eq. (22) we will now compare our BD simulation results to the macroscopic rate equation approach. As all macroscopic rates are fully specified by the microscopic parameters (Eqs. (19) and (20)), the comparison between the simulation and macroscopic theory does not involve any free parameter. For the BD simulations we randomly position = 500 single proteins in a periodic box of size Vbox and record the time-course of the fragment concentrations. All BD simulations were performed at a timestep resolution of Δ = 0.01 ns and the resulting concentrations were averaged over 40 independent trajectories. For the comparison with the rate equation approach we use an initial concentration of = ( = 0) = /. Here is the accessible simulation volume which is roughly estimated from the box volume by V=Vbox(N1)4π(2R)3/3.

        In Fig. 9, the evolution of the normalized fragment concentrations c˜i=(ci+δi,5c5)i/c0 is shown. The normalized concentration of a fragment is defined as the concentration weighted by the number of proteins contained in the fragment and normalized by the initial concentration. The normalized concentration predicted by the rate equation approach is shown with solid lines, while the averaged simulation results are represented by dashed lines. In the left column of Fig. 9 (Figs. 9(a), 9(c), and 9(e)), a microscopic reaction rate of = 1 ns−1 and a box volume of Vbox=(55nm)3 are chosen. As the diffusive off-rates range from k1,1D,b8.2ns1 to k2,3D,b4.2ns1, this assembly process can be considered as reaction-limited ( < ). In the case of a reaction-limited process, the macroscopic off- and on-rates given in Eqs. (19) and (20) simplify to k˜i,j+ki,jD/ki,jD,bka=Vi,jka and k˜i,jkd. Thus, in this case the rates only depend on the microscopic reaction rates and the encounter volume Vi,j. In the right column of Fig. 9 (Figs. 9(b), 9(d), and 9(f)) a microscopic association rate of = 8 ns−1 is chosen. This 8-fold increase in the microscopic association rate is compensated by an 8-fold decrease in concentration. In this case the reaction is strongly affected by diffusion () and the macroscopic rates crucially depend on the different diffusive on- and off-rates ki,jD and ki,jD,b, which vary with the diffusive properties of the different fragments as can be seen in Table II and a clearly different assembly dynamic is expected. We will refer to this case as diffusion-influenced assembly.

        FIG. 9.

        Comparison of the rate equation approach Eq. (22) (solid lines) and the BD algorithm (dashed lines). The evolution of the normalized concentrations of the ring fragments c˜i=icic0 with = / are shown for different parameters. All simulations were started with = 500 individual proteins with a time resolution of Δ = 0.01 ns and the simulation results have been averaged over 40 independent trajectories. The left column corresponds to reaction-limited assembly with = 1 ns−1 and Vbox=(55nm)3. The right column corresponds to diffusion-influenced assembly with = 8 ns−1 and Vbox=(110nm)3 in which the increase in is balanced by the decrease in concentration. In (a) and (b), the result for a reversible reaction with a dissociation rate of = 0.0001 ns−1 and without internal bond formation is (kintraa=0.0 ns−1) shown. In (c) and (d), the final ring is additionally stabilized by the formation of an internal bond with rate kintraa=0.001 ns−1 (see Eq. (15)). In (e) and (f), the evolution of the cluster size distribution is shown for irreversible binding ( = 0 ns−1).

         

        In Figs. 9(a) and 9(b), the normalized fragment concentrations for the assembly with reversible bonds ( = 0.0001 ns−1) but without an additional stabilization of the full ring structure (kintraa=0) is shown. Comparing the predictions from the rate equation approach (Eq. (22)) (solid lines) with the averaged simulation results (dashed lines) we observe excellent agreement between our simulation results and the results from the macroscopic rate equation approach in the reaction-limited regime as well as in the strongly diffusion-influenced regime. The small deviations between the simulation results and the rate equation approach in Fig. 9(a) can be explained by our rough estimate for the excluded volume. This agreement between the two approaches, especially for the diffusion influenced regime is remarkable and demonstrates that, by evaluating the relevant diffusive properties of the fragments, we are indeed able to relate our microscopic reaction parameters to macroscopic reaction rates that correctly reproduce the dynamics of the system. Comparing the equilibrium distribution of the cluster fragments in Figs. 9(a) and 9(b), we see that indeed the 8-fold decrease in concentration is balanced by the 8-fold increase in concentration, as it was expected due to our previous reasoning. However, comparing the kinetics of assembly we see a clear difference between Figs. 9(a) and 9(b). In general, the approach towards equilibrium is slower in Figs. 9(b) than in 9(a). This shows that we are indeed in a different assembly regime and the change in concentration is not compensated by the change in in the kinetics of the assembly process. At lower concentration (Fig. 9(b)) clusters need a longer time to find each other diffusively which is reflected in the slower approach towards equilibrium.

        In Figs. 9(c) and 9(d), the normalized fragment concentrations with an additional stabilization of the final ring structure kintraa=10kd=0.001 ns−1 is shown. As discussed previously, the internal bond formation can be understood as an intrinsic process which stabilizes the full ring (Eq. (23)). The value chosen for the internal association rate corresponds to a free energy of ≈ −2.3 associated with the bond formation according to Eq. (15). In Figs. 9(c) and 9(d), the normalized fragment concentration of the full ring (black line) is shown irrespective of the number of bonds in the ring. Comparing the simulation results and the rate equation approach we again observe very good agreement between them. Similar to the case without internal bond formation (Figs. 9(a) and 9(b)), the assembly dynamics at lower concentration (Fig. 9(d)) is slower than at higher concentration (Fig. 9(c)), while the equilibrium steady state remains the same in both cases. In contrast to the reversible bond formation without additional stabilization of the ring (Figs. 9(c) and 9(e)) the normalized concentration of the full ring is significantly increased by the internal bond formation. Due to mass conservation the increase in the concentration of the full ring structure is balanced by a decrease in the concentration of smaller fragments showing how an additional internal state can alter the equilibrium cluster size distribution.

        Finally, we investigate the assembly dynamics for irreversible bond formation ( = 0 ns−1) for two different microscopic reaction rates corresponding to the regime of reaction-limited reactions (Fig. 9(e)) and diffusion-influenced reactions (Fig. 9(f)). Again the simulation results and the rate equation approach agree remarkably well with each other. By comparing Figs. 9(e) and 9(f), we again verify that the approach towards the steady state is much slower in the case of lower concentration shown in Fig. 9(f). However, in contrast to the previously discussed reversible bond formation, the steady state is different in both cases with different fraction of ring fragments containing 3, 4, and 5 proteins. This reflects a fundamental difference between the steady state observed for reversible bond formation and the steady state observed for irreversible bond formation. In the case of reversible bond formation the steady state is an equilibrium state in which dissociation and association events balance each other. This state only depends on the microscopic reaction rates and and the size of the reactive volume (Eq. (13)). In contrast the steady state for irreversible bond formation is a trapped steady state without any underlying dynamics. In this case of irreversible bond formation the ring fragments cannot disassemble and fragments of size ⩾ 5/2 cannot further assemble if all smaller fragments have been used up. Thus, in the case of trapping the final state depends on the kinetics of the assembly. The difference in the two steady states observed in Figs. 9(e) and 9(f) also shows that the assembly dynamics in Fig. 9(f) is not only slowed down compared to Fig. 9(e), but that the different diffusive properties of the cluster fragments lead to a change in the ratio of the concentration of intermediates during the assembly process. In general, larger fragments are diffusing slower than smaller fragments, and hence the rate for a reaction between two larger fragments is more strongly affected by diffusion than the reaction rate for two smaller fragments. In Fig. 9(f), this results in a lower concentration of pentameric rings, as the reaction probability of two small fragments is less affected by diffusion than the reaction probability of a smaller and a larger fragment. This leads to a stronger depletion of small fragments and hence a higher concentration of trapped ring fragments with 3 or 4 proteins.

         

        Biological structures like the actin cytoskeleton13–15 or the nuclear pore complex7–9 are highly dynamic with their constituents being continuously exchanged. To understand the biological functionality of these fascinating systems, a detailed understanding of the assembly dynamics is crucial. Moreover, advances in the fabrication of colloidal particles with directed interactions (“patchy particles”) allow to design a plethora of differently shaped building blocks with tunable interactions that can self-assemble into new materials.18–23 To increase the yield of the desired structures encoded in the local particle interactions, kinetic trapping in unfavorable configurations during the assembly process has to be prevented. Recently, it has been shown that particle interactions during the assembly process can be actively controlled.24,25 Thus, by a state- or time-dependent switching of reactivity the assembly process can be actively steered to reduce kinetic trapping. However, to fully exploit the possibilities of controlling the assembly process, a detailed understanding of its full dynamics with the formation of assembly intermediates is necessary.

        Here we presented a novel simulation approach which is ideally suited to investigate the dynamics of large protein assemblies with well-defined architectural properties. Proteins are represented by (multiple) non-overlapping, hard spheres equipped with reactive patches. In contrast to most previous studies on protein assembly,25,56–63,66,67 which rely on force fields to describe protein interactions, we combine overdamped Brownian motion with the concept of reversible, stochastic reactivity for patchy particles. Each cluster is treated as rigid object with its diffusive properties being evaluated on-the-fly. If the reactive patches of two clusters form an encounter by force-free diffusive motion, a bond between the clusters is established with a predefined rate. Local rules are used to describe the resulting rearrangements, which are assumed to be very fast on the time scale of the BD simulations. For potential-based simulations, these rearrangements would proceed due to the corresponding forces, which do not exist in our approach. Similar to the association process, the dissociation process also occurs with a predefined rate. Here we have derived the rules for the placement of the two partners which are required to satisfy detailed balance.

        Together, these rules render our approach very efficient because complexes formed during the assembly are diffusing as rigid objects (thus reducing the number of degrees of freedom to be propagated to six) and because interactions are treated locally (as opposed to evaluating potentials). Nevertheless, the patchy particle approach includes detailed information on the architecture of protein assemblies, thus allowing us to study the spatial-temporal dynamics of specific biological systems of interest. Moreover, we have been able to show how the microscopic reaction parameters used in our approach correspond to the macroscopic reaction rates commonly used to describe the dynamics of the concentration of assembly intermediates. This allows a direct comparison between macroscopic, experimentally measurable concentrations and simulation results. To the best of our knowledge, this is the first time that a BD algorithm combining stochastic reactivity of anisotropic patches with reversible dynamics that fulfills detailed balance is presented. Testing our algorithm against mean field prediction for the case of bimolecular reactions revealed the importance of satisfying detailed balance.

        As an example for a multi-component cluster, we investigated the assembly of a pentameric ring. Here we compared our simulation results to a rate equation approach for the different fragment concentrations. By extracting the diffusive on- and off-rates from our simulations we find excellent agreement between the rate equation approach and our simulation results even for the case of strongly diffusion-influenced assembly without any free parameter. We showed that the assembly dynamics can be significantly changed by the change in the diffusive properties of the assembly intermediates. In the case of irreversible bond formation the system becomes trapped with incompatible ring fragments. In this case not only the assembly kinetics but also the finally reached steady state depends on the different diffusive properties of the assembly intermediates. For complex assembly geometries trapping is a common motif. Thus, in this case not only correctly predicted equilibrium structure distribution of the assembly intermediates is relevant but also the correctly predicted kinetics of the assembly process as the system might become trapped before reaching the expected equilibrium steady state. By demonstrating how the macroscopic reaction rates can be calculated from our simulation, we have developed a framework in which a qualitative prediction of the changes in the macroscopic reaction rates, based on the changes of the diffusive properties of the assembly intermediates, is possible.

        Because our approach aims to describe the assembly of large protein structures, the details of individual bond formations are not resolved within this framework. In general, the use of local rules requires well-defined target geometries and the instantaneous local rearrangement accompanying stochastic bond formation follows from the assumption that the time scales of assembly and local rearrangements are well separated. Given the coarse-grained nature of our approach, the instantaneous rearrangement during binding is a valid approximation if the encounter volume specifies a narrow region around the desired configuration. In the case of large encounter volumes the assumption of fast and small local rearrangements can break down and the use of local rules needs to be considered with care. In particular, in the case of higher protein densities or assembly close to a membrane this can lead to a situation in which physically reasonable rearrangements are suppressed by our steric rules. In classical MD or BD simulations with potentials, these rearrangements would be realized due to ensuing mechanical forces. Typical biological examples for situations which are out of the scope of our current approach are the bending of a membrane by BAR-proteins,81,82 the emergence of polymorphic virus structures83–85 or the assembly of actin networks from preexisting large fibers.86 In principle, such a situation could be resolved by using hybrid schemes that interface our approach with potential-based simulations for local rearrangements. For the time being, however, we focus on the assembly of protein complexes with well-defined architectures, such as native viruses,87 centrioles,6,88–90 or actin filaments in networks growing by incorporation of actin monomers from solution.14,86,91

        In the future, our approach can be used to study the assembly of such complex protein clusters. After we have successfully verified our approach by comparing averaged simulation data with macroscopic mean field results, we now can investigate the stochastic variance inherent to these self-assembly processes. Our approach is ideally suited to study the effect of time- or state-dependent changes in reactivity, as has been recently shown in a qualitative study on the effect on hierarchical assembly in virus capsids.26 Moreover, this approach might also help to design artificial systems that do not get kinetically trapped in undesired structures. It can also be coupled with hydrodynamic schemes, in particular with reactive multi-particle collision dynamics.92

         

        H.C.R.K. acknowledges financial support by the Cusanuswerk. U.S.S. is a member of the Heidelberg cluster of excellence CellNetworks. We thank Nils Becker and Joris Paijmans for helpful discussions.

         
        1. R. R. Gabdoulline and R. C. Wade,   72, 1917 (1997).http://dx.doi.org/10.1016/S0006-3495(97)78838-6
        2. W. Roos, R. Bruinsma, and G. Wuite,   6, 733 (2010).http://dx.doi.org/10.1038/nphys1797
        3. A. Zlotnick and S. Mukhopadhyay,   19, 14 (2011).http://dx.doi.org/10.1016/j.tim.2010.11.003
        4. D. A. Compton,   69, 95 (2000).http://dx.doi.org/10.1146/annurev.biochem.69.1.95
        5. E. Karsenti, F. Nédélec, and T. Surrey,   8, 1204 (2006).http://dx.doi.org/10.1038/ncb1498
        6. P. Gönczy,   13, 425 (2012).http://dx.doi.org/10.1038/nrm3373
        7. F. Alber, S. Dokudovskaya, L. M. Veenhoff, W. Zhang, J. Kipper, D. Devos, A. Suprapto, O. Karni-Schmidt, R. Williams, B. T. Chait ,   450, 683 (2007).http://dx.doi.org/10.1038/nature06404
        8. M. A. D’Angelo and M. W. Hetzer,   18, 456 (2008).http://dx.doi.org/10.1016/j.tcb.2008.07.009
        9. E. J. Tran and S. R. Wente,   125, 1041 (2006).http://dx.doi.org/10.1016/j.cell.2006.05.027
        10. B. Geiger and A. Bershadsky,   13, 584 (2001).http://dx.doi.org/10.1016/S0955-0674(00)00255-6
        11. A. D. Bershadsky, C. Ballestrem, L. Carramusa, Y. Zilberman, B. Gilquin, S. Khochbin, A. Y. Alexandrova, A. B. Verkhovsky, T. Shemesh, and M. M. Kozlov,   85, 165 (2006).http://dx.doi.org/10.1016/j.ejcb.2005.11.001
        12. H. Wolfenson, Y. I. Henis, B. Geiger, and A. D. Bershadsky,   66, 1017 (2009).http://dx.doi.org/10.1002/cm.20410
        13. C. C. Beltzner and T. D. Pollard,   283, 7135 (2007).http://dx.doi.org/10.1074/jbc.M705894200
        14. K. Guo, J. Shillcock, and R. Lipowsky,   131, 015102 (2009).http://dx.doi.org/10.1063/1.3159003
        15. T. D. Pollard,   36, 451 (2007).http://dx.doi.org/10.1146/annurev.biophys.35.040405.101936
        16. K. Guo, J. Shillcock, and R. Lipowsky,   133, 155105 (2010).http://dx.doi.org/10.1063/1.3497001
        17. G. Schreiber, G. Haran, and H.-X. Zhou,   109, 839 (2009).http://dx.doi.org/10.1021/cr800373w
        18. C. A. Mirkin, R. L. Letsinger, R. C. Mucic, and J. J. Storhoff,   382, 607 (1996).http://dx.doi.org/10.1038/382607a0
        19. L. Di Michele and E. Eiser,   15, 3115 (2013).http://dx.doi.org/10.1039/c3cp43841d
        20. Y. Wang, Y. Wang, D. R. Breed, V. N. Manoharan, L. Feng, A. D. Hollingsworth, M. Weck, and D. J. Pine,   491, 51 (2012).http://dx.doi.org/10.1038/nature11564
        21. C. Knorowski and A. Travesset,   15, 262 (2011).http://dx.doi.org/10.1016/j.cossms.2011.07.002
        22. S. Sacanna, W. Irvine, P. M. Chaikin, and D. J. Pine,   464, 575 (2010).http://dx.doi.org/10.1038/nature08906
        23. S. Sacanna, M. Korpics, K. Rodriguez, L. Colón-Meléndez, S.-H. Kim, D. J. Pine, and G.-R. Yi,   4, 1688 (2013).http://dx.doi.org/10.1038/ncomms2694
        24. M. E. Leunissen, R. Dreyfus, F. C. Cheong, D. G. Grier, R. Sha, N. C. Seeman, and P. M. Chaikin,   8, 590 (2009).http://dx.doi.org/10.1038/nmat2471
        25. L. Di Michele, F. Varrato, J. Kotar, S. H. Nathan, G. Foffi, and E. Eiser,   4, 2007 (2013).http://dx.doi.org/10.1038/ncomms3007
        26. J. E. Baschek, H. C. Klein, and U. S. Schwarz,   5, 22 (2012).http://dx.doi.org/10.1186/2046-1682-5-22
        27. M. v. Smoluchowski,   92, 9 (1917).
        28. F. C. Collins and G. E. Kimball,   4, 425 (1949).http://dx.doi.org/10.1016/0095-8522(49)90023-9
        29. K. Šolc and W. Stockmayer,   5, 733 (1973).http://dx.doi.org/10.1002/kin.550050503
        30. D. Shoup, G. Lipari, and A. Szabo,   36, 697 (1981).http://dx.doi.org/10.1016/S0006-3495(81)84759-5
        31. O. Berg,   47, 1 (1985).http://dx.doi.org/10.1016/S0006-3495(85)83870-4
        32. A. Shushin,   110, 12044 (1999).http://dx.doi.org/10.1063/1.479140
        33. M. Schlosshauer and D. Baker,   106, 12079 (2002).http://dx.doi.org/10.1021/jp025894j
        34. M. Schlosshauer and D. Baker,   13, 1660 (2004).http://dx.doi.org/10.1110/ps.03517304
        35. N. Agmon,   81, 2811 (1984).http://dx.doi.org/10.1063/1.447954
        36. N. Agmon and A. Szabo,   92, 5270 (1990).http://dx.doi.org/10.1063/1.458533
        37. S. H. Northrup, S. A. Allison, and J. A. McCammon,   80, 1517 (1984).http://dx.doi.org/10.1063/1.446900
        38. S. H. Northrup and H. P. Erickson,   89, 3338 (1992).http://dx.doi.org/10.1073/pnas.89.8.3338
        39. H. X. Zhou,   94, 8794 (1990).http://dx.doi.org/10.1021/j100388a010
        40. H.-X. Zhou,   64, 1711 (1993).http://dx.doi.org/10.1016/S0006-3495(93)81543-1
        41. R. R. Gabdoulline and R. C. Wade,   306, 1139 (2001).http://dx.doi.org/10.1006/jmbi.2000.4404
        42. G. Zou and R. D. Skeel,   85, 2147 (2003).http://dx.doi.org/10.1016/S0006-3495(03)74641-4
        43. R. Alsallaq and H.-X. Zhou,   15, 215 (2007).http://dx.doi.org/10.1016/j.str.2007.01.005
        44. R. Alsallaq and H.-X. Zhou,   92, 1486 (2007).http://dx.doi.org/10.1529/biophysj.106.096024
        45. S. Qin, X. Pang, and H.-X. Zhou,   19, 1744 (2011).http://dx.doi.org/10.1016/j.str.2011.10.015
        46. M. Klann and H. Koeppl,   13, 7798 (2012).http://dx.doi.org/10.3390/ijms13067798
        47. S. S. Andrews and D. Bray,   1, 137 (2004).http://dx.doi.org/10.1088/1478-3967/1/3/001
        48. S. S. Andrews,   2, 111 (2005).http://dx.doi.org/10.1088/1478-3975/2/2/004
        49. J. S. van Zon and P. R. Ten Wolde,   123, 234910 (2005).http://dx.doi.org/10.1063/1.2137716
        50. M. J. Morelli and P. R. Ten Wolde,   129, 054112 (2008).http://dx.doi.org/10.1063/1.2958287
        51. K. Takahashi, S. Tănase-Nicola, and P. R. ten Wolde,   107, 2473 (2010).http://dx.doi.org/10.1073/pnas.0906885107
        52. D. T. Gillespie,   81, 2340 (1977).http://dx.doi.org/10.1021/j100540a008
        53. J. Hattne, D. Fange, and J. Elf,   21, 2923 (2005).http://dx.doi.org/10.1093/bioinformatics/bti431
        54. D. Fange, O. G. Berg, P. Sjöberg, and J. Elf,   107, 19820 (2010).http://dx.doi.org/10.1073/pnas.1006565107
        55. D. Fange, A. Mahmutovic, and J. Elf,   28, 3155 (2012).http://dx.doi.org/10.1093/bioinformatics/bts584
        56. D. Rapaport,   101, 186101 (2008).http://dx.doi.org/10.1103/PhysRevLett.101.186101
        57. M. F. Hagan and D. Chandler,   91, 42 (2006).http://dx.doi.org/10.1529/biophysj.105.076851
        58. H. D. Nguyen, V. S. Reddy, and C. L. Brooks III,   7, 338 (2007).http://dx.doi.org/10.1021/nl062449h
        59. D. Rapaport,   86, 051917 (2012).http://dx.doi.org/10.1103/PhysRevE.86.051917
        60. I. G. Johnston, A. A. Louis, and J. P. Doye,   22, 104101 (2010).http://dx.doi.org/10.1088/0953-8984/22/10/104101
        61. C. Horejs, M. K. Mitra, D. Pum, U. B. Sleytr, and M. Muthukumar,   134, 125103 (2011).http://dx.doi.org/10.1063/1.3565457
        62. A. W. Wilber, J. P. Doye, A. A. Louis, E. G. Noya, M. A. Miller, and P. Wong,   127, 085106 (2007).http://dx.doi.org/10.1063/1.2759922
        63. A. W. Wilber, J. P. Doye, and A. A. Louis,   131, 175101 (2009).http://dx.doi.org/10.1063/1.3243580
        64. J. Largo, F. W. Starr, and F. Sciortino,   23, 5896 (2007).http://dx.doi.org/10.1021/la063036z
        65. P. F. Damasceno, M. Engel, and S. C. Glotzer,   337, 453 (2012).http://dx.doi.org/10.1126/science.1220869
        66. X. Liu, J. C. Crocker, and T. Sinno,   138, 244111 (2013).http://dx.doi.org/10.1063/1.4811656
        67. F. Sciortino, C. De Michele, and J. F. Douglas,   20, 155101 (2008).http://dx.doi.org/10.1088/0953-8984/20/15/155101
        68. J. D. Halverson and A. V. Tkachenko,   87, 062310 (2013).http://dx.doi.org/10.1103/PhysRevE.87.062310
        69. M. Eigen, “Diffusion control in biochemical reactions,” in   , Studies in the Natural Sciences Vol. 4, edited by S. L. Mintz, and S. M. Widmayer (Plenum Press, New York, 1974), pp. 37–61.
        70. D. Shoup and A. Szabo,   40, 33 (1982).http://dx.doi.org/10.1016/S0006-3495(82)84455-X
        71. C. Korn and U. Schwarz,   126, 095103 (2007).http://dx.doi.org/10.1063/1.2464080
        72. J. Schluttig, D. Alamanova, V. Helms, and U. S. Schwarz,   129, 155106 (2008).http://dx.doi.org/10.1063/1.2996082
        73. J. Schluttig, C. B. Korn, and U. S. Schwarz,   81, 030902 (2010).http://dx.doi.org/10.1103/PhysRevE.81.030902
        74. J. García de la Torre, M. L. Huertas, and B. Carrasco,   78, 719 (2000).http://dx.doi.org/10.1016/S0006-3495(00)76630-6
        75. B. Berger, P. W. Shor, L. Tucker-Kellogg, and J. King,   91, 7732 (1994).http://dx.doi.org/10.1073/pnas.91.16.7732
        76. T. L. Hill,   (Dover Publications, New York, 1986) (Unabridged, corr. republ. of the 2. print., Reading, MA, 1962).
        77. M. Pogson, R. Smallwood, E. Qwarnstrom, and M. Holcombe,   85, 37 (2006).http://dx.doi.org/10.1016/j.biosystems.2006.02.004
        78. M. T. Klann, A. Lapin, and M. Reuss,   5, 71 (2011).http://dx.doi.org/10.1186/1752-0509-5-71
        79. J. Paijmans, “The fundamental lower bound of the noise in transcriptional regulation,” Master's thesis, University of Amsterdam, Amsterdam (July 2012).
        80. N. Kern and D. Frenkel,   118, 9882 (2003).http://dx.doi.org/10.1063/1.1569473
        81. A. Arkhipov, Y. Yin, and K. Schulten,   95, 2806 (2008).http://dx.doi.org/10.1529/biophysj.108.132563
        82. A. Frost, V. M. Unger, and P. D. Camilli,   137, 191 (2009).http://dx.doi.org/10.1016/j.cell.2009.04.010
        83. H. D. Nguyen and C. L. Brooks III,   8, 4574 (2008).http://dx.doi.org/10.1021/nl802828v
        84. O. M. Elrad and M. F. Hagan,   8, 3850 (2008).http://dx.doi.org/10.1021/nl802269a
        85. H. D. Nguyen, V. S. Reddy, and C. L. Brooks III,   131, 2606 (2009).http://dx.doi.org/10.1021/ja807730x
        86. L. Blanchoin, R. Boujemaa-Paterski, C. Sykes, and J. Plastino,   94, 235 (2014).http://dx.doi.org/10.1152/physrev.00018.2013
        87. D. L. D. Caspar and A. Klug,   (Cold Spring Harbor Laboratory Press, 1962), Vol. 27, pp. 49–50.
        88. D. Kitagawa, I. Vakonakis, N. Olieric, M. Hilbert, D. Keller, V. Olieric, M. Bortfeld, M. C. Erat, I. Flückiger, P. Gönczy, and M. O. Steinmetz,   144, 364 (2011).http://dx.doi.org/10.1016/j.cell.2011.01.008
        89. M. van Breugel, M. Hirono, A. Andreeva, H.-a. Yanagisawa, S. Yamaguchi, Y. Nakazawa, N. Morgner, M. Petrovich, I.-O. Ebong, C. V. Robinson, C. M. Johnson, D. Veprintsev, and B. Zuber,   331, 1196 (2011).http://dx.doi.org/10.1126/science.1199325
        90. P. Guichard, V. Hachet, N. Majubu, A. Neves, D. Demurtas, N. Olieric, I. Fluckiger, A. Yamada, K. Kihara, Y. Nishida, S. Moriya, M. O. Steinmetz, Y. Hongoh, and P. Gönczy,   23, 1620 (2013).http://dx.doi.org/10.1016/j.cub.2013.06.061
        91. C. Erlenkämper and K. Kruse,   139, 164907 (2013).http://dx.doi.org/10.1063/1.4825248
        92. R. Kapral,   138, 020901 (2013).http://dx.doi.org/10.1063/1.4773981
         
        true
        true
        This is a required field
        Please enter a valid email address

        Oops! This section does not exist...

        Use the links on this page to find existing content.

        752b84549af89a08dbdd7fdb8b9568b5 journal.articlezxybnytfddd
        Scitation: Studying protein assembly with reversible Brownian dynamics of patchy particles
        http://aip.metastore.ingenta.com/content/aip/journal/jcp/140/18/10.1063/1.4873708
        10.1063/1.4873708
        SEARCH_EXPAND_ITEM