Allosteric control in icosahedral capsid assembly
AAllosteric control in icosahedral capsid assembly
Guillermo R. Lazaro and Michael F. Hagan ∗ Martin Fisher School of Physics, Brandeis University, Waltham, MA, 02454 (Dated: October 23, 2018)During the lifecycle of a virus, viral proteins and other components self-assemble to form an ordered proteinshell called a capsid. This assembly process is subject to multiple competing constraints, including the need toform a thermostable shell while avoiding kinetic traps. It has been proposed that viral assembly satisfies theseconstraints through allosteric regulation, including the interconversion of capsid proteins among conformationswith different propensities for assembly. In this article we use computational and theoretical modeling to explorehow such allostery affects the assembly of icosahedral shells. We simulate assembly under a wide range ofprotein concentrations, protein binding affinities, and two different mechanisms of allosteric control. We findthat, above a threshold strength of allosteric control, assembly becomes robust over a broad range of subunitbinding affinities and concentrations, allowing the formation of highly thermostable capsids. Our results suggestthat allostery can significantly shift the range of protein binding affinities that lead to successful assembly, andthus should be accounted for in models that are used to estimate interaction parameters from experimental data.
I. INTRODUCTION
The assembly of a virus outer protein shell (capsid) requiresa delicate balance among thermodynamic and kinetic con-straints. The proteins must assemble quickly to evade proteol-ysis and detection by the host, and their capsid must be suffi-ciently thermostable to survive intact under potentially harshconditions while searching for a new infection target. Yet,self-assembly of ordered structures usually requires weak, re-versible interactions among the components, since strong in-teractions lead to kinetic traps [1–3]. Many viruses must alsocontrol the time and place of assembly, so that the capsid canselect specific components from amidst a crowded cellular en-vironment. A number of strategies have been proposed bywhich viruses control their assembly process to ensure pro-ductive and timely capsid formation ( e.g. [4–8]). One suchstrategy is allosteric regulation [8, 9], in which capsid pro-teins sample an ensemble of conformational states with dif-ferent propensities for assembly, with the relative populationsof different states influenced by binding of proteins or othermolecules. In this article, we theoretically and computation-ally examine how allostery at the level of protein-protein in-teractions can lead to self-regulation of assembly kinetics.Learning mechanistic information such as allosteric regula-tion from experiments alone is challenging because most as-sembly intermediates are transient, and thus not readily ob-served. For example, a wealth of information has been ob-tained from in vitro experiments in which capsid assembly ki-netics are monitored by size exclusion chromatography (SEC)or x-ray or light scattering (e.g.) [10–17]. However, sinceintermediates are usually undetectable, these bulk techniquesprimarily report on the concentrations of assembled capsidsand unassembled subunits. Recently, techniques which mon-itor individual capsids or can detect transient intermediateshave begun to address this limitation [18–27]. However, eventhese techniques provide structural data at limited resolutionand cannot characterize the full ensemble of intermediates. ∗ [email protected] Therefore, theoretical models are needed to obtain a com-plete understanding of capsid assembly from such experimen-tal data.Theoretical models have already played an important rolein relating experimental data to assembly pathways and thedriving forces that control them. For example, binding affini-ties have been estimated by fitting the ratio of assembled cap-sids to unassembled subunits measured at long times to theequilibrium law of mass action [28, 29]. Assembly kineticshave been analyzed using models in which capsid formationis viewed as the assembly of rigid subunits (lacking internaldegrees of freedom) into polyhedral shells. These models canbe formulated as a master equation and solved numerically( e.g. [11, 30–33]) or analyzed by stochastically generatingtrajectories consistent with the master equation ( e.g. [34–40]).Despite simplifications used to make these methods tractable,they capture many features of experimental assembly kinetics.Fitting their results against light scattering data has enabledestimates of physical parameters such as subunit-subunit bind-ing affinities and rate constants [11, 14, 41–43]. These re-sults suggest that capsid subunit-subunit binding affinities aregenerically weak, on the order of 4 kcal/mol [3, 28]. Whenbinding affinities exceed this limit, assembly is limited by ki-netic traps, in which the formation of long-lived disordered orpartially assembled structures inhibit capsid formation.Despite these important insights, current models do notcapture all aspects of experimental data. Typically, mod-els cannot quantitatively reproduce kinetics at all timescalesacross a range of concentrations [41]. Similarly, capsid disas-sembly exhibits a surprising degree of hysteresis consideringthe measured weakness of subunit binding affinities [44, 45].For example, assembled capsids are highly stable in infinitedilution; e.g. hepatitis B virus (HBV) capsids exhibit virtu-ally no subunit dissociation even on a timescale of months[46]. While this observation could be accounted for by a post-assembly maturation process that increases capsid stability, asyet there is no evidence for this in HBV.One feature which has not yet been incorporated in most as-sembly models is that capsid proteins sample multiple confor-mational states, with different propensities for assembly. Forexample, structural studies on HBV [9], brome mosaic virus a r X i v : . [ q - b i o . S C ] A p r (BMV) and HIV [47] find that their capsid proteins adoptconformations in solution that are incompatible with insertioninto a capsid, suggesting that the protein’s primary conforma-tion in solution is inactive for assembly (assembly-inactive),and that capsid formation requires a transition to an assembly-active conformation. These observations have led to the hy-pothesis that assembly of icosahedral viruses may be subjectto allosteric regulation [9, 28, 41, 47–50].It has long been appreciated that protein conformationalswitching can play a key role in the kinetics of assemblingfilamentous structures. Asakura [51] showed that the elon-gation rate of the bacterial flagellum is limited by the rate offlagellin monomer undergoing a conformational transition toits bound form. From this observation, he deduced that in-teraction with the flagellum triggers the conformational tran-sition in the monomer by an ‘induced fit’. Similar observa-tions were made for other filamentous assemblies, includingthe bacteriophage T4 tail [52] and tobacco mosaic virus [53–56]. Caspar coined the term ‘autostery’ [8, 57] to describesuch an induced fit process, in which a protein that existsin an assembly-inactive conformation in solution is driven toswitch to an assembly-active conformation by interacting withother copies of itself within an assemblage. By controllingthe rate of nucleation in comparison to elongation, autosterycould provide mechanism for self-regulation of assembly.The role of an assembly-active/inactive transition in the as-sembly kinetics of icosahedral viruses has received less atten-tion; more work has focused on the roles of conformationalswitching in overcoming the geometrical constraints imposedby an icosahedral geometry [58–61] and structural polymor-phism [19, 62–65]. Although Packianathan et al. [9] wasprimarily an experimental study, they also used a rate equa-tion approach to compare the behaviors of assembly with noallosteric regulation, with induced fit (autostery), or allosterywithout induced fit. They concluded that the induced fit mech-anism leads to productive assembly and could increase hys-teresis associated with disassembly, but that allostery withoutinduced fit did not lead to productive assembly, because suchstrong subunit affinities were required that kinetic trapping re-sulted.In this article we perform a more extensive theoretical andcomputational investigation of the effects of allostery on as-sembly kinetics and their sensitivity to kinetic trapping. Incontrast to Ref. [9], we find that both mechanisms of allostericregulation (with or without induced fit) can drive productiveassembly, although induced fit allows for productive assem-bly over a wider parameter range. Under moderate param-eter values, allostery does not enhance assembly robustness— the width of the range of subunit concentrations or bind-ing affinities leading to productive assembly is not increased.However, under sufficiently strong allosteric control (mean-ing that the population of unassembled subunits is stronglyshifted toward the assembly-inactive conformation), high as-sembly yields are achieved over a broad range of parameters,including high binding affinities. Our results highlight the im-portance of accounting for allostery in models used to estimateparameters from experimental data.This paper is organized as follows. In the next section, we incorporate allostery into a computational model and a mas-ter equation description of assembly. In section ‘Scaling es-timates for the effect of allostery on assembly timescales’ wedevelop scaling estimates for the effects of allostery on assem-bly timescales and assembly robustness. Then, in the Resultswe test these scaling estimates against numerical results fromthe computational and master equation models. In the Dis-cussion we summarize the key observations, identify potentialfurther extensions to the models, and discuss implications forestimating parameters by fitting against experimental data. Fi-nally, the appendix provides further details about the modelsand compares the computational and master equation models. II. MODELSA. Allostery Models
To represent allostery in our models, unassembled sub-units interconvert between assembly-inactive and active con-formations, with equilibrium constant K A = exp( − g A /k B T ) and g A > the unfavorable free energy associated with amonomer adopting the active state. The case of no allosterycorresponds to the limit K A → ∞ . For simplicity, here andthroughout, we focus on the limit where the active/inactive in-terconversion rate is fast relative to assembly timescales. Wediscuss effects of a slow conformational interconversion stepin the Discussion.The two allostery behaviors that we consider throughout thepaper are illustrated by the schematics in Fig. 1. In the first,there is allostery but no induced fit, meaning that only subunitsin the active conformation can associate with any size interme-diate. In the second, subunits in the inactive conformation canassociate with intermediates equal to or larger than a nucleuscomprising n nuc subunits. For notational consistency, we referto these cases respectively as ‘No Induced Fit’ and ‘InducedFit’, and the case without allostery (all subunits are active) as‘No Allostery’.Structural evidence suggests that the active-inactive transi-tion in the HBV protein is closely linked to formation of astable critical nucleus [9]. The latter is associated with thegeometry of an icosahedral shell. Since smaller intermediateshave fewer interactions per subunit than large ones, stabilityincreases with intermediate size. The critical nucleus refers tothe smallest intermediate (or ensemble of intermediates) fromwhich assembly into a complete capsid is more probable thancomplete dissociation. In vitro experiments on several viruseshave determined critical nucleus sizes corresponding to smallpolygons ( e.g. a pentamer or trimer of dimers) [2]. Similarly,in the computational model described below we find that un-der certain conditions the critical nucleus corresponds to thesmallest polygon that can form, which is a pentagon as shownin Fig. 2C (see also Refs. [66, 67]). Therefore, we considerthe critical nucleus as the minimum seed capable of driving aconformation change. I.e., inactive subunits can bind to an in-termediate at or above the critical nucleus size, but only activesubunits can associate to pre-critical intermediates (Fig. 1).We now describe the two modeling approaches we use to Assembly subunit
Induced Fit K A Inactive Active Critical nucleus ElongationNucleation ... No Induced Fit ...
FIG. 1. The two mechanisms of allostery that we consider in this article. The irreducible assembly unit (a protein dimer in this schematic) inter-converts between assembly-inactive and assembly-active conformations, with equilibrium constant K A . In both mechanisms, only assembly-active subunits can combine to form intermediates smaller than a critical nucleus. In the ‘Induced Fit’ mechanism, both assembly-active andassembly-inactive subunits can bind to larger intermediates, whereas in the ‘No Induced Fit’ mechanism, only assembly-active subunits canbind throughout the assembly process. study the effects of allostery on assembly. B. Computational model (A) (B) (C)
FIG. 2. Computational model geometry. (A)
Geometry of two in-teracting subunits with the bond vectors depicted as arrows and at-tractors colored teal. The angle between each of the subunit bondvectors is o , and the interactions are described in the Appendix(Eqs. (16)). The dihedral angle φ bij is not shown. (B) The completecapsid, which contains 20 subunits arranged with icosahedral sym-metry. (C)
Critical nucleus for binding of assembly-inactive subunitsin the ‘Induced Fit’ model.
We consider a model for the assembly of icosahedral shellsused in Ref. [68], in which subunits are modeled as rigid bod-ies, with excluded volume interactions represented by spher-ically symmetric repulsive forces, and the complementarysubunit-subunit interactions that drive assembly representedby directional attractions. The lowest energy states in themodel correspond to ‘capsids’, which consist of N = 20 sub- units (each of which could represent a protein trimer) in a shellwith icosahedral symmetry. Because the spatial positions andorientations of all subunits are explicitly tracked, there are noassumptions about assembly pathways or the structures thatemerge from assembly.Following the approach developed by Schwartz et al. [60],the subunit excluded volume is spherically symmetric andthree attractive patches (bond vectors) are rigidly fixed to thesubunit, with each pair of bond vectors forming an angle of ◦ (see Fig. 2 and Eq. (16)). There is a favorable interactionbetween subunits when (1) the ends of bond vectors nearlyoverlap, (2) the bond vectors are nearly anti-parallel, and (3)the secondary bond vectors are nearly coplanar. Twenty sub-units realizing these conditions results in the minimum energytarget structure (a complete capsid) shown in Fig. 2. The in-teraction strength is tuned by the parameter (cid:15) b . All interactionpotentials are pairwise.The subunit-subunit interaction follows additional rules de-pending on the type of allostery being modeled. For ‘NoAllostery’, all pairs of subunits meeting the interaction cri-teria listed above experience attractive interactions. For ‘NoInduced Fit’, subunits stochastically switch between inactiveand active conformations. Only pairs of active subunits sat-isfying the binding criteria experience attractive interactions.Since there is no autostery in this case, the effects of allosteryon assembly arise only due to the cooperativity of capsid as-sembly. For ‘Induced Fit’, we define the threshold for au-tostery activity to be the formation of at least one polygon(any closed cycle of subunit-subunit interactions, such as apentamer). Thus, in this case interactions are possible for anypair of active subunits, and also between inactive subunits andsubunits within a partial capsid which has at least one com-plete polygon. Full details of the model are given in the ap-pendix. Simulation parameters.
The parameters of the model arethe energy associated with the attractive potential, (cid:15) b , and thespecificity of the directional attractions, which is controlledby the angular parameters θ c and φ c . Subunit positions andorientations were propagated using overdamped Brownian dy-namics according to a second order predictor-corrector algo-rithm [69, 70], with the unit of time t = σ /D , where D is the subunit diffusion coefficient and σ is the subunit diam-eter. We simulated systems with 500 subunits in a periodicbox with side length 17, where all distances are measuredin units of the subunit diameter σ . For each parameter set,results were averaged over 20 or more independent simula-tions. The orientational specifity parameters were θ c = 0 . and φ c = π . These parameters tend to disfavor the forma-tion of incorrect subunit-subunit interaction geometries andthus inhibit formation of malformed capsids, allowing us tostudy the high affinity limit. However, such configurations doarise at higher binding energies as discussed below. To obtaindimensionless units, we rescale energies by k B T and timesby t . The subunit conformational switching rate was . /t .Simulations were initialized by generating random positionsand orientations for subunits, with subunit positions that ledto subunit-subunit overlap (positive potential energies in ex-cess of 1 k B T ) rejected. C. Master equation model for capsid assembly
We also consider a master equation description of poly-hedral shell assembly, which is sufficiently computationallytractable to allow modeling assembly over broad parameterranges. Specifically, we extend the ‘nucleation and growth’model described in Ref.[66] to include multiple subunit con-formations. This model was based on the work of Zlotnickand coworkers [11, 30, 31], and consists of a system of cou-pled rate equations that describe the time evolution of concen-trations of empty capsid intermediates: dc dt = − k c + 2¯ k c + N (cid:88) n =2 − k n c n c + ¯ k n c n dc n dt = k n − c c n − − k n c c n n = 2 . . . N − ¯ k n c n + ¯ k n +1 c n +1 (1)where c n is the concentration of intermediates with n sub-units, and k n and ¯ k n are respectively association and dissoci-ation rate constants for intermediate n . The initial conditionis c (0) = c , c n (0) = 0 for n > . The extensions of Eq. (1)to describe allostery are given in the appendix.There are several important assumptions underlying themaster equation: malformed capsids are not considered [60,62–64, 71–74], assembly proceeds along a single pathway[34, 40, 75], only single subunits can bind or unbind, and onlyone k n and ¯ k n are considered for each size n (averaged over all intermediates of that size). However, these assumptionsare not present in the BD simulations described above,andwe find close agreement between the two approaches (see ap-pendix Fig. 11). Moreover, Ref. [66] showed that extendingEqs. (1) to relax these simplifications does not qualitativelychange their predictions. Most importantly, rate equations ofthis form capture many features of experimental assembly ki-netics data ( e.g. [11, 76]).The association and dissociation rate constants are re-lated by detailed balance ¯ k n = k n exp(∆ g n /k B T ) /v , with ∆ g n = G n − G n − the change in free energy due to asso-ciation of a subunit. Association free energies ∆ g n , whichcan be fit to experimental data using the law of mass ac-tion [28, 77, 78], include hydrophobic and electrostatic inter-actions [79] and depend on p H and salt concentration [28].Specifying the assembly model requires defining the set of in-termediates n and the transition rates { k n , ¯ k n } .The model we use here is based on those of Zlotnick andcoworkers [11, 30, 31], in which the subunit-subunit asso-ciation free energy for intermediate n is proportional to thenumber of new subunit-subunit contacts n c, n formed by addi-tion of a subunit to that intermediate [80]. Specifying { ∆ g n } thus requires defining the geometry of each intermediate. Thisusually begins with specifying the geometry of a capsid andits subunits in terms of a polyhedron (for example, see Fig. 2or Fig. 1 in Ref. [31]), and assuming that assembly proceedsalong a single path. The path can be comprised of the low-est energy intermediate for each size n [30] or correspond toan ‘average’ pathway [31] in which all subunits, except dur-ing the initial and final stages of assembly, make the sameaverage number of contacts. We choose the latter definition,since it is simpler and both definitions lead to a qualitativelysimilar behavior. Specifically, the association rate constant k is independent of intermediate size and association free en-ergies are given by ∆ g n = g b before nucleation ( n < n nuc )and ∆ g n = g elong during elongation ( n nuc ≤ n < N − ),where n nuc is the critical nucleus size. Finally, inserting thelast subunit makes the maximum number of possible contactsand thus enjoys the most favorable association free energy g N .Because this model was previously explored extensively inthe absence of conformation changes [66], we focus on oneset of interaction parameters (except where noted otherwise):capsid size N = 120 corresponding to 120 dimer subunits inhepatitis B virus [28], critical nucleus size n nuc = 5 (a pen-tamer of dimers), association rate constant k = 10 M − s − [76], and free energy parameters g b = 7 k B T ( ≈ kcal/mol)[28], g elong = 2 g b and g N = 2 g elong , which imply that addinga subunit becomes on average twice as favorable after nucle-ation and four times as favorable for the final subunit. III. SCALING ESTIMATES FOR THE EFFECT OFALLOSTERY ON ASSEMBLY TIMESCALES
To gain an intuitive understanding of how allostery can af-fect assembly, in this section we derive simple scaling esti-mates for the timescales associated with the two assemblymechanisms shown in Fig. 1, based on the master equationmodel (Eqs. 1). Although we introduce a number of simplifi-cations, in the next section we show that the resulting scalingestimates apply at least qualitatively when these simplifica-tions are relaxed in the computational and theoretical models.We closely follow Ref. [66], except that we extend the analy-sis to include allostery.We consider a system of capsid protein subunits with to-tal concentration c that assemble into capsids with N sub-units. The word subunit refers to the basic assembly unit,which could be a protein dimer or larger oligomer [2]. As inthe master equation model, we break the assembly of a capsidinto ‘nucleation’ and ‘elongation’ phases. For simplicity weassume that the association rate constant k is independent ofintermediate size, so that for the ‘No Allostery’ reaction ratesof association to each intermediate are kc with c the concen-tration of free subunits. We assume the limit of fast conforma-tional interconversion, so for the ‘No Induced Fit’ case, asso-ciation rates are given by kf A c with f A = K A / (1 + K A ) .For ‘Induced Fit’, association rates are kf A c for intermediatesize n < n nuc and kc for n ≥ n nuc .We now write the time required for an individual capsid toassemble as τ = τ nuc + τ elong with τ nuc and τ elong the averagetimes for nucleation and growth, respectively. The elongationtimescale can be estimated by the mean first passage time fora biased random walk with a reflecting boundary conditionsat n nuc and absorbing boundary conditions at N , with forwardand reverse hopping rates given by the subunit association anddissociation rates respectively [66]. For early in the reaction,when c ≈ c , this results in τ elong ∼ = N/kc f m A A (2)where m A = 1 for ‘No Induced Fit’ and m A = 0 for ‘InducedFit’. We see that the elongation time is equal for ‘No Al-lostery’ and ‘Induced Fit’ since in the latter case all subunitscan bind to post-nucleated intermediates. Ref. [66] showedthat the duration of the lag phase in light scattering is propor-tional to τ elong , thus predicting that the lag time should scaleinversely with subunit concentration. This prediction was re-cently confirmed by experiments on HBV assembly [81].The mean nucleation time at the beginning of the reactioncan be estimated from the statistics of a random walk biasedtoward disassembly [66]. Including conformation dynam-ics results in τ minnuc ≈ k − exp ( G ˆ n /k B T ) c − n nuc f − n nuc A , where ˆ n = n nuc − so that G ˆ n is the interaction free energy of thestructure just below the critical nucleus. This estimate can beunderstood by noting that the pre-critical nucleus is presentwith concentration c ˆ n ∼ = exp( G ˆ n ) c ˆ n f ˆ n A , and the rate of ac-tive subunits associating to the precritical nucleus is given by kc f A .However, because free subunits are depleted by assembly,the nucleation rate never reaches this value, and net nucle-ation asymptotically approaches zero as the concentration ofcompleted capsids approaches its equilibrium value. Thus, weestimate the median assembly time τ (the time at which thereaction is 50% complete) by treating the system as a two-state reaction with n nuc -th order kinetics, which yields [66] τ ∼ = A / P N N k exp ( G ˆ n /k B T ) c − ˆ n f − ˆ n − A (3) with A / = ˆ n − n , and P N as the equilibrium fraction ofsubunits in complete capsids. The factor of N − in Eq. 3accounts for the fact that N subunits are depleted by each as-sembled capsid.When capsid growth is fast compared to nucleation, the ex-pressions Eq. 2 and Eq. 3 respectively predict the duration ofthe lag phase and the median assembly time. However, theserelations begin to fail as intermediate concentrations build upabove a crossover concentration c c at which the initial elonga-tion and nucleation times are equal c c ∼ = N − n e G ˆ n ˆ nk B T f − ˆ n − m A +1ˆ n A (4)Significant kinetic trapping then sets in above a threshold con-centration set by τ elong = τ : c kt ∼ = (cid:0) A / N − (cid:1) n − e G ˆ n (ˆ n − k B T f − ˆ n − m A +1ˆ n − A . (5)While the above analysis identifies a maximum concentra-tion above which the reaction will become kinetically trapped,we can also identify a minimum concentration below whichthe median assembly time becomes longer than the maximumtimescale of the experiment or simulation τ max c min ∼ = (cid:0) A / /τ max N k (cid:1) − n e βG ˆ n / ˆ n f − ˆ n +1ˆ n A (6)which applies to both ‘No Induced Fit’ and ‘Induced Fit’.A similar analysis can be performed for fixed subunit con-centration and varying g b , corresponding to changing p H orsalt concentration. For example, if we assume that pre-criticalintermediates with ˆ n subunits have ˆ n − subunit-subunit in-teractions (Fig. 1), then G ˆ n = (ˆ n − g b and the minimumaffinity for assembly in finite time τ max is shifted according to g shiftb = g b + ˆ n + 1ˆ n − f A . (7) A. The influence of conformation dynamics on assemblyrobustness
Based on the above scaling estimates, we now examinewhether introducing allostery makes assembly more robust — i.e. , at a finite timescale τ max relevant to an experiment or acell, does allostery increase the range of concentrations overwhich productive assembly occurs? Specifically, we assumethat kinetic trapping precludes productive assembly within τ max , and consider the ratio between the minimum concentra-tion leading to significant nucleation, and the maximum con-centration leading to kinetic trapping, δc = c kt /c min . Sincenucleation with allostery must occur within τ max , we increasethe binding affinity according to Eq. (7), resulting in δc ∼ = exp( βg shiftb / ˆ n ) f m A ˆ n − A . (8)Since m A = 0 for ‘Induced Fit’, allostery with induced fit hasno effect on the range of concentrations over which assemblyoccurs (for moderate interaction strengths), while in the ab-sence of induced fit, allostery renders assembly kinetics moresensitive to concentration (since f A < ). However, eitherallostery mechanism can significantly increase the maximumthermostability of a capsid that can be achieved without ki-netic trapping. For a fixed concentration c , equating Eqs. (2)and (3) shows that the minimum binding free energy g ktb be-low which kinetic trapping occurs is decreased (higher bind-ing affinity) according to g ktb = g ktb ( f A = 1) + ˆ n + 1 − m A ˆ n − f A . (9)The maximum kinetically accessible thermostability of acomplete capsid is then controlled by the free energy per sub-unit, g sub = n B / g ktb with n B the number of contacts eachsubunit makes with its neighbors. For the case we considerbelow with n nuc = 5 and n B = 4 , g sub increases by a fac-tor n B m c /g A / ≈ . g A . Thus, even a modest activationenergy could substantially increase the caspid stability. Thiseffect, together with the asymptotic approach to equilibriuman assembly reaction discussed below, could contribute to ob-servations of unexpectedly large hysteresis between capsid as-sembly and disassembly. Strong allostery drives robust assembly in the high-affinity limit.
While the above analysis shows that allosterydoes not lead to more robust assembly for moderate inter-action strengths, we observe dramatically different behaviorwith strong interactions, g b (cid:29) k B T . It is well-known that un-regulated assembly fails due to kinetic trapping in this limit.Since subunit-subunit interactions are effectively irreversible,the effective critical nucleus size is reduced to a dimer regard-less of the capsid geometry, and elongation becomes slowerthan nucleation for any feasible parameters. However, equat-ing Eqs. (2) and (3) shows that this trap can be avoided by asufficiently high g A . We expect trapping to be avoided whenthe parameter p trap (cid:46) , with p trap ∼ = f − m A A N . (10)This result holds independent of subunit concentration andcapsid geometry, provided the capsid terminates at a finite sizeand the subunit-subunit affinity is strong enough to stabilizethe capsid. It is important to note that assembly would be sen-sitive to formation of defective capsids in this limit, which weneglect in this scaling analysis. However, as shown below wedo observe this limit in Brownian dynamics simulations wheredefective assembly is allowed. IV. NUMERICAL RESULTSA. Effect of Allostery on Assembly Robustness
In this section, we test the scaling predictions against resultsfrom the Brownian dynamics (BD) simulations and masterequation model. For the BD simulations, we focus on quan-tities involving variation of the subunit binding energy pa-rameter (cid:15) b , as varying this parameter is more computationallytractable than varying subunit concentrations over the rangeneeded to test scaling. Specifically, we test how the assembly yield and median assembly time depend on (cid:15) b and the activa-tion energy g A , and we test the prediction that a sufficientlylarge value of g A enables assembly even in the limit of highsubunit-subunit binding affinity. The latter test is particularlyimportant since the BD simulations allow for the formationof defective capsids. We then vary concentration using themaster equation model.We begin by examining the time dependence of the frac-tion of subunits in complete capsids ( P N ) as the binding en-ergy is varied. A capsid is defined as a structure containing20 subunits, each with strong interactions with three neigh-bors. This definition counts only configurations which cor-respond to small fluctuations around the icosahedral groundstate. We note that experimental measurements of assem-bly kinetics commonly employ SEC, which monitors P N , orlight scattering, which under certain conditions monitors themass-averaged molecular weight of assemblies [11]. Ref. [66]showed that P N and assembly molecular weight closely trackeach other below the crossover concentration c c (Eq. (4)), andquantities such as the lag phase and nucleation time followthe same scaling laws when calculated from either observable.Therefore, here we present only P N .Fig. 3A shows P N as a function of time for several valuesof (cid:15) b with the ‘No Allostery’ model. We see that initially as (cid:15) b increases the lag phase shortens, P N rises in time more rapidly,and asymptotes at a higher value. However, at the highest (cid:15) b shown, (cid:15) b = 30 , the fraction capsid quickly saturates ata low value, indicative of kinetic trapping. To illustrate thispoint, Fig. 3B shows P N ( t ) for two values of (cid:15) b , overlaid witha plot of the fraction of subunits in intermediates, P int ( t ) = (cid:80) N − n =2 nc n ( t ) /c . We see that for (cid:15) b = 14 the intermediatespeak near the end of the lag phase, and then rapidly fall ascapsids are produced. In contrast, for (cid:15) b = 30 , nearly allsubunits are trapped in intermediates. Thus capsids are onlyslowly produced when larger intermediates scavenge subunitsfrom smaller ones. Since we find that it correlates well withother measures of kinetic trapping, we use P int to characterizetrapping as a function of parameter values.For comparison, Figs. 3C,D show the time-dependence of P N and P int calculated from the master equation with varyingsubunit concentration [82]. We see that the kinetics and therelationship between P int and the onset of kinetic trapping areconsistent with the BD results.Next, we consider how assembly robustness changes whenallostery is introduced. Fig. 4A shows P N measured at longbut finite times as a function of (cid:15) b for the three conformationcases with g A = 4 . We see that in both cases, allostery shiftsthe onset of assembly to dramatically higher binding affin-ity. More significantly, assembly stays remarkably productiveup to extremely high values of (cid:15) b for both allostery mecha-nisms. This observation corresponds to the high affinity limitdiscussed above (Eq. (10)). Analysis of simulation trajecto-ries shows that the eventual decline in P N at high (cid:15) b arises dueto the formation of defective capsids.To show how assembly robustness depends on the strengthof allostery, we compare P N calculated for the ‘No Allostery’with ‘Induced Fit’ at two values of g A in Fig. 4B. To aid incomparing these cases, we have ‘undone’ the shift in the bind- t / 10 t P N (A) b =16 b =14 b =30 t / 10 t P N , P i n t =30 b =14 (B) P N t [sec] μ M20 μ M40 μ M60 μ M(C) P N , P i n t t [sec] (D) 100 μ M40 μ M FIG. 3. The dependence of assembly kinetics on parameter values for ‘No Allostery’. (A)
The fraction of subunits in complete capsids P N observed in Brownian dynamic simulations (BD) is shown as a function of time for indicated values of the binding energy parameter (cid:15) b . (B) The fraction capsid (solid lines) is compared to the fraction of subunits in intermediates P int (dashed lines) for small and large values of (cid:15) b .For (A) and (B), each simulation is run until . × t , and each data point corresponds to an average over 20 independent simulations. (C) The fraction capsid as a function of time measured for the master equation is shown at indicated total subunit concentrations c , with bindingaffinity g b = − (all energies are in units of k B T ). (D) The fraction capsid is compared to the intermediate fraction for the master equation,for concentrations below and above the trapping point c kt = 60 µ M (Eq. 5). ing energies suggested by Eq. (7) according to (cid:15) effb = (cid:15) b + ˆ n + 1ˆ n − f A , (11)so that the onset of productive assembly occurs roughly at thesame value of (cid:15) effb for each case. We note that the curves do notline up perfectly because this mapping is only approximate,since the free energy includes binding entropy factors whosemagnitude depend logarithmically on (cid:15) b (see Eq. (12) in theappendix) and the degeneracy of available binding sites on agiven structure [71, 83, 84]. We see that the smaller activationenergy g A = 2 leads to assembly over approximately the samerange of (cid:15) effb as for the ‘No Allostery’ case, supporting ourearlier conclusion that, below the high affinity limit, allosterydoes not significantly increase assembly robustness.To further understand the onset of the high affinity limit androbust assembly, we examine the dependence of kinetic trap-ping on parameter values. First, to support the earlier state-ment that P int is a good metric for the extent of kinetic trap-ping, Fig. 5 compares P N and the intermediate fraction P int measured at long times for ‘No Allostery’; we see that a risein P int correlates with the decline in P N .Next we focus on high binding affinity, (cid:15) b = 30 , and show P int in Fig. 6A for different values of g A , plotted against thetrapping parameter p trap (Eq. (10)). We see that for ‘No In-duced Fit’ the onset of trapping occurs when p trap is of orderone as expected, but trapping for ‘Induced Fit’ is shifted tolower values of p trap . Analysis of simulation trajectories sug- gests that this difference arises due to the formation of mal-formed capsids, which are are more prevalent in the ‘InducedFit’ trajectories. We can understand this observation by not-ing that the probability of a defective capsid increases withthe rate of subunit addition, which is larger for ‘Induced Fit’(since inactive subunits can bind to growing capsids, and be-cause a given value of p trap corresponds to higher f A for ‘In-duced Fit’ in comparison to ‘No Induced Fit’). As further ev-idence for the importance of defective capsids, Fig. 6B showsthe equivalent plot from the rate equation model, which doesnot allow for malformed capsids. In that case, we see that re-sults from both allostery mechanisms fall approximately onthe same curve. We have also tested the predicted scalingagainst capsid size N in that figure by including results fromrate equation calculations with N varying from 20 to 240.Finally, we consider the effect of allostery on assemblytimescales. Fig. 7A shows the median assembly time τ (defined by P N ( τ ) = 1 / ) calculated from the master equa-tion as a function of initial subunit concentration for eachconformation case, with modest interaction parameter values( g b = − , g A = 2 ). This plot allows testing of several scal-ing predictions. Firstly, we see that both in the presence andabsence of allostery, the τ at low concentrations scales as c n nuc − as predicted by Eq. (3). As expected, median timesare equal for ‘No Induced Fit’ and ‘Induced Fit’ at low con-centrations where the reaction is nucleation-limited, since theautostery does not affect nucleation in our model. Secondly,the dashed vertical lines show the point of kinetic trapping c kt / k B T
10 20 30 40 50 60 P N b Induced FitNo AllosteryNo Induced Fit / k B T
10 15 20 25 30 P N beffb No Allostery Induced Fit g A =2Induced Fit g A =4 FIG. 4. (top) Capsid fraction P N as a function of the binding energyparameter (cid:15) b in Brownian dynamics (BD) simulations for the threeconformation cases, with g A = 4 k B T . (bottom) The capsid fractionis shown as a function of the shifted binding energy (cid:15) effb (Eq. 11) for‘No Allostery’ and ‘Induced Fit’ with two activation energy values.For both panels, other parameters are as in Fig. 3. b / k B T
10 20 30 40 50 60 P , P i n t N b P int P N FIG. 5. Capsid fraction P N and intermediate fraction P int as a func-tion of the binding energy (cid:15) b in Brownian dynamics (BD) simulationsfor ‘No Allostery’. predicted by Eq. (5) for each case. We see that the medianassembly time takes off above these threshold concentrations,indicating the presence of trapping.Notice that even a modest conformation energy g A = 2 shifts the region of productive assembly to unrealistically highsubunit concentrations, due to the scaling of nucleation timeswith f A . However, this effect can be counteracted by increas-ing the subunit binding affinity. To illustrate this point, and totest the predicted effect of conformation specificity on assem-bly robustness (Eq. (8)), we also calculated median assemblytimes with the binding affinity for the allostery cases increasedaccording to Eq. (7) (Fig. 7B). As anticipated, the nucleationtimes become equal at low concentrations, and the ‘InducedFit’ and ‘No Allostery’ cases enjoy the same range of produc-tive assembly, while the ‘No Induced Fit’ case becomes ki-netically trapped at a lower concentration due to its increasedelongation time relative to nucleation. To further test the ex- f A m N -6 -4 -2 P i n t Induced Fit No Induced Fit f A m N -2 -1 P i n t NC, N =20NC, N =30NC, N =60NC, N =120NC, N =240NIF, N =20NIF, N =30NIF, N =60NIF, N =120NIF, N =240IF, N =20IF, N =30IF, N =60IF, N =120IF, N =240 FIG. 6. (top)
Kinetic trapping in BD simulations. The magnitudeof trapping P int measured in BD simulations is shown in the limit ofvery high subunit binding strength ( (cid:15) b = 40 ). Results are shown as afunction of the high affinity trapping parameter p trap = f m A N , with m = 1 or respectively for the ‘No Induced Fit’ and ‘Induced Fit’cases (see Eq. (10)) and N = 20 . The parameter f A was tuned byvarying g A /k B T from 1 to 6. Other parameters are as in Fig. 3. (bot-tom) Relationship between capsid size and sensitivity to trapping inthe high-affinity limit from the master equation. The magnitude oftrapping P int calculated for g b = − is shown for the three confor-mation cases and capsid sizes N ∈ [20 , . The parameter f A wastuned by varying g A /k B T from 0 to 10. pected scaling with activation energy, we show the trappingmeasure P int as a function of subunit concentration normal-ized by c kt for all three cases over a range of g A in Fig. 7C.We see that in all cases trapping takes off as the thresholdconcentration is crossed.Fig. 8 shows the variation of τ with the shifted bindingenergy for the three allostery cases calculated from BD simu-lations. We see that the nucleation-dominated regime, markedby an exponential decrease in τ with increasing (cid:15) effb , roughlyoverlaps for the three cases, suggesting that the mapping is ap-proximately correct. Notice that with allostery, τ becomesconstant at large (cid:15) effb ; this corresponds to the high affinity limit. B. Implications of allostery for parameter estimations fromexperimental data
In this section we compare assembly kinetics with and with-out allostery, and consider ramifications for parameter esti-mations made using models that do not account for allostery.Fig. 9 compares the time-dependence of P N at several con-centrations for the ‘No Allostery’ and ‘Induced Fit’ casescalculated from the master equation, with the subunit bind-ing affinity in the ‘Induced Fit’ case adjusted according toEq. (7). In particular, with g b = − for ‘No Allostery’, weobtain g b = − . for ‘Induced Fit’ with g A = 2 . We see c [ M] / P N eq (A) / Induced FitNo Allostery No Induced Fit c [ M] / P N eq (B) / No AllosteryInduced FitNo Induced Fit c / c kt -2 P i n t NANIF g A =0NIF g A =2NIF g A =4NIF g A =6NIF g A =6IF g A =0IF g A =1IF g A =1.5IF g A =2 (C) FIG. 7. The effect of allostery on assembly timescales calculatedfrom the master equation. (A)
The median assembly time τ isshown as a function of initial subunit concentration for the three al-lostery cases, with the subunit binding affinity g b = − and theactivation energy g A = − . For each case, the vertical dashed lineindicates the point of kinetic trapping c kt calculated by Eq. (5). (B) The median assembly time is shown as a function of concentrationfor the three cases, but with the subunit binding affinity shifted ac-cording to Eq. (7), g b = − for NC and g b = − . for NA andA, with g A = 2 . The vertical lines indicate the predicted locationof c kt for these modified subunit affinities. (C) The maximum frac-tion of subunits found in intermediates, P int , is shown as a functionof initial subunit concentration normalized by the trapping threshold c kt (Eq. 5), for indicated allostery cases and values of the activationenergy g A . that this relationship, derived to relate assembly kinetics inthe nucleation-dominated regime, leads to assembly kineticswhich also closely match throughout early times.The most significant difference between the kinetics withand without allostery appears at a very early times. As thekinetics transition from the lag phase to the rapid productionof capsids, the takeoff is more rapid for the ‘Induced Fit’ case.A similar observation was made by Chen et al. [14] and shownto be consistent with light scattering from assembling BMVviruses. This observation can be understood by noting that,due to the stronger binding affinity in the ‘Induced Fit’ case,kinetics in the elongation phase are closer to irreversible andthus the distribution of lag times is more sharply peaked.Despite its modest impact on the form of kinetics, allostery b / k B T
20 30 40 50 / / P N eq No Allostery No Induced FitInduced Fit b FIG. 8. The median assembly time τ measured in BD simulationsis shown for the three conformation cases as a function of (cid:15) effb , with g A = 4 . Other parameters are as in Fig. 3 except simulations wererun until . × t . P N t [sec] μ M40 μ M60 μ M(A) -0.01 0 0.01 0.02 0.03 0.04 0.05 0 20 40 60 80 100 P N t [sec] Induced FitNo Allostery(B)
FIG. 9. The time-dependence of assembly is compared between the‘No Allostery’ (solid lines) and ‘Induced Fit’ (dashed lines) cases for c = 20 µ M. (A) The fraction capsid is shown for indicated initialsubunit concentrations. (B)
The early time course of assembly isshown for both cases at c = 20 µ M, showing that ‘Induced Fit’exhibits a sharper take off at the end of the lag phase in comparisonto ‘No Allostery’. can dramatically skew quantitative parameter values estimatedfrom experimental data. Fig. 9 demonstrates that a fit againstassembly kinetics not accounting for allostery would underes-timate g b , by an amount proportional to g A (see Eq. (7)). Weemphasize that the shift factor in g shiftb is independent of sub-unit concentration. Thus, testing data fits against multiple sub-unit concentrations does not necessarily identify the presenceof allostery. A binding affinity underestimated at one subunitconcentration (due to not accounting for allostery) would beconsistent with data at other subunit concentrations.A second approach to estimate subunit binding affinities isto measure the ratio of capsid to free subunits at very long0 -2 No AllosteryInduced Fit c [ μ M] t =10 sec t =10 secInduced Fitequilibrium No Allosteryequilibrium P N FIG. 10. The difficulties of detecting allostery from finite-time ex-periments. The fraction of subunits in capsids calculated from themaster equation as a function of total subunit concentration is shownfor and seconds, for ‘No Allostery’ (o symbols) and ‘In-duced Fit’ ( (cid:5) symbols). The equilibrium values of P N are shown foreach case as dashed lines. Parameters are g b = − for ‘No Allostery’and g b = − . and g A = 2 for ‘No Allostery’. times, and fit the results to the equilibrium law of mass ac-tion [28]. As has been pointed out previously [2, 45, 66], P N only asymptotically approaches its equilibrium value, andthus fits assuming equilibrium will tend to underestimate g b .We find that unlooked for allostery can dramatically enhanceunderestimation, since the lack of distinguishability seen inFig. 9 persists for very long times. In Fig. 10 we show thefraction capsid P N as a function of subunit concentration pre-dicted from the master equation with and without allostery,for the same interaction parameters as in Fig. 9. We see thatthe curves overlap perfectly at 1000 seconds, and still nearlyoverlap at seconds (about 4 months). By this time, the ‘NoAllostery’ system has nearly reached its equilibrium, and thusthe equilibrium subunit binding affinity estimated from theseresults would be reasonably accurate. However, for ‘InducedFit’ the binding affinity is underestimated by . k B T . Largervalues of g A lead to more severe underestimation with similarlack of distinguishability between the two cases. V. DISCUSSION AND OUTLOOK
It has been well established that the assembly of rigid sub-units into ordered structures requires weak, reversible inter-actions, as interactions which are strong in comparison to thethermal energy lead to kinetic traps [1–3, 71, 74]. Here, wehave used computational and theoretical models to investigatehow this constraint changes when the subunits have internaldegrees of freedom, such as conformational states, allowingthem to change their capacity for interaction during the as-sembly process.We find that a sufficient bias in the free subunit populationtoward the inactive conformation allows productive assem-bly at very high subunit concentrations or binding affinities.In particular, allostery differentially regulates the rates of nu-cleation and elongation, thus suppressing the kinetic trap that arises when free subunits are depleted before capsid elonga-tion is completed. To our knowledge, such robust assemblyhas not observed in vitro, but this effect could be importantfor reactions such as the assembly of the mature HIV cap-sid which occurs at high concentration within the budded vi-ral particle [85]. However, this mechanism does not providecomplete protection against kinetic traps that arise due to de-fective capsid assembly, as we observed in the Brownian dy-namics simulations at very high (cid:15) b . The prevalence of suchdefects, and hence the actual range available for productiveassembly, will depend on the orientational specificity of thesubunit-subunit interactions [71, 84], which has not yet beenevaluated for specific capsid proteins.Given that increasing g A can qualitatively change assemblyrobustness, it would be of interest to measure the subunit con-formational equilibrium as a function of protein sequence andsolution conditions. We have recently used all-atom simula-tions to estimate a free energy difference of about k B T be-tween two quasi-equivalent conformations (meaning differentconformations found at positions with different local symme-try in the capsid) of the MS2 coat protein. However, we notethat the accuracy of such calculations is necessarily limited byforce field accuracy and the quality of sampling. The simula-tions also found that RNA binding could significantly shiftsthe populations. While similar calculations are possible forthe active/inactive transition, they will be complicated by thefact that the inactive ‘conformation’ can be an ensemble withsignificant structural diversity [9].Our findings have several implications for interpretingmechanisms and estimating interaction parameters from ex-perimental data. Firstly, we find that the kinetics of an as-sembly reaction with allostery are quite difficult to distinguishfrom those of a reaction with no conformation dependence atmoderate parameters (Fig. 9). One commonly used test for thequality of a model is to fit interaction parameters to kineticsat one subunit concentration, and then test their predictionsagainst kinetics measured at other concentrations. However,the binding affinity adjustment Eq. (7) used to match param-eters between the ‘No Allostery’ and ‘Induced Fit’ cases inFig. 9 is independent of concentration, and we find that thekinetics match quite closely even into the regime where ki-netic trapping starts to set in. Thus, a binding affinity under-estimated at one subunit concentration (due to not account-ing for allostery) would also appear consistent with data atother subunit concentrations. While it might be argued thatsuch a close match arises because of the simplicity of therate equation model, parameter estimation from experimen-tal data generally relies on models with similar levels of ap-proximations to enable computational tractability. Moreover,we observe a similar matching of kinetics in the BD simula-tions. Thus, our results suggest that strong emphasis shouldbe placed on matching the very early phases of assembly ki-netics during parameter estimation. The importance of fittingthe lag phase was also suggested by Chen et al. [14], who fur-thermore demonstrated this phase can be monitored by lightscattering with millisecond resolution. Techniques sensitiveto individual capsids ( e.g. [26]) will allow further investiga-tion of early-time kinetics.1As we note above, by shifting productive assembly tohigher binding affinities, allostery increases the maximumcapsid thermostability that is kinetically accessible. Thus, al-lostery, in the form of conformational transitions during as-sembly, may offer an alternative strategy to post-assemblyconformational changes or covalent modifications used bysome bacteriophages to stabilize their capsids ( e.g. [86, 87]).Finally, there are several effects we did not investigate herewhich could lead to additional control over assembly. Wehave focused on the limit in which subunit conformationaltransitions are fast in comparison to assembly timescales.Our approaches are easily generalized to other conformationaltimescales. A preliminary investigation showed that decreas-ing the conformational transition rate allows higher assemblyyields in the high affinity limit for the ‘Induced Fit’ case (sim-ilar to the case of bacterial flagella assembly [51]). While wefocus here on allostery at the level of protein-protein inter-actions, there are many examples in the literature suggestingthat interaction with non-protein components, such as RNA,lipid membranes, or small molecule assembly effectors canexert additional allosteric control on protein conformations(see [6]). Understanding how these multiple regulatory mech-anisms cooperate to control the time, place, and rate of as-sembly will lead to a more complete understanding of virallife cycles, and also may identify new strategies for designinghuman-made assembly systems. VI. APPENDIXA. Comparison between Computational and Master EquationModels
In the main text our analysis of the master equation focuseson capsid sizes and parameters consistent with experiments onHBV [11]. Here we consider parameters roughly correspond-ing to those of the BD simulations to further evaluate the de-gree of similarity between master equation and BD results. Tominimize reliance on data fitting, we calculated the subunit-subunit binding free energy g b as a function of the simulationwell-depth (cid:15) b according to g b ( n ) = − ∆ n c, n (cid:15) b − T ( s b + s c ( n )) where ∆ n c, n is the number of subunit-subunit contacts addedto form a cluster of size n (assuming the lowest energy con-figuration at each size n ), s b ( n ) denotes the translational androtational binding entropy penalties and s c denotes the changein ‘configurational’ entropy associated with degeneracy of thelowest energy configuration for each cluster size n . The bind-ing entropy can be estimated from a saddle point approxima-tion of the partition function for a subunit dimer as [71, 83] s b /k B ≈ −
32 log β∂ U att ( r ) ∂r (cid:12)(cid:12)(cid:12)(cid:12) r =2 / σ −
12 log ( β(cid:15) b ) π θ φ . (12)We note that this approximation provides an accurate descrip-tion of the dimerization free energy but neglects additionalentropy penalties incurred by subunits forming more than onecontact. The number of contacts { ∆ c ( n ) } and configura-tional entropy { s c ( n ) } values for building an icosahedron are
10 15 20 25 30 P N beff / k B T (A) Master equationBD simulations
10 20 30 40 P N beff / k B T (B) Induced Fit, g A =4Induced Fit, g A =2No Allostery FIG. 11. (A)
The fraction capsid P N predicted by BD simulationsand the master equation as a function of the binding energy param-eter (cid:15) b for ‘No Allostery’, for t = 10 t . Parameters for the BDsimulations are as in Fig. 3, and the master equation parameters areset correspondingly (see the text in Appendix VI A). (B) The masterequation results for fraction capsid are shown as a function of theshifted binding energy (cid:15) effb (Eq. 7) for ‘No Allostery’ and ‘InducedFit’ with two activation energy values (compare to the bottom panelof Fig. 4). The only adjustable parameter for the master equation inthese results is the association rate constant k , which was set (by eye)to k = 0 . σ /t . given in Table S2 of Roldao et al. [88]. The shifted bind-ing energy (Eq. (11)) is modified to (cid:15) effb = (cid:15) b + ˆ n +1ˆ n − log f A +3 log (cid:0) (cid:15) effb /(cid:15) b (cid:1) .The subunit association rate k is the only adjustable param-eter we used in this comparison. While k could be estimateddirectly from simulations, its value changes for each clustersize due to different occluded volumes neighboring subunits.For simplicity, we use one average value k = 0 . σ /t ,which we estimated (by eye) by comparing master equationand simulation results for P N ( t ) for several values of (cid:15) b .Although the master equation kinetics do not perfectlymatch BD results, their agreement is reasonable consideringthe approximations in our estimate of g b and the simplifica-tions inherent in the master equation. For example, in Fig. 11we show the master equation results for P N calculated at longtime t as a function of the binding energy parameter andactivation energy. The most significant difference is that theBD simulations result in a small number of completed cap-sids at large (cid:15) b even for ‘No Allostery’, whereas the masterequation does not. We believe this difference arises due tobinding of oligomers in the BD simulations [71] which is notaccounted for in the master equation. We note however thatthe effect of oligomer binding would diminish at smaller sub-unit concentrations.2 B. Rate equation model with conformation changes
Extension of Eqs. (1) to include interconversion of free sub-units between inactive and active conformations results in dc dt = − k c + 2¯ k c + N (cid:88) n =2 − k n c n c + ¯ k n c n (13) + k A c − ¯ k A c dc n dt = (cid:0) k n − c + k ∗ n − c ∗ (cid:1) c n − − ( k n c + k ∗ n c ∗ ) c n − (cid:0) ¯ k n + ¯ k ∗ n (cid:1) c n + (cid:0) ¯ k n +1 + ¯ k ∗ n +1 (cid:1) c n +1 for n = 2 . . . N (14) dc ∗ dt = N (cid:88) n =2 − k ∗ n c n c ∗ + ¯ k ∗ n c n − k A c ∗ + ¯ k A c (15)where c ∗ and c are respectively the concentrations ofinactive and active free subunits, k A and ¯ k A are theactive/inactive interconversion rate constants related by k A / ¯ k A = exp ( − g A /k B T ) , and k ∗ n and ¯ k ∗ n are the rate con-stants for association and dissociation of inactive subunits tonucleated partial capsids. For simplicity, we take the nucleussize for inactive subunits binding equal to the energetic crit-ical nucleus size, and the association rate constant for inac-tive and active subunits to be equal above this size. Specifi-cally, for ‘Induced Fit’, k ∗ n = Θ( n − n nuc ) f n with Θ( n ) theHeaviside function and f n the association rate constant foractive subunits, while for ‘No Induced Fit’, k ∗ n = 0 ∀ n . Al-ternate choices for these quantities do not qualitatively affectthe results. Finally, dissociation rate constants are given bydetailed balance, ¯ k ∗ n = k ∗ n exp [(∆ g n − g A ) /k B T ] /v , wherethe free energy change upon association of an inactive sub-unit includes the activation energy g A . The initial condition is c ∗ (0) = (1 − f A ) c , c (0) = f A c , c n (0) = 0 for n > . C. Brownian dynamics model details
We use the patchy-sphere model presented in Ref. [68];our description here closely follows that reference. The min-imum energy structure is a complete capsid of 20 subunits,which have a spherical excluded volume with three attractivepatches, or bond vectors, that are separated by ◦ and rotaterigidly with the subunit. The attractive interaction betweentwo complementary bond vectors on respective subunits i and j is maximized when (1) the distance between the attractors r bij is minimized, (2) the angle θ bij between bond vectors isminimized, and (3) the dihedral angle φ bij calculated from twosecondary bond vectors, which are not involved in the primaryinteraction, is minimized. A schematic of the subunit interac-tions is shown in Fig. 2. Minimizing φ bij creates an interactionthat resists torsion and enforces angular specificity commen-surate with a complete capsid. The potentials are given by Eqs(16) [68, 89] U = U rep ( R ij ) + (cid:88) b U att ( r bij ) S ( θ bij , φ bij ) U rep ( R ij ) = L ( R ij , σ, σ ) U att ( r ij ) = (cid:15) b L (( r ij + 2 σ ) , σ, r c ) S ( θ, φ ) = 14 Θ( θ − θ c )Θ( φ − φ c )(cos( πθ/θ c ) + 1) (cos( πφ/φ c ) + 1) (16) with L p a generalized truncated and shifted Lennard-Jonesfunction: L p ( x, x c , σ ) ≡ (cid:18)(cid:16) xσ (cid:17) − p − (cid:16) xσ (cid:17) − p/ − (cid:16) x c σ (cid:17) − p + (cid:16) x c σ (cid:17) − p/ (cid:19) Θ( x − x c ) (17)In Eq. (16) the index b sums over pairs of complementarybond vectors, Θ( x ) is the Heaviside step function and R ij isthe subunit center-to-center distance. Conformational dynamics.
In our BD simulations withconformational dynamics, free subunits stochastically switchbetween inactive and active conformations. The only differ-ence between the two conformations is their interaction part-ners. For ‘No Induced Fit’, inactive subunits do not expe-rience attractive interactions with any subunit (but they stillexperience the repulsive excluded volume with all types ofsubunits), while pairs of active subunits experience the attrac-tive interactions described above. ‘Induced Fit’ has a simi-lar matrix of interactions, except that in active subunits ex-perience attractive interactions with any subunit in a partialcapsid which has at least one completed polygon (meaning aclosed cycle of interactions among subunits). In particular, in-active subunits experience no attractions with other inactivesubunits, free active subunits, or subunits in partial capsidswith no completed polygon, but do experience attractions topartial capsids with at least one completed polygon. Pairsof active subunits experience attractions as described above.Since we focus on the limit of fast conformational dynamics(relative to assembly timescales), subunits underwent confor-mational sampling with frequency . t − . At this frequency,the conformation of each free subunit was stochastically set toinactive or active with respective probabilities − f A and f A . Acknowledgements
This work was supported by the NIH, Award NumberR01GM108021 from the National Institute Of General Med-ical Sciences and the Brandeis Center for Bioinspired SoftMaterials, an NSF MRSEC, DMR-1420382. Computationalresources were provided by the NSF through XSEDE com-puting resources and the Brandeis HPCC which is partiallysupported by the Brandeis MRSEC.3 [1] S. Whitelam and R. L. Jack, Ann Rev Phys Chem , 143(2015).[2] M. F. Hagan, Adv. Chem. Phys. , 1 (2014).[3] A. Zlotnick, Virology , 269 (2003).[4] A. Rao, Annu. Rev. Phytopathol. , 61 (2006).[5] R. F. Garmann, A. Gopal, S. S. Athavale, C. M. Knobler, W. M.Gelbart, and S. C. Harvey, RNA , 877 (2015).[6] J. D. Perlmutter and M. F. Hagan, Annu. Rev. Phys. Chem. , 217 (2015), http://dx.doi.org/10.1146/annurev-physchem-040214-121637.[7] P. G. Stockley, N. A. Ranson, and R. Twarock, Future Virology , 531 (2013).[8] D. L. D. Caspar, T. Rosenstiel, and B. Medical, Biophys. J. ,103 (1980).[9] C. Packianathan, S. P. Katen, C. E. Dann, and A. Zlotnick, J.Virol. , 1607 (2010).[10] P. E. Prevelige, D. Thomas, and J. King, Biophys. J. , 824(1993).[11] A. Zlotnick, J. M. Johnson, P. W. Wingfield, S. J. Stahl, andD. Endres, Biochemistry , 14644 (1999).[12] A. Zlotnick, R. Aldrich, J. M. Johnson, P. Ceres, and M. J.Young, Virology , 450 (2000).[13] G. L. Casini, D. Graham, D. Heine, R. L. Garcea, and D. T.Wu, Virology , 320 (2004).[14] C. Chen, C. C. Kao, and B. Dragnea, J. Phys. Chem. A ,9405 (2008).[15] C. Berthet-Colominas, M. Cuillel, M. H. J. Koch, P. Vachette,and B. Jacrot, Eur. Biophys. J. with Biophys. Lett. , 159(1987).[16] S. Kler, R. Asor, C. Li, A. Ginsburg, D. Harries, A. Oppen-heim, A. Zlotnick, and U. Raviv, J. Am. Chem. Soc. , 8823(2012).[17] S. Kler, J. C.-Y. Wang, M. Dhason, A. Oppenheim, and A. Zlot-nick, ACS Chem. Biol. , 2753 (2013).[18] P. G. Stockley, O. Rolfsson, G. S. Thompson, G. Basnak,S. Francese, N. J. Stonehouse, S. W. Homans, and A. E.Ashcroft, J. Mol. Biol. , 541 (2007).[19] G. Basnak, V. L. Morton, O. Rolfsson, N. J. Stonehouse, A. E.Ashcroft, and P. G. Stockley, J. Mol. Biol. , 924 (2010).[20] C. Uetrecht, I. M. Barbu, G. K. Shoemaker, E. van Duijn, andJ. R. Heck Albert, Nat Chem , 126 (2011).[21] V. Baumgaertel, B. Mueller, and D. C. Lamb, Viruses , 777(2012).[22] N. Jouvenet, S. M. Simon, and P. D. Bieniasz, J. Mol. Biol. , 501 (2011).[23] A. Borodavka, R. Tuma, and P. G. Stockley, Proc. Natl. Acad.Sci. U. S. A. , 15769 (2012).[24] G. Tresset, C. Le Coeur, J.-F. Bryche, M. Tatou, M. Zeghal,A. Charpilienne, D. Poncet, D. Constantin, and S. Bressanelli,J. Am. Chem. Soc. , 15373 (2013).[25] E. E. Pierson, D. Z. Keifer, L. Selzer, L. S. Lee, N. C. Contino,J. C. Y. Wang, A. Zlotnick, and M. F. Jarrold, J. Am. Chem.Soc. , 3536 (2014).[26] K. Zhou, L. Li, Z. Tan, A. Zlotnick, and S. C. Jacobson, J. Am.Chem. Soc. , 1618 (2011).[27] Law-Hine, D. and Sahoo, A. K. and Bailleux, V. and Zeghal,M. and Prevost, S. and Maiti, P. K. and Bressanelli, S. and Con-stantin, D. and Tresset, G., J. Phys. Chem. Lett. , 3471 (2015).[28] P. Ceres and A. Zlotnick, Biochemistry , 11525 (2002).[29] A. Zlotnick, J. Z. Porterfield, and J. C.-Y. Wang, Biophys. J. , 1595 (2013). [30] A. Zlotnick, J. Mol. Biol. , 59 (1994).[31] D. Endres and A. Zlotnick, Biophys. J. , 1217 (2002).[32] P. Moisant, H. Neeman, and A. Zlotnick, Biophys. J. , 1350(2010).[33] P. van der Schoot and R. Zandi, Phys. Biol. , 296 (2007).[34] Q. Zhang, S. Gupta, T. Emrick, and T. Russell, J. Am. Chem.Soc. , 3898 (2006).[35] T. Keef, C. Micheletti, and R. Twarock, J. Theor. Biol. ,713 (2006).[36] M. Hemberg, S. Yaliraki, and M. Barahona, Biophys. J. ,3029 (2006).[37] E. C. Dykeman, P. G. Stockley, and R. Twarock, Phys Rev E , 022717 (2013).[38] G. R. Smith, L. Xie, B. Lee, and R. Schwartz, Biophys. J. ,310 (2014).[39] F. Jamalyaria, R. Rohlfs, and R. Schwartz, J. Comput. Phys. , 100 (2005).[40] B. Sweeney, T. Zhang, and R. Schwartz, Biophys. J. , 772(2008).[41] A. Zlotnick and S. Mukhopadhyay, Trends Microbiol. , 14(2011).[42] M. S. Kumar and R. Schwartz, Phys. Biol. , 045005 (2010).[43] L. Xie, G. Smith, X. Feng, and R. Schwartz, Biophys. J. ,1545 (2012).[44] S. Singh and A. Zlotnick, J. Biol. Chem. , 18249 (2003).[45] W. H. Roos, R. Bruinsma, and G. J. L. Wuite, Nat. Phys. , 733(2010).[46] C. Uetrecht, N. R. Watts, S. J. Stahl, P. T. Wingfield, A. C.Steven, and A. J. R. Heck, Phys. Chem. Chem. Phys. , 13368(2010).[47] L. Deshmukh, C. D. Schwieters, A. Grishaev, R. Ghirlando,J. L. Baber, and G. M. Clore, J. Am. Chem. Soc. , 16133(2013), pMID: 24066695, http://dx.doi.org/10.1021/ja406246z.[48] F. Birnbaum and M. Nassal, J. Virol. , 3319 (1990).[49] S. J. Stray, C. R. Bourne, S. Punna, W. G. Lewis, M. G. Finn,and A. Zlotnick, Proc. Natl. Acad. Sci. U. S. A. , 8138(2005).[50] P. T. Wingfield, S. J. Stahl, R. W. Williams, and A. C. Steven,Biochemistry , 4919 (1995).[51] S. Asakura, J. Mol. Biol. , 237 (1968).[52] Y. Kikuchi and J. King, J. Mol. Biol. , 695 (1975).[53] P. J. G. Butler, Philosophical Transactions of the Royal Societyof London Series B-Biological Sciences , 537 (1999).[54] A. Klug, Philos. Trans. R. Soc. Lond. B. Biol. Sci. , 531(1999).[55] D. L. Caspar and K. Namba, Adv. Biophys. , 157 (1990).[56] D. J. Kraft, W. K. Kegel, and P. van der Schoot, Biophys. J. , 2845 (2012).[57] D. L. Caspar, Curr. Biol. , 30 (1991).[58] D. L. D. Caspar and A. Klug, Cold Spring Harbor Symp. Quant.Biol. , 1 (1962).[59] B. Berger, P. W. Shor, L. Tuckerkellogg, and J. King, Proc.Natl. Acad. Sci. U. S. A. , 7732 (1994).[60] R. Schwartz, P. W. Shor, P. E. Prevelige, and B. Berger, Bio-phys. J. , 2626 (1998).[61] R. Schwartz, R. L. Garcea, and B. Berger, Virology , 461(2000).[62] O. M. Elrad and M. F. Hagan, Nano Lett. , 3850 (2008).[63] H. D. Nguyen, V. S. Reddy, and C. L. Brooks, Nano Lett. ,338 (2007). [64] H. D. Nguyen, V. S. Reddy, and C. L. Brooks, J. Am. Chem.Soc. , 2606 (2009).[65] V. L. Morton, E. C. Dykeman, N. J. Stonehouse, A. E. Ashcroft,R. Twarock, and P. G. Stockley, J. Mol. Biol. , 298 (2010).[66] M. F. Hagan, Biophys. J. , 1065 (2010).[67] R. Zandi, P. van der Schoot, D. Reguera, W. Kegel, andH. Reiss, Biophys. J. , 1939 (2006).[68] M. R. Perkett and M. F. Hagan, J. Chem. Phys. , 214101(2014).[69] A. Branka and D. Heyes, Phys. Rev. E , 2381 (1999).[70] D. Heyes and A. Branka, Mol. Phys. , 1949 (2000).[71] M. F. Hagan and D. Chandler, Biophys. J. , 42 (2006).[72] A. W. Wilber, J. P. K. Doye, A. A. Louis, E. G. Noya, M. A.Miller, and P. Wong, J. Chem. Phys. , 085106 (2007).[73] S. D. Hicks and C. L. Henley, Phys. Rev. E , 031912 (2006).[74] D. Rapaport, Phys. Rev. Lett. , 186101 (2008).[75] D. Endres, M. Miyahara, P. Moisant, and A. Zlotnick, ProteinSci. , 1518 (2005).[76] J. M. Johnson, J. H. Tang, Y. Nyame, D. Willits, M. J. Young,and A. Zlotnick, Nano Lett. , 765 (2005).[77] P. Ceres, S. J. Stray, and A. Zlotnick, J. Virol. , 9538 (2004).[78] A. Zlotnick, J. Mol. Biol. , 14 (2007).[79] W. K. Kegel and P. van der Schoot, Biophys. J. , 3905 (2004). [80] This assumption neglects the fact that rotational entropy penal-ties are most likely not proportional to the number of contacts(see Ref. [83] for further discussion).[81] L. Selzer, S. P. Katen, and A. Zlotnick, Biochemistry , 5496(2014), selzer, Lisa Katen, Sarah P. Zlotnick, Adam.[82] Since the master equations can be integrated to arbitrarily longtimes, we report the maximum value of P int observed at anytime for a given parameter set.[83] M. F. Hagan, J. Chem. Phys. , 114902 (2009).[84] M. F. Hagan, O. M. Elrad, and R. L. Jack, J. Chem. Phys. ,104115 (2011).[85] W. I. Sundquist and H.-G. Krausslich, Cold Spring Harbor per-spectives in medicine , a006924 (2012).[86] E. R. May, J. Feng, and C. L. Brooks, Biophys. J. , 606(2012).[87] D. Veesler and J. E. Johnson, Annual Review of Bio-physics , 473 (2012), http://dx.doi.org/10.1146/annurev-biophys-042910-155407.[88] A. Roldao, M. C. M. Mellado, J. C. Lima, M. J. T. Carrondo,P. M. Alves, and R. Oliveira, Plos Comp Biol , e1002367(2012).[89] M. F. Hagan, Phys. Rev. E77