Studying protein assembly with reversible Brownian dynamics of patchy particles
SStudying Protein Assembly with Reversible Brownian Dynamicsof Patchy Particles
Heinrich C. R. Klein
Institute for Theoretical Physics, Heidelberg University, 69120 Heidelberg, Germany
Ulrich S. Schwarz ∗ Institute for Theoretical Physics, Heidelberg University, 69120 Heidelberg, Germany andBioQuant, Heidelberg University, 69120 Heidelberg, Germany (Dated: November 14, 2018)Assembly of protein complexes like virus shells, the centriole, the nuclear porecomplex or the actin cytoskeleton is strongly determined by their spatial structure.Moreover it is becoming increasingly clear that the reversible nature of protein assem-bly is also an essential element for their biological function. Here we introduce a com-putational approach for the Brownian dynamics of patchy particles with anisotropicassemblies and fully reversible reactions. Different particles stochastically associateand dissociate with microscopic reaction rates depending on their relative spatial po-sitions. The translational and rotational diffusive properties of all protein complexesare evaluated on-the-fly. Because we focus on reversible assembly, we introduce ascheme which ensures detailed balance for patchy particles. We then show how themacroscopic rates follow from the microscopic ones. As an instructive example, westudy the assembly of a pentameric ring structure, for which we find excellent agree-ment between simulation results and a macroscopic kinetic description without anyadjustable parameters. This demonstrates that our approach correctly accounts forboth the diffusive and reactive processes involved in protein assembly. ∗ [email protected] a r X i v : . [ q - b i o . S C ] M a y I. INTRODUCTION
Assembly of biomolecules into supramolecular complexes is at the heart of many bio-logical processes and the dynamic interplay of the different components leads to biologicalfunctionality. The most important type of self-assembling biomolecules are proteins. Al-though they often have a globular shape, protein assemblies can be strongly anisotropicdue to the localized binding interactions between the proteins. Like biological systems ingeneral, protein assemblies are operative on many different scales. Their size ranges fromthe nanometer-scale (for example for two-component complexes like Barnase and Barstar[1]) through tens of nanometers (for example for viruses, which typically consist of smallmultiples of 60 components [2, 3]) up to the micrometer scale (for example the mitotic spin-dle [4–6]). Even in steady state most biological complexes remain highly dynamic, withassociation events being balanced by dissociation events. Prominent examples are nuclearpore complexes [7–9] or focal adhesions[10–12]. Another important example for the dynamicnature of biological complexes are actin filaments [13–15], which are called living polymers 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 < M − s − up to > M − s − . [17] The highly dynamic nature of biological assem-blies as well as the large range of involved spatial and temporal scales renders this physicalproblem challenging but fascinating.Interestingly, recent advances in the fabrication of functionalized colloids with directionalinteractions ( patchy particles ) make it possible to design elementary building blocks of mi-crometer sizes which can be used to self-assemble complexes with new functionality. Exper-imental techniques ranging from DNA-mediated self-assembly [18–21] to entropic depletioninteractions [22, 23] have been rapidly advancing during the last decade, thus providing aplethora of possibilities to fabricate colloidal particles whose shape and interactions can becontrolled in detail. Moreover external stimuli such as temperature, light or pH can be usedto control the inter-particle interactions during the assembly process . [24, 25] These tech-niques allow for a state- or time-dependent switching of the interactions and can be usedto steer the assembly process. Controlling the particle interactions during the assemblyprocess can prevent kinetic trapping resulting in a higher yield of the desired structure, ashas recently been shown in a computational study for virus assembly. [26] To fully exploitthe potential of these techniques, a detailed understanding of the dynamics of the assemblyprocess (distribution of intermediates, relevant time scales) is of crucial importance.Understanding the mechanisms governing chemical reactions has a long tradition in the-oretical physics and chemistry. The most powerful analytical technique in this context is theFokker-Planck or Smoluchowski equation, which has been used early to study bimolecularassociation based on diffusive motion. This approach was pioneered by Smoluchowski whofirst calculated the maximum diffusion-limited reaction rate for a fully reactive sphericalparticle . [27] An important generalization of this result was derived by Collins and Kimball[28] who studied the effect of finite reactivity by introducing a radiation boundary conditionrelating the concentration at contact to the reactive flux. Another essential extension isthe case of particles with anisotropic reactivity ( patchy particles ) . [29–34] However, mostof these generalizations focus on the calculation of association rates and only a few aim atdescribing reversible reactions. [35, 36]To study the full dynamics of protein assembly one needs to consider both associationand dissociation processes. Moreover, even for globular proteins the intermediates formedduring the assembly process are often of highly non-spherical geometry, thus rendering ananalytical treatment of assembly in the framework of the Fokker-Planck equation very dif-ficult. In this case computer simulations provide a valuable alternative. While detailedmolecular dynamics (MD) simulations have been successfully used to investigate the behav-ior of a single molecule in great detail, studying the dynamics of large protein complexes onbiologically relevant time and length scales is prohibited by the high computational costsof these simulations. Therefore coarse-grained models are required for this case. Browniandynamics (BD) simulations, the numerical counterpart to the Fokker-Planck equation, havebeen 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 associationrates in bimolecular reactions based on realistic protein shapes and binding interactions. Tostudy reaction dynamics in a large system consisting of many proteins, various simulationframeworks have been developed. [46] They range from space- and time-continuous BD sim-ulations as for example in the Smoldyn framework [47, 48] through event-driven tools likeGreens functions reaction dynamics (GFRD) [49–51], which is based on the analytical solu-tion of the Fokker-Planck equation, up to a space-discretized version of Gillespie’s chemicalmaster equation approach [52], as for example used in MesoRD [53–55]. While these simula-tion frameworks have been successfully used to study reaction kinetics on a large scale, theycannot be used to study the details of assembly processes as all of these frameworks lack adetailed 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 forwhich these elements are essential, coarse-grained MD simulations [56–59] and Monte Carlo(MC) techniques [60–63] have been used. Similarly the interaction of colloidal particles hasbeen investigated on various scales with different simulation techniques [21], including MDsimulations [64], MC simulations [65, 66] and BD studies [67, 68].Here we introduce a simulation framework which allows us to study the spatial andstochastic aspects of protein assembly to large detail and nevertheless is computationallyrelatively cheap. In our approach, we combine BD of realistic protein shapes with stochas-tic reactivity for both, association and dissociation. Proteins are described as assembliesof spheres equipped with reactive patches. Their motion in real and orientation space isbased on the overdamped Langevin equation with an anisotropic diffusion tensor. Inspiredby the notion of an encounter complex [1, 17, 28, 36, 43, 69, 70], we decompose the reactionprocess 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 firstshow how these microscopic reaction rates can be related to macroscopic, experimentallymeasurable rates. We then verify that our algorithm correctly reproduces the macroscop-ically expected equilibrium behavior for bimolecular reactions. Finally we investigate theassembly of a pentameric ring structure and compare our simulation results to a macroscopicrate equation approach. Here we again find excellent agreement between the stochastic sim-ulations and the macroscopic description if the macroscopic rates are calculated without anyfree 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 asrealistic diffusion properties to understand the dynamics governing protein assembly. FIG. 1. Model definition. a) A dumbbell-shaped protein modeled by two hard spheres. Theprotein is covered with a patch of radius R p whose center coincides with the center of one of thesteric spheres and is thus shifted by (cid:126)c p relative to the center of the protein (contact point of thespheres). The patch has an opening angle of θ p around the orientation vector (cid:126)o pointing alongthe long axis of the dumbbell. b) The encounter between two of the proteins of a). Here (cid:126)r is thecenter-to-center(ctc)-vector between the proteins and (cid:126)r p is the center-to-center vector between thepatches. c) Illustration of the reorientation processes upon reaction. In a first step the ctc vectorsof the proteins are aligned with the desired ctc-vectors. In a next step the projection of the torsionvectors (cid:126)t and (cid:126)t on a plane perpendicular to (cid:126)r ctc are aligned. Finally the clusters are shifted sothat the desired relative distance between the proteins is achieved. II. MODEL AND METHODSA. 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: theundirected diffusive motion of the two binding partners A and B until they stochasticallyreach the encounter state A · B and the reaction from the encounter state to the boundcomplex C . In terms of the free energy landscape the encounter complex represents abarrier separating the flat diffusive energy landscape from the reaction funnel. Thus, itis a transient state which can be characterized by the onset of highly correlated motionbetween the binding partners. [43] A barrier in the free energy landscape might arise forexample due to a necessary rearrangement of the binding site or due to a reorganization ofthe water layer. In some cases there might not exist an explicit barrier, but even then theencounter complex is a helpful theoretical concept. Considering the encounter complex asthe “watershed” between diffusion and reaction processes, the bimolecular reaction can besplit according to the following reaction scheme: A + B k + (cid:10) k − C (1) A + B k D (cid:10) k D,b A · B k a (cid:10) k d C . (2)Here k + and k − represent the overall kinetic rates of the reaction. In Eq. 2 the bindingand unbinding process is split into a diffusive part represented by the rates k D and k D,b and a reaction part represented by k a and k d , respectively. Separating the binding processin free diffusive motion without steering forces and localized binding requires short-rangedinteractions. Thus, our approach for the simulation of protein association corresponds toa regime of high salt concentration screening long-ranged electrostatic forces. Indeed thisis a reasonable assumption for physiological salt conditions. Moreover, our approach canbe used to study the assembly of µ m-sized colloids functionalized with reactive patcheswhere the size of the reactive region is small compared to the particle size. As we considera situation in which the interaction between the proteins is short-ranged and only affectsthe reaction probabilities in the encounter state, individual clusters undergo free diffusivemotion without an additional drift term. Moreover, here we do not consider hydrodynamicinteractions 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 multipleproteins held together by bonds. Two clusters can react with each other by bond formationbetween adhesive patches. A cluster consisting of multiple proteins can decay into smallerclusters by bond dissociation. Proteins are described by sets of non-overlapping, hard spheresapproximating their shape. The proteins are equipped with reactive patches which reflectthe localized binding sites. A patch is defined starting from a sphere of radius R p whosecenter (cid:126)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 (cid:126)o . As a simple example for anon-globular protein, in Fig. 1a we schematically show a dumbbell-shaped protein consistingof two hard spheres in contact. The protein is equipped with one reactive patch of radius R p whose center coincides with the center of one of the steric spheres and is thus shifted by (cid:126)c p from the protein center located at the contact point of the steric spheres. The orientationvector (cid:126)o points along the long axis of the dumbbell. Using the above definitions we formulatethe following set of equations to define the encounter between two clusters mediated by apair of patches, compare Fig. 1b: r p ≤ R p, + R p, (3)acos (cid:0) (cid:126)r · (cid:126)o | (cid:126)r || (cid:126)o | (cid:1) ≤ θ p, , acos (cid:0) − (cid:126)r · (cid:126)o | (cid:126)r || (cid:126)o | (cid:1) ≤ θ p, . (4)Eq. 3 means that the distance r p between the centers of the two patches has to be sufficientlysmall for an encounter to occur and can be implemented easily in a simulation. Eq. 4 ismore complex. It involves the center-to-center vector (cid:126)r and not only reflects that the patchesare anisotropic, but also assumes that an encounter only occurs if the two partners have afavorable orientation relatively to each other which is defined by the two parameters θ p, and θ p, . In general an anisotropic protein can be described by multiple spheres and thecenter of the patches can be chosen independent of the steric spheres, allowing for a variabledescription of proteins and their localized interactions. C. Particle motion
In our approach particle motion is described by the six-dimensional, overdamped Langevinequation (Brownian motion) describing translational and rotational diffusion of an arbitrar-ily shaped but rigid object [71, 72]: ∂ t X t = g t with: (cid:104) g t (cid:105) = 0, (cid:104) g t g t (cid:48) (cid:105) = 2 k B T M δ ( t − t (cid:48) ). (5) P (r) / ( (cid:47) r ) r (nm)Remain at old positionChoose new position FIG. 2. Rejection methods. The distance distribution of two hard spheres of radius R = 1nm ina periodic box of volume V box = (8nm) is shown for two different rejection methods. In the firstrejection method (green) new positions of the particles are chosen (starting from their positionsbefore collision) until no overlap is observed. In the second rejection method (red) the spheresare set back to their position before the collision. The simulations were performed at a time-stepresolution of ∆ t = 0 . Here X t is a six-dimensional position vector describing the position and orientations at timet and g t is Gaussian white noise. M is the mobility matrix which is calculated on-the-fly forany cluster shape [71–73] following the method of de la Torre [74]. Because intermediatescontinuously grow and shrink due to association and dissociation, the diffusive properties ofthe 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 schemeis not appropriate. Here we treat steric collisions by a rejection method similar to MCsimulations. If a steric overlap between two clusters is created during a move step, this movestep is rejected. One method to deal with this situation would be to repeatedly choose a newposition of the two clusters, starting from their original location and orientation, until themove step is accepted. This method leads to incorrect simulation results as is demonstratedfor the example of two hard spheres in a periodic box in Fig. 2. Here a histogram of thedistance distribution between the two spheres, scaled by the volume of a spherical shell, isshown in a range in which boundary effects are not observed. By choosing a new positionuntil successful acceptance of the move step (green histogram), an effective repulsion betweenthe two spheres is introduced, leading to a depletion zone at small particle distances insteadof the expected uniform distribution. If, instead of choosing a new position, we set the twoclusters involved in a steric collision back to their original position (and orientation) beforethe move step, we recover the expected uniform equilibrium distribution (shown in red inFig. 2). Thus, this method ensures that the correct equilibrium distribution is realized. D. Local rules
In our approach, clusters are rigid assemblies consisting of one or multiple proteins. As-suming that the assembly geometry is determined by unique local interactions, a set oflocal rules is necessary to describe the new relative orientation and position of two reactingclusters.[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 twoproteins directly involved in the binding process as well as their relative orientation (tor-sion vectors). Upon reaction the two clusters instantaneously flip into the desired relativeconfiguration. This rule follows from the assumption that the short-ranged forces involvedin the final binding step lead to a fast rearrangement on the time scale of our BD simula-tion. The reorientation process is schematically shown in Fig. 1c. In a first step the clustersare shifted and rotated so that the predefined ctc-vector matches the real ctc-vector. Thisprocedure 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 projectionof the two torsion vectors on a plane perpendicular to the ctc-vector. The necessary rota-tion and translation of the clusters is distributed between them according to their diffusiveweights. This means that for a small and a large partner, essentially only the small partner ismoved, as one expects for physical reasons; for two similarly sized partners, both are movedto a similar extent. Since all interactions are local (e.g. local patches or constraints onorientation/torsion), the encounter complex corresponds to a region in configuration spacesurrounding the desired relative position and orientation. If this region is sufficiently small,the instantaneous flipping into the correct configuration is comparable to the resolution ofthe BD simulation. Reorientation of the clusters during the binding process can lead toa steric overlap either with another cluster or between the merging clusters. If this is the0case, the reaction is rejected and the old positions and orientations of the clusters before themove step are resumed. The steric collision of two merging clusters reflects that during theassembly 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 fiveproteins, 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 stochasticreactions with two corresponding microscopic rates. If the reactive patches of two clustersform an encounter, specified by the constraints in Eq. 3 and Eq. 4, they can react with abond-specific rate k (cid:48) a . Similarly an existing bond can dissociate with a bond-specific rate k (cid:48) d . Here we assume that the waiting times for association or dissociation of a bond betweentwo clusters in encounter are Poisson-distributed, P ( t, k (cid:48) i ) = k (cid:48) i exp( − k (cid:48) i t ), i = { a, d } . In thiscase the probability that no association or dissociation has occurred after a timestep ∆ t isgiven by S (∆ t, k (cid:48) i ) = 1 − (cid:82) ∆ t P ( t, k (cid:48) i ) dt and hence the probabilities for bond dissociation orbond formation are given by: [50] P assoc = 1 − S (∆ t, k (cid:48) a ) ≈ k (cid:48) a ∆ t , P dissoc = 1 − S (∆ t, k (cid:48) d ) ≈ k (cid:48) d ∆ t . (6)The approximations used in Eq.6 are valid if k (cid:48) d ∆ t (cid:28) k (cid:48) a ∆ t (cid:28)
1, which is the casethroughout this work.We now show how the microscopic reaction rates k (cid:48) a and k (cid:48) d can be related to macroscopicreaction rates. In a macroscopic framework the reaction scheme in Eq. 2 can be interpretedas a system of ordinary differential equations describing the changes in the concentrations c A , c B , c A · B and c C of the different species involved in the reaction . [69] In this frameworkthe encounter is understood as a single, intermediate state connecting the bound complex C and the unbound A and B particles. The encounter can either react to the final complexwith the first order rate k a or decay into two separated particles with the first order rate k D,b . An encounter can be formed by the first order decay of cluster C with a rate k d or bythe diffusive association of an A and B particle described by the second order rate k D . Here k D has the dimension m / s if we consider particle concentrations in 3 dimensions. Using asteady state approximation of the encounter complex, the following relations are obtained1 FIG. 3. Illustration of the reactive volume V (cid:63) . The configuration space is defined by the relativeposition (cid:126)r of the particles and the 2 × (cid:126) Ω. a) In the case ofone reactive patch V (cid:63) corresponds to the region in configuration space in which two particles areconsidered to be in encounter (Eq. 3 and Eq. 4). b) For particles with multiple reactive patchesdifferent regions V (cid:63)i,j in configuration space correspond to an encounter mediated by one particularpatch combination i , j . For identical microscopic rates k a and non-overlapping V (cid:63)i,j the total reactivevolume V (cid:63) between the two particles is given by the sum over the encounter volumes associatedwith each patch combination. k + = k D k a k a + k D,b , k − = k D,b k d k a + k D,b (7) ⇒ K eq = k + k − = k D k D,b k a k d . (8)In the case of a diffusion-limited reaction k a (cid:29) k D,b . [69] In this case the overall forwardand reverse rates simplify into the purely diffusion-limited forward rate k + ≈ k D and thebackward rate k − becomes k − ≈ k D,b k d /k a . In the case that the reaction is reaction-limited( k D,b (cid:29) k a ), the overall forward rate becomes k + ≈ k D k a /k D,b and the overall backward rateis k − ≈ k d .In the framework of the Fokker-Planck equation two spherical particles are considered tobe in an encounter if they are in contact. [27, 28] Finite reactivity for particles in encounterwas introduced by Collins and Kimball [28] in the form of a radiation boundary conditionin which the rate κ a (termed κ in Ref. [70] and k in Ref. [28]) relates the concentrationat contact to the reactive flux. By comparing the escape probabilities from the encounterin the Fokker-Planck framework and the macroscopic rate equation framework, Shoup andSzabo showed that κ a can be related to macroscopic rates by identifying κ a = k D k a /k D,b and K eq = κ a /k d . [70] Agmon and Szabo derived the same relation for the equilibrium constant( K eq = κ a /k d ) by considering an isolated pair of reactive spheres using a back-reactionboundary condition. [36]Here we show how the microscopic reaction rate used to describe reactions in our sim-ulation approach can be related to the reaction scheme in Eq. 2 and to the equilibriumconstant (Eq. 8). Inspired by the work of Shoup and Szabo [70] we calculate the equilibriumconstant for the formation of the encounter. In our approach the encounter is defined as aregion in configuration space (Eq. 3 and Eq. 4) around the desired relative position in thebound complex instead of a boundary as commonly used in the Fokker-Planck picture. Theconfiguration space of two rigid unbound particles in a periodic box of volume V can bedescribed by their relative position vector (cid:126)r , their center-of-mass vector (cid:126)R and 2 × (cid:126)ω and (cid:126)ω . Assuming free diffusive motion without reactions all configurationsin the two particle configuration space are equally probable and they only interact by stericrepulsion. By using the thermodynamic extremum principle for the Gibbs free energy, wecan relate the equilibrium constant for the formation of an encounter to the partition sumsof the free molecules z A and z B and to the partition sum of the encounter complex z A · B K enceq = k D k D,b = ( z A · B /V )( z A /V )( z B /V ) (9) z A = z B = (cid:90) V d x (cid:90) d ω = V × π × π (10) z A · B = (cid:90) V d R (cid:90) d r (cid:89) i =1 (cid:0) (cid:90) d ω i ) (11) ⇒ K enceq = (cid:90) d r (cid:89) i =1 (cid:0) π (cid:90) d ω i (cid:1) =: V (cid:63) . (12)For the single molecules all positions and orientations (neglecting the excluded volume) areaccessible, resulting in a factor of V from the integration over the translational degrees offreedom and a factor of 8 π from the orientational degrees of freedom of the rigid molecules.To determine the partition sum of the encounter complex z A · B , we first integrate out thecenter-of-mass coordinate (cid:126)R of the two molecules resulting in a factor of V . The boundariesof the remaining integrals over the relative position coordinate (cid:126)r and the orientation coor-dinates (cid:126)ω and (cid:126)ω are determined by the patch definition and the particle geometries. Theyfollow from Eq. 3 and Eq. 4 and have to be calculated for every specific case, as will bediscussed in more detail below. The resulting reactive volume V (cid:63) is the central concept torelate microscopic and macroscopic rates. From Eq. 12 we see that it equals the equilibriumconstant K enceq .If we now allow for reactions, not every configuration in the encounter region is equallypopulated as some particles will already react before they explore the inner region of V (cid:63) .However, if on average the encounter region is well explored before an association, we canapproximate the real distribution within V (cid:63) by a uniform distribution. This approximationis equivalent to the underlying assumption in the work of Pogson et al [77] and Klann et al.[78] who introduced a reactive volume (and an intrinsic association rate [78]) based on theassumption that the particles are randomly distributed in the box at all times.Alsallaq and Zhou also used thermodynamic arguments to determine the equilibriumconstant for the bound complex. [44] However, in contrast to this work we consider allconfigurations within the generalized reactive volume V (cid:63) as encounter configurations. Byintroducing finite reactivity with the microscopic reaction rate k (cid:48) a with which a bound com-plex can be formed from the encounter region (see Fig. 3a), our approach can account forreactions which are reaction- or diffusion-limited, as will be shown below in the results4section. To study assembly, this is of great importance as intermediates emerging duringthe assembly process can have very different diffusion properties and a reaction that wasreaction-limited for small clusters can become diffusion-limited for larger clusters.Thus, by identifying the equilibrium constant for the formation of the encounter withthe reactive volume V (cid:63) (Eq. 12) and relegating bond association and dissociation to thereaction rates, we can identify the reaction rates k (cid:48) a and k (cid:48) d in the microscopic frameworkwith the corresponding reaction rates k a and k d in the macroscopic framework (Eq. 8) andthe equilibrium constant is given by: K eq = k D k D,b k a k d = V (cid:63) k a /k d . (13)In accordance with the work by Shoup and Szabo [70] the rate κ a = V (cid:63) k a is given by theproduct of the equilibrium constant for the formation of the encounter complex and thereaction rate for the formation of the final complex. This rate can be identified with thereaction rate used by Morelli and ten Wolde [50] ( k a in ref. [50]).While the reactive volume V (cid:63) = k D /k D,b is sufficient to predict the correct equilibriumconstant (Eq. 13), the kinetics of the reactions depend on the absolute values of the diffusiverates k D and k D,b . Thus in order to determine the macroscopic rates k + and k − defined inEq. 7, we need to evaluate k D . This can be done within our simulation approach based on analgorithm by Zhou [39]. In this algorithm k D can be calculated from the survival probability S ( t ) of two particles starting in encounter by [39, 43]: k D = lim t →∞ κ a S ( t )1 − S ( t ) . (14)Given the diffusive on-rate k D , the diffusive off-rate k D,b can be simply calculated by k D,b = k D /V (cid:63) (Eq. 9).Any simulation algorithm with underlying reversible dynamics needs to satisfy detailedbalance. To establish detailed balance we follow the approach by Morelli and ten Wolde [50]who inferred from a detailed balance consideration that the relative position distributionof 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 spherescovered with a reactive shell, it has been argued that detailed balance can be introduced byplacing the two hard spheres according to a radial uniform distribution within the reactiveshell followed by one additional move step [79]. Here we generalize this idea as follows:5Upon the dissociation of two clusters, a new configuration is chosen uniformly within theencounter region V (cid:63) followed by a diffusive move step of the two clusters. To establishthe uniform distribution in V (cid:63) we exploit the time invariance symmetry of the diffusionpropagator. This symmetry ensures that in a confined volume of the configuration spacea uniform distribution is established as the steady state distribution. As V (cid:63) describes aconfiguration around the desired relative configuration of the proteins, we can generate auniform distribution without prior knowledge of the exact shape of the reactive volume byperforming various “pseudo-diffusion” steps of the two clusters starting from their predefinedrelative configuration. These “pseudo-diffusion” steps are confined to the region V (cid:63) in therelative accessible configuration space of the two partners. We call this procedure “pseudo-diffusion” as the two dissociating clusters are propagated diffusively within V (cid:63) , however, thetimestep used to establish the uniform distribution has no physical meaning. This procedurefor generating a uniform distribution within the reactive area does not depend on the particleor patch geometry. Thus it is especially suited to study assembly dynamics as parts of thereactive volume might become sterically blocked during the assembly process. These stericeffects are automatically taken into account with our method. The effect of steric blockingand its consequences for assembly will be discussed in the results section. Moreover, thismethod would also remain valid when including long range electrostatic interactions. In thiscase 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 V (cid:63) .In general the encounter volume V (cid:63) has to be understood as the volume of all two particleconfigurations resulting in an encounter between two clusters (see Fig. 3b). This means thatfor clusters with multiple patches the total encounter volume V (cid:63) = (cid:83) V (cid:63)p i ,p j where V (cid:63)p i ,p j isthe encounter volume between a pair of patches. In the case of a dissociation of clusterswith multiple patches a uniform configuration in the subvolume V (cid:63)p i ,p j of the patches thatformed 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 complexformation is given by the sum of the different equilibrium constants for different bonds.6 F. Intrinsic association
Proteins with multiple reactive patches can assemble in structures containing loops. Inthis case an already connected cluster can contain open bonds with two patches being inpermanent encounter due to the geometry of the cluster. These patches can bind with theinternal association rate k intra a , where we again assume Poisson-distributed waiting times forthe bond formation. The bond formation in this case is fundamentally different from theabove discussed process as no diffusion is involved. Since clusters are rigid objects, internalbond formation does not change the shape of the cluster. Thus it can be considered as aninternal process stabilizing an already existing structure. Given a dissociation rate of k d for a specific bond, we can relate the internal association rate to the free energy E of bondformation. For a closed loop containing n bonds, there are n different configurations withone open bond in the loop and the fraction of the open and closed loop in equilibrium isgiven by: p open p closed = nk d k intra a = Z open Z closed = ne − βE ⇒ k intra a = k d e − β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 enhancedinternal movement. This will be reflected in an increase in the entropy which would shift theratio of open and closed ring structures towards the open state (smaller k intra a ). In principleit should be possible to calculate the change in energy and entropy upon dissociation of abond using detailed MD simulations. Such simulations could be used to specify the internalassociation rate. G. Outline of the algorithm
In this part we will describe the work flow of our algorithm. Once a starting configurationis seeded the following steps are repeated iteratively. Firstly every particle is translated androtated according to its diffusive properties in the framework of free Brownian motion withina timestep ∆ t . With the new positions and orientations all clusters with steric overlaps aswell as all clusters participating in a reaction are tagged. A reaction between two clustersoccurs with the bond specific probability k a ∆ t if they are in encounter. In a next step allsteric overlaps (including box overlaps if the box is non-periodic) are corrected by setting7the involved particles back to their original position before the move step. This might leadto 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 andorientations (see Fig. 1c). If a reaction leads to a steric overlap either between the two merg-ing clusters or with other clusters, the reaction is rejected and the particles reassume theirconfiguration before the move step. Collisions with other clusters after the reorientationsstep is a higher order effect in the concentration of the particles. In dilute systems, especiallyfor reactive configuration volumes around the desired relative configuration, this happensrarely. In a next step open, internal bonds can be closed with probability P accintra = k intra a ∆ t .Finally each existing bond can dissociate with the bond specific probability k d ∆ t . If thisleads to two unconnected cluster fragments, the diffusive properties of the fragments arecalculated. Subsequently a uniform distribution within the encounter volume V (cid:63) is realizedby “pseudo-diffusive” motion preserving the encounter followed by one unconstrained dif-fusion step of the two clusters. All simulations have been performed at room temperature T = 293 K and with the viscosity of water η = 1mPas. III. RESULTS AND DISCUSSIONA. The reactive volume
As discussed in the previous section, we can relate the microscopic rates to the macro-scopic equilibrium constant using the concept of a generalized reactive volume V (cid:63) in con-figuration space (Eq. 13). Here we introduce a model system in which it is possible toanalytically calculate V (cid:63) and show that the results from our MC scheme (Eq. 17) agree wellwith the analytical results. In this model system a protein is described by one hard sphereof radius R equipped with one polar patch as depicted in Fig. 3a. The patch is centered atthe origin of the protein. It has a radius of R p and an opening angle of θ p ∈ [0 , π ] aroundthe z-axis of the proteins. With our definition of the encounter complex for two such par-ticles (Eq. 3 and Eq. 4), we can analytically calculate V (cid:63) by decomposing the radial and8 E n c oun t e r v o l u m e V * ( n m ) Patch angle (cid:101) p MCAnalyt. calc
FIG. 4. V (cid:63) for two identical hard spheres of radius R = 1nm equipped with one polar patch ofradius R p = 1 . θ p . The analytical result from Eq. 16 is shownwith a dashed line while the MC estimates (according to Eq. 17) are represented by red pointswith error bars. orientational part, similar to the potential model proposed by Kern and Frenkel[80]: V (cid:63) = V (cid:63) rad × V (cid:63) ori (16) V (cid:63) rad = 43 π (( R p ) − R ), V (cid:63) ori = (cid:0) π (cid:90) π dφ (cid:90) θ p sin( θ ) dθ (cid:1) = 14 (1 − cos( θ p )) .The square in the orientation part arises from the fact that we have to consider two indepen-dently rotating particles and only if the orientation of both particles is correct, an encounteris reached (see Eq. 4). For θ p = π the orientation constraint vanishes and the sphericalresult 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 inte-gration scheme. In this scheme we sample different configurations by randomly positioningtwo clusters with random orientations in a periodic box of volume V box . We repeat this pro-cedure N times and can estimate V (cid:63) by counting the fraction of trials n = N encounter /N total that result in an encounter configuration. Then we simply have V (cid:63)MC = nV box . (17)In Fig. 4 the configuration space volume V (cid:63) for two identical, spherical proteins equippedwith a polar patch ( R = 1nm and R p = 1 . θ p . Comparing the analytical result (black line) and MC simulation (Eq. 17) we seethat our simple MC scheme correctly estimates V (cid:63) . 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 V (cid:63)p i ,p j of two specific patches p i and p j and the total encounter volume V (cid:63) of two proteins. The total configuration volume V (cid:63) is given by the union of all volumes V (cid:63)p i ,p j for the encounter of two different patches i , j that can bind to each other. If the encounter volumes of different patch combinations aredisjoint, the total configuration space volume for the encounter of two clusters is given bysumming over all V (cid:63)p i ,p j . For multiple, identical patches that can bind with each other, thetotal reactive volume V (cid:63) is given by multiplying the configuration volume V (cid:63)p i ,p j of one patchcombination with a combinatorial prefactor accounting for the multiplicity of the differentpatch configurations if the corresponding volumes are not overlapping in configuration space. B. Bimolecular reactions
To verify our algorithm and to show the importance of detailed balance, we first analyzea system with a small number of proteins only. In detail we consider a situation in which oneparticle of type A and N B particles of type B are placed in a periodic box of volume V box .In this case we can relate the concentration c C of the complex C to the fraction of time t b in which particle A is bound to particle B ( p b = t b / ( t u + t b )), while the concentration c A ofparticle A can be related to the fraction of time t u in which A is unbound ( p b = t u / ( t u + t b )).The unbound A particle is surrounded by a concentration of c B = N B /V particles of type B . Here V = V box − V excl is the freely accessible volume. For proteins of identical radius( R A = R B = R ) the excluded volume can be approximated by V exc = 4 πN B (2 R ) /
3. Thuswe can relate the probability that A is bound to the equilibrium constant of the reactionby: [36, 50] K eq = c C c B c A = p bound N B V (cid:0) − p bound (cid:1) ⇒ p bound ( N B ) = K eq N B K eq N B + V . (18)Eq. 18 follows from elementary thermodynamic arguments and is valid irrespective of thedetails of the binding interactions. In the following subsections we will show that the correctequilibrium bound probability is reproduced with our algorithm for spherically reactiveand patchy particles when appropriately taking the generalized encounter volume V (cid:63) into0 P bound N B a) Mean field theory (cid:54) t=0.01ns (balance) (cid:54) t=0.001ns (balance) (cid:54) t=0.0001ns (balance) (cid:54) t=0.01ns (contact) (cid:54) t=0.001ns (contact) (cid:54) t=0.0001ns (contact) 0 0.02 0.04 0.06 0.08 0.1 2 2.1 2.2 2.3 2.4 2.5 P (r) / ( (cid:47) r ) Radius (nm)b) r in , k a =10ns -1 r out , k a =10ns -1 r in , k a =0.1ns -1 r out , k a =0.1ns -1 FIG. 5. Equilibrium bound probabilities and radial distribution for spherically reactive particles.a) The probability of one A particle being bound to a B particle is shown as a function of the numberof B particles according to Eq. 18 (black line) for K eq = V box . The following parameters have beenused for the simulations: R A = R B = 1nm, R p,A = R p,B = 1 . θ p,A = θ p,B = π leading to V (cid:63) ≈ = 11 . (Eq. 16). For a given k a = 10ns − , the dissociation rate k d is chosen so that V (cid:63) k a k d = V box = 8000nm (Eq. 13). Circles depict the simulation results obtained with our algorithmfor different time resolutions while stars correspond to simulations in which detailed balance isviolated by placing particles in contact after dissociation. b) The normalized radial distribution∆ t before reaction (circles) and after dissociation including the additional unconstrained move step(stars) is shown for two different association rates k a = 10ns − (red) and k a = 0 . − (green) witha time resolution of ∆ t = 0 . V (cid:63) (red points) and the expected uniformdistribution in a spherical shell (black line). p bound from Eq. 18.
1. Spherically reactive particles
Here we compare our BD algorithm with the mean field result (Eq. 18) for the case ofspherically reactive A and B particles ( R A = R B = 1nm, R p,A = R p,B = 1 . θ p,A = θ p,B = π , V (cid:63) = 11 . ) and show that our algorithm fulfills detailed balance in thiscase. Given an association rate of k a = 10ns − we chose k d so that K eq = V (cid:63) k a k d = V . InFig. 5a the probability of particle A being bound to particle B is shown as a function ofthe number of B particles N B and compared to the mean field result (black line). Placingparticles uniformly in the reactive shell by “pseudo-diffusive” motion and allowing for oneunconstrained move step of the two proteins involved in the dissociation according to thealgorithm described above (circles) leads to very good agreement between the mean fieldtheory and the simulation results. Only for a very large timestep ∆ t = 0 . . p bound and the simulation results cannot be reconciled with the mean field theory. Althoughreducing the timestep from ∆ t = 0 . t = 0 . k a we observe the same radialdistribution before reaction and after dissociation and thus our algorithm satisfies detailedbalance for this quasi one-dimensional case. In the inset of Fig. 5b we show a histogram of thedistance between two dissociating particles after the “pseudo-diffusive” motion constraint to V (cid:63) without the additional, unconstrained move step of the two particles. Here we see thatwe are indeed able to generate a uniform distribution within V (cid:63) by exploiting the symmetry2 P bound N B Mean field theory (cid:54) t=0.01ns (balance) (cid:54) t=0.001ns (balance) (cid:54) t=0.0001ns (balance) (cid:54) t=0.01ns (contact) (cid:54) t=0.001ns (contact) (cid:54) t=0.0001ns (contact)
FIG. 6. Equilibrium bound probabilities for patchy particles. The probability of one A particlebeing bound to a B particle is shown as a function of the number of B particles according toEq. 18 (black line) for K eq = V box . The mean field theory is independent of the specific particledetails. The following parameters have been used for the simulations: R A = R B = 1nm, R p,A = R p,B = 1 . θ p,A = θ p,B = π/ V (cid:63) ≈ = 0 . (Eq. 16). For a given k a = 10ns − the dissociation rate k d is chosen so that V (cid:63) k a k d = V box = 8000nm . Circles show the results of ouralgorithm for various timesteps while stars correspond to simulations in which detailed balance isviolated by placing particles in contact and with perfect alignment of the patches after dissociation. of the diffusion propagator as the histogram of the simulated positions agrees very well withthe expected uniform distribution.
2. Patchy Particles
After verifying that our algorithm reproduces the expected results for spherically reactiveparticles, we will now investigate the effect of anisotropic reactivity and show that ourgeneralized definition of the reactive volume enables us to compare our simulation resultswith the mean field theory, which is the same as for spherically reactive particles. Here weconsider the same setup as in the previous part with one particle of type A and N B particlesof type B in a periodic box ( R A = R B = 1nm and R p,A = R p,B = 1 . θ p = π/
4) allowing only for anassociation of two particles if the patches are sufficiently aligned with the ctc-vector of the3 (cid:113) a) before associationk a =0.1 0 0.01 0.02 0.03 0.04 0.05 b) after dissociationk a =0.1 2.05 2.1 2.15 2.2 2.25 2.3r (nm) 1.6 1.8 2 2.2 2.4 2.6 2.8 3 (cid:113) c) before associationk a =10. 0.04 0.035 0.03 0.025 0.02 0.015 0.01 0.005 2.05 2.1 2.15 2.2 2.25 2.3r (nm) 0 0.01 0.02 0.03 0.04 0.05 d) after dissociationk a =10. FIG. 7. Color-coded histogram of the particle positions before reaction and after dissociationas a function of the relative distance r and the angle between the orientations of the two patches φ . This histogram has been weighted by the generalized bin volume (4 πr dr × π sin( φ ) dφ ) .Several contour lines are shown for better comparability. In a) and c) the distribution of particlepositions before a reaction is shown for k a = 0 . − and k a = 10ns − , respectively. In b) and d)the corresponding distributions of particle positions after dissociation is shown. The distributionshave been recorded with a timestep resolution of ∆ t = 0 . proteins (Eq. 4). In this case the reactive volume reduces to V (cid:63) = 0 . according toEq. 16. Given the reactive volume V (cid:63) and the association rate k a = 10ns − , we again choose k d in such a way that K eq = V . The comparison between the simulation and mean fieldtheory is shown in Fig. 6. By using the appropriate scaling with the generalized volume V (cid:63) 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 proposeddetailed balance algorithm (circles). When placing the particles in contact (alignment ofpatch orientation with ctc-vector and radial contact), we drastically overestimate the boundprobability (stars). Compared to the spherical case (Fig. 5) this overestimation of p bound ismuch stronger. Thus, in the case of localized reactivity it is even more important to carefullyconsider the relative position and orientation of dissociating proteins. Moreover our resultsconfirm the relation between microscopic rates and the macroscopic equilibrium constant4given in Eq. 13.In Fig. 7 we again show the distribution of the relative position and orientation of twoparticles equipped with a polar patch before a reaction and after dissociation. Here φ corresponds to the angle between the orientation vectors of the patches and r is the distancebetween the centers of the two particles. Due to the symmetry of the particles aroundthe z-axis these two parameters provide a good description of the relevant configurationspace. The normalized probability distribution p ( r, φ ) / (4 πr × πsin ( φ )) is shown color-coded for two microscopic association rates k a = 0 . − (Fig. 7 upper part) and k a =10ns − (Fig. 7 lower part) at a time resolution of ∆ t = 0 . − . In the left part ofFig. 7 the distribution before reaction is shown while in the right part the correspondingdistribution after dissociation is depicted. As expected the distribution has its maximumfor particles in contact and anti-parallel aligned patch orientations ( r = 2 . φ = π ). Although the distributions fluctuate more for a smaller association constant (Fig. 7aand 7b), as can be seen by looking at the contour lines, very good agreement betweenthe distribution before association and after dissociation is observed for both associationparameters. 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 affectthe distributions suggesting that our algorithm works well for a large spectrum of microscopicparameters. This first example of non-spherically reactive particles shows that it is essentialto appropriately scale the macroscopic parameters using the concept of a generalized reactivevolume in configuration space. Moreover, careful consideration of the relative position andorientation of the dissociating particle is necessary to obtain quantitative agreement betweenBD simulations and mean field theory or experimental results. C. Beyond bimolecular reaction: Assembly of a pentameric ring
In this section we analyze the assembly of a pentameric ring structure consisting of 5proteins (Fig. 8b) as an example of an assembly process involving more than two partners.As elementary building blocks we consider spherical proteins of radius R = 1nm, eachequipped with two reactive patches (Fig. 8a). The patches are centered at the origin of theprotein and both have a radius of of R p = 1 . θ p = π/
5. Thegeometry of the ring is reflected in the different orientation vectors (cid:126)o associated with each5
FIG. 8. Simulation of the assembly of a pentameric ring structure. a) Elementary sphericalbuilding block of radius R = 1nm for the pentameric ring equipped with two patches. The anglebetween the centers of the patches is determined by the desired ring structure. Each patch hasan opening angle of θ p = π/ R p = 1 . patch, and the relative position and orientation of the proteins are encoded in the bondstructure.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 demon-strate how the macroscopic reaction rates can be calculated from the microscopic reactionparameters and diffusive properties of the intermediates. The parameter free comparison ofour simulation results with the rate equation predictions shows excellent agreement. Thisdemonstrates that we are indeed able to predict macroscopic reaction rates from our simu-lations taking into account changes in diffusive properties of the assembly intermediates aswell as steric effects arising during the assembly process.
1. Macroscopic rates
In Eq. 7 the reaction rates k + and k − are given in the case of a bimolecular reaction. Wenow specify the macroscopic rates k i,j + and k i,j − for the reaction between two ring fragments6 Size 1 2 3 41 0.40 0 .
40 0 .
40 0 .
242 0 .
40 0 .
40 0 .
36 -3 0 .
40 0 .
36 - -4 0 .
24 - - -TABLE I. Generalized reactive volumes V (cid:63)i,j in nm for the encounter of different ring fragmentscalculated by MC sampling according to Eq. 17 and analytically according to Eq. 16 for V (cid:63) , ,respectively. Size 1 2 3 41 3.29 2.69 2.15 1.172 2.69 2.03 1.52 -3 2.15 1.52 - -4 1.17 - - -TABLE II. Diffusive on-rates k D for different ring fragments in nm / ns. The rates have beencalculated from the survival probabilities of fragments starting in a encounter according to thealgorithm proposed by Zhou[39] with a minimum timestep resolution of ∆ t = 0 . with i and j proteins. These rates are given by: k i,j + = ˜ k i,j + (cid:122) (cid:125)(cid:124) (cid:123) k i,jD k a k i,jD,b + k a (1 − δ i,j ) (19) k i,j − = k i,jD,b kdk i,jD,b + k a (cid:124) (cid:123)(cid:122) (cid:125) ˜ k i,j − (2 − δ i,j ). (20) k i,j + describes the association of two ring fragments of size i and j . As the diffusive propertiesand the accessible reaction volume change with the size of the fragments k i,j + is different fordifferent combination of fragments, albeit the microscopic reaction rate k a remains constant.For the reaction of identical particles the macroscopic rate is reduced by a prefactor of 1 / k a . This reflects the fact that for reactions of identicalcomponents ( A + A → C ) the observed macroscopic rate is by a factor of 1 / N A ( N A − / k i,j − describes the macroscopic dissociation rate of a cluster of size i + j into two ring fragments of size i and j , respectively. To determine the macroscopicreaction rates given in Eq. 19 and Eq. 20 we need to calculate the diffusive on- and off-rates k i,jD and k i,jD,b = k i,jD /V (cid:63)i,j which change for different combination of fragments. Here V (cid:63)i,j is theconfiguration volume associated with the reaction of two fragments of size i and j .The reactive volume for two different fragment configurations i and j can be calculated byMC sampling as has been described in Eq. 17. The resulting volumes V (cid:63)i,j for two fragmentsof size i and j are shown in Table I. For the reaction of two monomers the reactive volume V (cid:63) , can also be evaluated by using Eq. 16 for the encounter of one specific pair of patches. Takinginto account the multiplicity of bond interactions 2 results in a total encounter volume of V (cid:63) = 0 . . While fragments not combining into a full ring ( i + j <
5) have the sameconfiguration volume as the monomeric fragments, the configuration volume for fragmentsthat can combine into a full ring is reduced. This reflects the fact that, independent of theparticle geometry, the encounter volume and thus the number of configurations in this volumeremains unchanged, unless parts of the reactive volume become sterically blocked, which isthe case for fragments that can assemble into the full ring structure. Moreover, for thesefragments some configurations with a simultaneous overlap of different patch combinationsexist. 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 arenaturally included in our simulation algorithm and for the simulations only knowledge of thereduced volume for the basic building blocks V (cid:63) , is required. However, in order to compareour simulation results to a macroscopic rate equation approach, we need to determine thechanges in the encounter volumes for different assembly intermediates.Moreover, in order to capture the reaction kinetics with the above defined rates, we needto determine the diffusive on-rates k i,jD . In order to calculate the diffusive on-rate k i,jD wefollow a scheme which was originally proposed by Zhou. [39]. Here, two ring fragmentsof size i and j are initially placed in a random configuration in encounter (e.g. a randomconfiguration within V (cid:63)i,j is chosen). Starting from this configuration the two fragmentsare propagated diffusively within our simulation framework. If the two fragments form anencounter, they can react with the probability k a ∆ t . In this case the run is terminated. Byrecording the survival probability S i,j ( t ) (the fraction of runs which survive until time t ) k i,jD κ i,ja = k a V (cid:63)i,j . In order to calculate S i,j ( t → ∞ )we use a microscopic reaction rate of k a = 1ns − and an adaptive timestep scheme with aminimum timestep of ∆ t min = 0 . i and j is given by k i,jD,b = k i,jD /V (cid:63)i,j .The absolute values for the different ring fragments are given in Tab. II. As fragments with( i + j ) > k i,jD = 0.
2. Macroscopic rate equation approach
After having defined the macroscopic reaction rates for the association and dissociationof different ring fragments (Eq. 19 and Eq. 20) we will now introduce a system of ODEs todescribe the changes in concentration c i of a ring fragment of size i . This description assumesthat the system is homogeneous and well mixed at all times. Stochastic fluctuations arisingfrom the finite number of proteins are neglected. While for ring fragments containing lessthan 5 proteins only one state exists, two states exist for the full ring: an open ring with4 bonds and a closed ring with 5 bonds. The closed ring can be considered as an internalstate as the last bond formation in our framework does not involve any diffusion and itis connected to the open state by the internal association rate k intra a (see Eq. 15). Theconcentration of the open state is denoted by c while the concentration of the closed stateis denoted by c (cid:63) . The changes in concentration for a ring consisting of N proteins can bedescribed by the following set of ODEs:9˙ c i = − N − i (cid:88) l =1 k i,l + (1 + δ i,l ) c i c l − i − (cid:88) l =1 l ≤ i − l k l,i − l − c i − k intra a δ i,N c N + i − (cid:88) l =1 l ≤ i − l k l,i − l + c l c i − l + N − i (cid:88) l =1 k i,l − (1 + δ i,l ) c l + i + N δ i,N k d c (cid:63)N (21)= − N − i (cid:88) l =1 ˜ k i,l + c i c l − i − (cid:88) l =1 ˜ k l,i − l − c i − k intra a δ i,N c N + 12 i − (cid:88) l =1 ˜ k l,i − l + c l c i − l + 2 N − i (cid:88) l =1 ˜ k i,l − c l + i + N δ i,N k d c (cid:63)N (22)˙ c (cid:63)N = − N k d c (cid:63)N + k intra a c N . (23)In our case N = 5. Eq. 21 shows the different processes which lead to a change in theconcentration of a fragment containing i proteins. The first three terms lead to a decrease inconcentration. The first term describes the change in concentration due to the reaction of afragment of size i with another fragment of size l . If i = l a twofold decrease in concentration c i is observed. The second term describes the change in concentration due to the decay ofa fragment of size i into two fragments of size l and i − l . To prevent the double countingof dissociation events the constraint l ≤ i − l for the summation is used here. In the thirdterm the internal bond formation from the open to the closed state is described. Similarlyin the last three terms all contributions which lead to an increase in concentration c i due toassociation and dissociation processes are described. In Eq. 23 the dynamics of the closedpentameric ring state is separately described. In the derivation of Eq. 22 from Eq. 21 thesymmetry of the matrices k i,l ± = k l,i ± has been exploited. The above described set of ODEsobeys mass conservation ( (cid:80) i ( c i + δ i, c (cid:63) ) i = const). Thus, one of the concentrations couldbe 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 fragmentconcentration in Eq. 22, we will now compare our BD simulation results to the macroscopicrate equation approach. As all macroscopic rates are fully specified by the microscopicparameters (Eq. 19 and Eq. 20), the comparison between the simulation and macroscopic0 a) N o r m a li z ed c on c en t r a t i on k a =1 ns -1 , V box =(55 nm) k d =0.0001 ns -1 k aintra =0.0 ns -1
0 5000 10000 15000 b) k a =8 ns -1 , V box =(110 nm) k d =0.0001 ns -1 k aintra =0.0 ns -1 c) N o r m a li z ed c on c en t r a t i on k d =0.0001 ns -1 k aintra =0.001 ns -1
0 10000 20000 30000 40000 d) k d =0.0001 ns -1 k aintra =0.001 ns -1 e) N o r m a li z ed c on c en t r a t i on Time (ns)k d =0.0 ns -1 k aintra =0.0 ns -1 f) Time (ns)k d =0.0 ns -1 k aintra =0.0 ns -1 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 = ic i c with c = N/V are shown for different parameters. All simulations were started with N = 500 individualproteins with a time resolution of ∆ t = 0 . k a = 1ns − and V box = (55nm) . The right column corresponds to diffusion-influenced assemblywith k a = 8ns − and V box = (110nm) in which the increase in k a is balanced by the decreasein concentration. In a) and b) the result for a reversible reaction with a dissociation rate of k d = 0 . − and without internal bond formation is ( k intra a = 0 . ns − ) shown. In c) and d) thefinal ring is additionally stabilized by the formation of an internal bond with rate k intra a = 0 . − (see Eq. 15). In e) and f) the evolution of the cluster size distribution is shown for irreversiblebinding ( k d = 0ns − ). theory does not involve any free parameter. For the BD simulations we randomly position N = 500 single proteins in a periodic box of size V box and record the time-course of thefragment concentrations. All BD simulations were performed at a timestep resolution of∆ t = 0 . c = c ( t = 0) = N/V . Here V is the accessible simulation volume which is roughlyestimated from the box volume by V = V box − ( N − π (2 R ) / c i = ( c i + δ i, c (cid:63) ) i/c is shown. The normalized concentration of a fragment is defined as the concentration c i weighted by the number of proteins contained in the fragment and normalized by the initialconcentration. The normalized concentration predicted by the rate equation approach isshown with solid lines, while the averaged simulation results are represented by dashed lines.In the left column of Fig. 9 (Fig. 9 a,c and e) a microscopic reaction rate of k a = 1ns − and a box volume of V box = (55nm) are chosen. As the diffusive off-rates range from k , D,b ≈ . − to k , D,b ≈ . − , this assembly process can be considered as reaction-limited( k a < k D,b ). In the case of a reaction-limited process, the macroscopic off- and on-rates givenin Eq. 19 and Eq. 20 simplify to ˜ k i,j + ≈ k i,jD /k i,jD,b k a = V (cid:63)i,j k a and ˜ k i,j − ≈ k d . Thus, in this casethe rates only depend on the microscopic reaction rates and of the encounter volume V (cid:63)i,j .In the right column of Fig. 9 (Fig. 9 b,d and f) a microscopic association rate of k a = 8ns − is chosen. This 8-fold increase in the microscopic association rate is compensated by an8-fold decrease in concentration. In this case the reaction is strongly affected by diffusion( k a ≥ k D,b ) and the macroscopic rates crucially depend on the different diffusive on- and off-rates k i,jD and k i,jD,b , which vary with the diffusive properties of the different fragments as canbe seen in Tab. II and a clearly different assembly dynamic is expected. We will refer tothis case as diffusion-influenced assembly.In Fig. 9a and Fig. 9b the normalized fragment concentrations for the assembly withreversible bonds ( k d = 0 . − ) but without an additional stabilization of the full ringstructure ( k intra a = 0) is shown. Comparing the predictions from the rate equation approach(Eq. 22) (solid lines) with the averaged simulation results (dashed lines) we observe excellentagreement between our simulation results and the results from the macroscopic rate equationapproach 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 inFig. 9a can be explained by our rough estimate for the excluded volume. This agreementbetween the two approaches, especially for the diffusion influenced regime is remarkable anddemonstrates that, by evaluating the relevant diffusive properties of the fragments, we areindeed able to relate our microscopic reaction parameters to macroscopic reaction rates that2correctly reproduce the dynamics of the system. Comparing the equilibrium distributionof the cluster fragments in Fig. 9a and Fig. 9b we see that indeed the 8-fold decrease inconcentration is balanced by the 8-fold increase in concentration, as it was expected dueto our previous reasoning. However, comparing the kinetics of assembly we see a cleardifference between Fig. 9a and Fig. 9b. In general, the approach towards equilibrium isslower in Fig. 9b than in Fig. 9a. This shows that we are indeed in a different assemblyregime and the change in concentration is not compensated by the change in k a in the kineticsof the assembly process. At lower concentration (Fig. 9b) clusters need a longer time to findeach other diffusively which is reflected in the slower approach towards equilibrium.In Fig. 9c and Fig. 9d the normalized fragment concentrations with an additional sta-bilization of the final ring structure k intra a = 10 k d = 0 . − is shown. As discussedpreviously, the internal bond formation can be understood as an intrinsic process which sta-bilizes the full ring (Eq. 23). The value chosen for the internal association rate correspondsto a free energy of E ≈ − . k B T associated with the bond formation according to Eq. 15.In Fig. 9c and Fig. 9d, the normalized fragment concentration of the full ring (black line) isshown irrespective of the number of bonds in the ring. Comparing the simulation results andthe rate equation approach we again observe very good agreement between them. Similarlyto the case without internal bond formation (Fig. 9a and Fig. 9b), the assembly dynamicsat lower concentration (Fig. 9d) is slower than at higher concentration (Fig. 9c), while theequilibrium steady state remains the same in both cases. In contrast to the reversible bondformation without additional stabilization of the ring (Fig. 9c and Fig. 9e) the normalizedconcentration of the full ring is significantly increased by the internal bond formation. Dueto mass conservation the increase in the concentration of the full ring structure is balancedby a decrease in the concentration of smaller fragments showing how an additional internalstate can alter the equilibrium cluster size distribution.Finally, we investigate the assembly dynamics for irreversible bond formation ( k d =0ns − ) for two different microscopic reaction rates k a corresponding to the regime of reaction-limited reactions (Fig. 9e) and diffusion-influenced reactions (Fig. 9f). Again the simulationresults and the rate equation approach agree remarkably well with each other. By comparingFig. 9e and Fig. 9f we again verify that the approach towards the steady state is muchslower in the case of lower concentration shown in Fig. 9f. However, in contrast to thepreviously discussed reversible bond formation, the steady state is different in both cases3with different fraction of ring fragments containing 3,4 and 5 proteins. This reflects afundamental difference between the steady state observed for reversible bond formation andthe steady state observed for irreversible bond formation. In the case of reversible bondformation the steady state is an equilibrium state in which dissociation and associationevents balance each other. This state only depends on the microscopic reaction rates k a and k D and the size of the reactive volume V (cid:63) (Eq. 13). In contrast the steady state forirreversible bond formation is a trapped steady state without any underlying dynamics. Inthis case of irreversible bond formation the ring fragments cannot disassemble and fragmentsof size i ≥ / IV. CONCLUSION
Biological structures like the actin cytoskeleton [13–15] or the nuclear pore complex [7–9]are highly dynamic with their constituents being continuously exchanged. To understandthe biological functionality of these fascinating systems, a detailed understanding of the as-sembly dynamics is crucial. Moreover, advances in the fabrication of colloidal particles withdirected interactions (“patchy particles”) allow to design a plethora of differently shapedbuilding 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 pre-vented. Recently it has been shown that particle interactions during the assembly process4can be actively controlled. [24, 25] Thus, by a state- or time-dependent switching of reac-tivity the assembly process can be actively steered to reduce kinetic trapping. However, tofully exploit the possibilities of controlling the assembly process, a detailed understandingof its full dynamics with the formation of assembly intermediates is necessary.Here we presented a novel simulation approach which is ideally suited to investigate thedynamics of large protein assemblies with well-defined architectural properties. Proteinsare 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 relyon force fields to describe protein interactions, we combine overdamped Brownian motionwith the concept of reversible, stochastic reactivity for patchy particles. Each cluster istreated as rigid object with its diffusive properties being evaluated on-the-fly. If the reactivepatches of two clusters form an encounter by force-free diffusive motion, a bond betweenthe clusters is established with a predefined rate. Local rules are used to describe theresulting rearrangements, which are assumed to be very fast on the time scale of the BDsimulations. For potential-based simulations, these rearrangements would proceed due tothe corresponding forces, which do not exist in our approach. Similar to the associationprocess, the dissociation process also occurs with a predefined rate. Here we have derivedthe 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 duringthe assembly are diffusing as rigid objects (thus reducing the number of degrees of freedom tobe propagated to six) and because interactions are treated locally (as opposed to evaluatingpotentials). Nevertheless, the patchy particle approach includes detailed information on thearchitecture of protein assemblies, thus allowing us to study the spatial-temporal dynamicsof specific biological systems of interest. Moreover, we have been able to show how the mi-croscopic reaction parameters used in our approach correspond to the macroscopic reactionrates commonly used to describe the dynamics of the concentration of assembly interme-diates. This allows a direct comparison between macroscopic, experimentally measurableconcentrations and simulation results. To the best of our knowledge, this is the first timethat a BD algorithm combining stochastic reactivity of anisotropic patches with reversibledynamics that fulfills detailed balance is presented. Testing our algorithm against meanfield prediction for the case of bimolecular reactions revealed the importance of satisfyingdetailed balance.5As an example for a multi-component cluster, we investigated the assembly of a pen-tameric ring. Here we compared our simulation results to a rate equation approach forthe different fragment concentrations. By extracting the diffusive on- and off-rates fromour simulations we find excellent agreement between the rate equation approach and oursimulation results even for the case of strongly diffusion-influenced assembly without anyfree parameter. We showed that the assembly dynamics can be significantly changed by thechange in the diffusive properties of the assembly intermediates. In the case of irreversiblebond formation the system becomes trapped with incompatible ring fragments. In this casenot only the assembly kinetics but also the finally reached steady state depends on the dif-ferent diffusive properties of the assembly intermediates. For complex assembly geometriestrapping is a common motif. Thus, in this case not only correctly predicted equilibriumstructure distribution of the assembly intermediates is relevant but also the correctly pre-dicted kinetics of the assembly process as the system might become trapped before reachingthe expected equilibrium steady state. By demonstrating how the macroscopic reaction ratescan be calculated from our simulation, we have developed a framework in which a qualita-tive prediction of the changes in the macroscopic reaction rates, based on the changes of thediffusive properties of the assembly intermediates, is possible.Because our approach aims to describe the assembly of large protein structures, the detailsof individual bond formations are not resolved within this framework. In general, the use oflocal rules requires well-defined target geometries and the instantaneous local rearrangementaccompanying stochastic bond formation follows from the assumption that the time scalesof assembly and local rearrangements are well separated. Given the coarse-grained natureof our approach, the instantaneous rearrangement during binding is a valid approximationif the encounter volume specifies a narrow region around the desired configuration. In thecase of large encounter volumes the assumption of fast and small local rearrangements canbreak down and the use of local rules needs to be considered with care. In particular inthe case of higher protein densities or assembly close to a membrane this can lead to asituation in which physically reasonable rearrangements are suppressed by our steric rules.In classical MD or BD simulations with potentials, these rearrangements would be realizeddue to ensuing mechanical forces. Typical biological examples for situations which are outof the scope of our current approach are the bending of a membrane by BAR-proteins[81, 82], the emergence of polymorphic virus structures [83–85] or the assembly of actin6networks from preexisting large fibers [86]. In principle, such a situation could be resolved byusing hybrid schemes that interface our approach with potential-based simulations for localrearrangements. For the time being, however, we focus on the assembly of protein complexeswith well-defined architectures, like native viruses[87], centrioles [6, 88–90] or actin filamentsin 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 proteinclusters. After we have successfully verified our approach by comparing averaged simulationdata with macroscopic mean field results, we now can investigate the stochastic varianceinherent to these self-assembly processes. Our approach is ideally suited to study the effectof time- or state-dependent changes in reactivity, as has been recently shown in a qualitativestudy on the effect on hierarchical assembly in virus capsids. [26] Moreover, this approachmight also help to design artificial systems that do not get kinetically trapped in undesiredstructures. It can also be coupled with hydrodynamic schemes, in particular with reactivemulti-particle collision dynamics [92].
ACKNOWLEDGMENTS
HCRK acknowledges financial support by the Cusanuswerk. USS is a member of theHeidelberg cluster of excellence CellNetworks. We thank Nils Becker and Joris Paijmans forhelpful discussions.7 [1] R. R. Gabdoulline and R. C. Wade, Biophysical Journal , 1917 (1997).[2] W. Roos, R. Bruinsma, and G. Wuite, Nature Physics , 733 (2010).[3] A. Zlotnick and S. Mukhopadhyay, Trends in Microbiology , 14 (2011).[4] D. A. Compton, Annual Review of Biochemistry , 95 (2000).[5] E. Karsenti, F. N´ed´elec, and T. Surrey, Nature Cell Biology , 1204 (2006).[6] P. G¨onczy, Nature Reviews Molecular Cell Biology , 425 (2012).[7] F. Alber, S. Dokudovskaya, L. M. Veenhoff, W. Zhang, J. Kipper, D. Devos, A. Suprapto,O. Karni-Schmidt, R. Williams, B. T. Chait, et al. , Nature , 683 (2007).[8] M. A. D’Angelo and M. W. Hetzer, Trends in Cell Biology , 456 (2008).[9] E. J. Tran and S. R. Wente, Cell , 1041 (2006).[10] B. Geiger and A. Bershadsky, Current Opinion in Cell Biology , 584 (2001).[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, European Journal of CellBiology , 165 (2006).[12] H. Wolfenson, Y. I. Henis, B. Geiger, and A. D. Bershadsky, Cell Motility and the Cytoskele-ton , 1017 (2009).[13] C. C. Beltzner and T. D. Pollard, Journal of Biological Chemistry , 7135 (2007).[14] K. Guo, J. Shillcock, and R. Lipowsky, The Journal of Chemical Physics , 015102 (2009).[15] T. D. Pollard, Annual Review of Biophysics and Biomolecular Structure , 451 (2007).[16] K. Guo, J. Shillcock, and R. Lipowsky, The Journal of Chemical Physics , 155105 (2010).[17] G. Schreiber, G. Haran, and H.-X. Zhou, Chemical Reviews , 839 (2009).[18] C. A. Mirkin, R. L. Letsinger, R. C. Mucic, and J. J. Storhoff, Nature , 607 (1996).[19] L. Di Michele and E. Eiser, Physical Chemistry Chemical Physics , 3115 (2013).[20] Y. Wang, Y. Wang, D. R. Breed, V. N. Manoharan, L. Feng, A. D. Hollingsworth, M. Weck,and D. J. Pine, Nature , 51 (2012).[21] C. Knorowski and A. Travesset, Current Opinion in Solid State and Materials Science , 262(2011).[22] S. Sacanna, W. Irvine, P. M. Chaikin, and D. J. Pine, Nature , 575 (2010). [23] S. Sacanna, M. Korpics, K. Rodriguez, L. Col´on-Mel´endez, S.-H. Kim, D. J. Pine, and G.-R.Yi, Nature Communications , 1688 (2013).[24] M. E. Leunissen, R. Dreyfus, F. C. Cheong, D. G. Grier, R. Sha, N. C. Seeman, and P. M.Chaikin, Nature Materials , 590 (2009).[25] L. Di Michele, F. Varrato, J. Kotar, S. H. Nathan, G. Foffi, and E. Eiser, Nature Communi-cations , 2007 (2013).[26] J. E. Baschek, H. C. Klein, and U. S. Schwarz, BMC Biophysics , 22 (2012).[27] M. v. Smoluchowski, Zeitschrift f¨ur physikalische Chemie , 9 (1917).[28] F. C. Collins and G. E. Kimball, Journal of Colloid Science , 425 (1949).[29] K. ˇSolc and W. Stockmayer, International Journal of Chemical Kinetics , 733 (1973).[30] D. Shoup, G. Lipari, and A. Szabo, Biophysical Journal , 697 (1981).[31] O. Berg, Biophysical Journal , 1 (1985).[32] A. Shushin, The Journal of Chemical Physics , 12044 (1999).[33] M. Schlosshauer and D. Baker, The Journal of Physical Chemistry B , 12079 (2002).[34] M. Schlosshauer and D. Baker, Protein Science , 1660 (2004).[35] N. Agmon, The Journal of Chemical Physics , 2811 (1984).[36] N. Agmon and A. Szabo, The Journal of Chemical Physics , 5270 (1990).[37] S. H. Northrup, S. A. Allison, and J. A. McCammon, The Journal of Chemical Physics ,1517 (1984).[38] S. H. Northrup and H. P. Erickson, Proceedings of the National Academy of Sciences , 3338(1992).[39] H. X. Zhou, Journal of Physical Chemistry , 8794 (1990).[40] H.-X. Zhou, Biophysical Journal , 1711 (1993).[41] R. R. Gabdoulline and R. C. Wade, Journal of Molecular Biology , 1139 (2001).[42] G. Zou and R. D. Skeel, Biophysical Journal , 2147 (2003).[43] R. Alsallaq and H.-X. Zhou, Structure , 215 (2007).[44] R. Alsallaq and H.-X. Zhou, Biophysical Journal , 1486 (2007).[45] S. Qin, X. Pang, and H.-X. Zhou, Structure , 1744 (2011).[46] M. Klann and H. Koeppl, International Journal of Molecular Sciences , 7798 (2012).[47] S. S. Andrews and D. Bray, Physical Biology , 137 (2004).[48] S. S. Andrews, Physical Biology , 111 (2005). [49] J. S. van Zon and P. R. Ten Wolde, The Journal of Chemical Physics , 234910 (2005).[50] M. J. Morelli and P. R. Ten Wolde, The Journal of Chemical Physics , 054112 (2008).[51] K. Takahashi, S. T˘anase-Nicola, and P. R. ten Wolde, Proceedings of the National Academyof Sciences , 2473 (2010).[52] D. T. Gillespie, The Journal of Physical Chemistry , 2340 (1977).[53] J. Hattne, D. Fange, and J. Elf, Bioinformatics , 2923 (2005).[54] D. Fange, O. G. Berg, P. Sj¨oberg, and J. Elf, Proceedings of the National Academy of Sciences , 19820 (2010).[55] D. Fange, A. Mahmutovic, and J. Elf, Bioinformatics , 3155 (2012).[56] D. Rapaport, Physical Review Letters , 186101 (2008).[57] M. F. Hagan and D. Chandler, Biophysical Journal , 42 (2006).[58] H. D. Nguyen, V. S. Reddy, and C. L. Brooks III, Nano Letters , 338 (2007).[59] D. Rapaport, Physical Review E , 051917 (2012).[60] I. G. Johnston, A. A. Louis, and J. P. Doye, Journal of Physics: Condensed Matter ,104101 (2010).[61] C. Horejs, M. K. Mitra, D. Pum, U. B. Sleytr, and M. Muthukumar, The Journal of ChemicalPhysics , 125103 (2011).[62] A. W. Wilber, J. P. Doye, A. A. Louis, E. G. Noya, M. A. Miller, and P. Wong, The Journalof Chemical Physics , 085106 (2007).[63] A. W. Wilber, J. P. Doye, and A. A. Louis, The Journal of Chemical Physics , 175101(2009).[64] J. Largo, F. W. Starr, and F. Sciortino, Langmuir , 5896 (2007).[65] P. F. Damasceno, M. Engel, and S. C. Glotzer, Science , 453 (2012).[66] X. Liu, J. C. Crocker, and T. Sinno, The Journal of Chemical Physics , 244111 (2013).[67] F. Sciortino, C. De Michele, and J. F. Douglas, Journal of Physics: Condensed Matter ,155101 (2008).[68] J. D. Halverson and A. V. Tkachenko, Physical Review E , 062310 (2013).[69] M. Eigen, Diffusion control in biochemical reactions. In Quantum Statistcal Mechanics in theNatural Sciences , edited by S. L. Mintz and S. M. Widmayer, Studies in the Natural Sciences,Vol. 4 (Plenum Press, New York and London, 1974) pp. 37–61.[70] D. Shoup and A. Szabo, Biophysical Journal , 33 (1982). [71] C. Korn and U. Schwarz, The Journal of Chemical Physics , 095103 (2007).[72] J. Schluttig, D. Alamanova, V. Helms, and U. S. Schwarz, The Journal of Chemical Physics , 155106 (2008).[73] J. Schluttig, C. B. Korn, and U. S. Schwarz, Physical Review E , 030902 (2010).[74] J. Garc´ıa de la Torre, M. L. Huertas, and B. Carrasco, Biophysical Journal , 719 (2000).[75] B. Berger, P. W. Shor, L. Tucker-Kellogg, and J. King, Proceedings of the National Academyof Sciences , 7732 (1994).[76] T. L. Hill, An Introduction to Statistical Thermodynamics (Dover Publications, New York,1986) (Unabridged, corr. republ. of the 2. print., Reading, Mass., 1962).[77] M. Pogson, R. Smallwood, E. Qwarnstrom, and M. Holcombe, Biosystems , 37 (2006).[78] M. T. Klann, A. Lapin, and M. Reuss, BMC Systems Biology , 71 (2011).[79] J. Paijmans, The Fundamental Lower Bound of the Noise in Transcriptional Regulation (Mas-ter thesis, University of Amsterdam, Amsterdam, 2012).[80] N. Kern and D. Frenkel, The Journal of Chemical Physics , 9882 (2003).[81] A. Arkhipov, Y. Yin, and K. Schulten, Biophysical Journal , 2806 (2008).[82] A. Frost, V. M. Unger, and P. D. Camilli, Cell , 191 (2009).[83] H. D. Nguyen and C. L. Brooks III, Nano Letters , 4574 (2008).[84] O. M. Elrad and M. F. Hagan, Nano Letters , 3850 (2008).[85] H. D. Nguyen, V. S. Reddy, and C. L. Brooks III, Journal of the American Chemical Society , 2606 (2009).[86] L. Blanchoin, R. Boujemaa-Paterski, C. Sykes, and J. Plastino, Physiological Reviews ,235 (2014).[87] D. L. D. Caspar and A. Klug, in Cold Spring Harbor Symposia on Quantitative Biology , Vol. 27(Cold Spring Harbor Laboratory Press, 1962) pp. 49–50.[88] D. Kitagawa, I. Vakonakis, N. Olieric, M. Hilbert, D. Keller, V. Olieric, M. Bortfeld, M. C.Erat, I. Fl¨uckiger, P. G¨onczy, and M. O. Steinmetz, Cell , 364 (2011).[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, andB. Zuber, Science , 1196 (2011).[90] P. Guichard, V. Hachet, N. Majubu, A. Neves, D. Demurtas, N. Olieric, I. Fluckiger, A. Ya-mada, K. Kihara, Y. Nishida, S. Moriya, M. O. Steinmetz, Y. Hongoh, and P. G¨onczy, Current Biology , 1620 (2013).[91] C. Erlenk¨amper and K. Kruse, The Journal of Chemical Physics , 164907 (2013).[92] R. Kapral, The Journal of Chemical Physics138