{'fraction_non_alphanumeric': 0.034886905087901855, 'fraction_numerical': 0.014988042426565894, 'mean_word_length': 5.01559411377114, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Pathways for virus assembly around nucleic acids

Understanding the pathways by which viral capsid proteins assemble around their genomes could identify key intermediates as potential drug targets. In this work we use computer simulations to characterize assembly over a wide range of capsid protein-protein interaction strengths and solution ionic strengths. We find that assembly pathways can be categorized into two classes, in which intermediates are either predominantly ordered or disordered. Our results suggest that estimating the protein-protein and the protein-genome binding affinities may be sufficient to predict which pathway occurs. Furthermore, the calculated phase diagrams suggest that knowledge of the dominant assembly pathway and its relationship to control parameters could identify optimal strategies to thwart or redirect assembly to block infection. Finally, analysis of simulation trajectories suggests that the two classes of assembly pathways can be distinguished in single molecule fluorescence correlation spectroscopy or bulk time resolved small angle x-ray scattering experiments.


II. INTRODUCTION
In many virus families, the spontaneous assembly of a protein shell (capsid) around the viral genome is an essential step in the viral life cycle [43]. These families include most viruses with single stranded RNA (ssRNA) genomes, as well as the Hepadnaviridae (e.g. hepatitis B virus, HBV). Understanding the mechanisms which underlie this cooperative assembly process could facilitate efforts to develop anti-viral drugs that block or derail the formation of infectious particles (for reviews see [91,128]) and promote efforts to reengineer them for biomedical delivery. In this article, we explore how the interactions between the molecular components determine the mechanism of assembly, and how these interactions can be altered by changing solution conditions or mutagenesis to modulate assembly pathways.
The most detailed knowledge of capsid-RNA interactions comes from structural analysis of assembled viral particles. Atomic resolution structures of capsids assembled around ssRNA have been obtained by x-ray crystallography and/or cryo-electron microscopy (cryo-EM) (e.g. [34,38,45,57,65,74,98,106,112,[114][115][116][117]). The packaged NAs are less ordered than their protein containers and thus have been more difficult to characterize. However cryo-EM experiments have identified that the nucleotide densities are nonuniform, with a peak near the inner capsid surface and relatively low densities in the interior [21,112,127]. While atomistic detail has not been possible in these experiments, all-atom models have been derived from equilibrium simulations [24,35,125]. In some cases, striking image reconstructions reveal that the packaged RNA adopts the symmetry of the overlying capsid (e.g. [4,32,98,112,114]). While it has been proposed that this order arises as a function of the assembly mechanism for several viruses [26,27,67], computational analysis of polyelectrolyte configurations inside capsids also indicate that capsid-polymer interactions can generically drive spatial organization of the packaged polymer [2,3,7,24,29,33,44,52,56,69,89,101,126]. Theoretical works have also characterized the relationship between the NA charge and structure and the length which is optimal for packaging [2, 7, 50, 83, 89, 99, 101-103, 113, 118].
In addition to this structural data on assembled capsids, an extensive combination of mass spectrometry, assembly kinetics experiments, constraints from assembled capsid structures, and mathematical modeling has delineated assembly pathways for several viruses, with a particular focus on the role of interactions between CPs and specific RNA sequences called packaging signals. Recent single molecule fluorescence correlation spectroscopy (smFCS) experiments indicate that, for these viruses, assembly around the viral genome is more robust and proceeds by a different mechanism as compared to around heterologous RNA [8]. Yet, in other cases capsid proteins show no preference for genomic RNA over heterologous RNA (e.g. HBV [90]), and cowpea chlorotic mottle virus (CCMV) proteins preferentially encapsidate heterologous RNA (from BMV) over the genomic CCMV RNA with equivalent length [19]. Furthermore, experimental model systems in which capsid proteins assemble into icosahedral capsids around synthetic polyelectrolytes or other polyanions [6,16,20,22,25,37,47,49,53,72,73,104,108]) demonstrate that specific RNA sequences are not required for capsid formation or cargo packaging. Thus, a complete picture of capsid assembly mechanisms requires understanding how assembly pathways depend on those features which are generic to polyelectrolytes, as well as those which are specific to viral RNAs.
In previous work on assembly around a simple model for a polymer, Elrad and Hagan [29] proposed that mechanisms for assembly around a cargo (i.e. RNA, polymer, or nanoparticle) can be classified on the basis of two extreme limits. In the first (originally proposed by McPherson [78] and then Refs [24,40]), strong protein-cargo interactions drive proteins to adsorb 'en masse' onto the cargo in a disordered manner, meaning there are few protein-protein interactions. Once enough subunits are bound, subunits undergo cooperative rearrangements (potentially including dissociation of excess subunits) to form an ordered capsid. This mechanism has been observed in recent simulations [29,40,75,88,89]. In the second limit, where protein-protein interactions dominate, a small partial capsid nucleates on the cargo, followed by a growth phase in which individual proteins or small oligomers sequentially add to the growing capsid. This class of pathways resembles the nucleation-and-growth mechanism by which empty capsids assemble [30], except that the polymer plays an active role by stabilizing protein-protein interactions and by enhancing the flux of proteins to the assembling capsid [29,51,62].
It is difficult to determine assembly mechanisms directly from experiments, due to the small size ( 10 of nm) and transience (∼ms) of most intermediates. Observations in vitro suggest that both mechanisms may be viable. Kler et al. [63] used time-resolved X-ray scattering (trSAXS) to monitor SV40 capsid proteins assembling around ssRNA. The scattering profiles at all time points during assembly could be decomposed into unassembled components (RNA + protein subunits) and complete capsid; the absence of any signal corresponding to a large disordered intermediate suggests this assembly follows the nucleation-and-growth (ordered) assembly mechanism [63]. Other observations suggest that viruses can assemble through the en masse mechanism. Refs. [13,36] found that in vitro assembly CCMV assembly was most productive when performed in two steps: (1) at low salt (strong protein-RNA interactions) and neutral pH (weak protein-protein interactions) the proteins undergo extensive adsorption onto RNA, then (2) pH is reduced to activate binding of protein-protein binding [36]. Similarly, a recent observation of capsid protein assembly around charge-functionalized nanoparticles found that assembly initially proceeded through nonspecific aggregation of proteins and nanoparticles, followed by the gradual extrusion of nanoparticles within completed capsids [76].
These experiments used viral proteins with relatively weak protein-protein interactions (CCMV and BMV) [129] and moderate salt concentrations (100 − 150mM). The Kler et al. [63,64] experiments considered SV40 proteins, which have strong protein-protein interactions [129], and high salt (250mM) Together, these in vitro experiments suggest that productive assembly could proceed by either the en masse or nucleation-and-growth mechanism.
In this work, we use dynamical simulations to investigate the extent to which the assembly mechanism can be controlled by tuning solution ionic strength and protein-protein attractions. We extend a model that was recently used to calculate the thermostability and assembly yields of viral particles as a function of protein charge and nucleic acid (NA) length and structure. Those previous simulations found quantitative agreement between predicted NA lengths that optimize capsid thermostability and viral genome length for seven viruses [89]. Here, we perform extensive new simulations of assembly, in which protein-protein interactions, the sequence of charges in capsid protein-NA binding domains, and the solution ionic strength are varied. We find that, by varying these control parameters, the assembly mechanism can be systematically varied between the two extreme limits described above. Our results suggest that knowledge of protein-protein and protein-NA binding affinities may be sufficient to predict which assembly mechanism will occur, and we estimate relative protein-NA binding interactions for three viruses (based on non-specific electrostatic interactions). These findings suggest that assembly mechanisms can be rationally designed through choice of solution conditions and mutagenesis of capsid protein-protein interfaces and protein-NA binding domains. Finally, by calculating hydrodynamic radii and SAXS profiles associated with assembly intermediates, we show that assembly mechanisms can be distinguished by experimental techniques recently applied to virus assembly, respectively smFCS [8]) and trSAXS [63]. While the NA is represented by a linear polyelectrolyte in most of the simulations, we obtain qualitatively similar results when considering a model for base-paired NAs developed in Ref. [89].

III. RESULTS
To study how capsid assembly around a polyelectrolyte depends on the strength of protein subunit-subunit and subunit-polyelectrolyte interactions, we performed Brownian dynamics simulations with a recently developed model [89] (Figure 1). This model was motivated by the observation [63,64] that purified simian virus 40 (SV40) capsid proteins assemble around ssRNA molecules to form capsids composed of 12 homopentamer subunits. The capsid is modeled as a dodecahedron composed of 12 pentagonal subunits (each of which represents a rapidly forming and stable pentameric intermediate, which then more slowly assembles into the complete capsid, as is the case for SV40 in vitro around ssRNA [63,64,71]). Subunits are attracted to each other via attractive pseudoatoms at the vertices (type 'A') and driven toward a preferred subunit-subunit angle by repulsive 'Top' pseudoatoms (type 'T') and 'Bottom' pseudoatoms (type 'B') (see Fig. 1 and the Methods section). These attractions represent the interactions between capsid protein subunits that arise from hydrophobic, van der Waals, and electrostatic interactions as well as hydrogen bonding [43]). The magnitude of these interactions varies between viral species [129] and can be experimentally tuned by pH and salt concentration [15,43,61]. Here, the attraction strength is controlled by the model parameter ε ss . In order to relate this potential energy to the free energy of dimerization, we have run a separate series of calculations, where in the absence of cargo, we find that the free energy of subunit of dimerization is g ss =5.0 − 1.5 * ε ss (see SI section 12). Throughout this article energies are given in units of the thermal energy, k B T .
Capsid assembly around nucleic acids and other polyelectrolytes is primarily driven by electrostatic interactions between negative charges on the encapsulated polyelectrolyte and positive charges on the inner surfaces of capsid proteins [43]. We consider a linear bead-spring polyelectrolyte, with a charge of e per bead and a persistence length comparable to that of ssRNA in the absence of base pairing. Positive charges on capsid proteins are located in flexible polymers affixed to the inner surface of the model subunit, which represent the highly charged, flexible terminal tails known as arginine rich motifs that are typical of positive-sense ssRNA virus capsid proteins (e.g., [98]). Electrostatics are modeled using Debye-Huckel interactions, where the Debye screening length (λ D ) is determined by the ionic strength I, equivalent to the concentration C salt of a monovalent salt, as λ D ≈ 0.3/I 1/2 with λ D in nm and I in molar units. Perlmutter et al. [89] showed that Debye-Huckel interactions compare well to simulations with explicit counterions for the parameter values under consideration; furthermore, inclusion of multivalent ions at experimentally relevant concentrations had only a small effect on properties such as the optimal genome length. The effect of base-pairing on the optimal length was also considered in that reference [89].

A. Kinetic phase diagram
We first consider the predominant assembly products (Fig. 2) and the packaging efficiency (Fig. 3) as a function of Debye length λ D and subunit-subunit interaction strength ε ss . The packaging efficiency is defined as the fraction of trajectories in which the polyelectrolyte is completely encapsulated by a well-formed capsid, which contains 12  [105]. Space-filling model of the generic subunit model (middle) and a pentamer from the PC2 model (right). Beads are colored as follows: blue=excluders, green=attractors, yellow=positive ARM bead, gray=neutral ARM bead, red=polyelectrolyte.
subunits each of which interact with 5 neighbors. We refer to this as a kinetic phase diagram [41,43] because we characterize products at a finite observation time of t f =2 × 10 8 time steps, which is long enough that assemblages do not vary significantly with time (except for under weak interactions, see below), but is not sufficient to equilibrate kinetic traps if there are large activation barriers [29,43]. We see that for the range of simulated ionic strengths (1 ≤ C salt ≤ 500 mM or 10 ≥ λ D ≥ 0.4 nm) assembly yields are highest for ε ss =5k B T and C salt =100 mM (the parameter values focused on in [89]), and that for moderate subunit-subunit interaction strengths (4 ≤ ε ss ≤ 6, k B T ) yield to remain high as the ionic strength is increased to about 300 mM (λ D ≈ 0.6 nm). For higher salt concentrations, yields are depressed by the appearance of long-lived on-pathway intermediates. As will be discussed further below, weakening the electrostatic interaction between the polymer and protein limits the ability of the polymer to promote assembly. Although we expect that these simulations would eventually result in complete capsids, the low yield at our finite measurement time reflects the fact that assembly is less efficient than for lower salt concentrations. At the highest ionic strength considered (C salt =500 mM), the most prevalent outcome is that no nucleation occurs. At lower salt concentrations (C salt ≥ 10 mM), rapid adsorption of a super-stroichiometric quantity of subunits results in malformed capsids. At larger-than-optimal protein-protein interaction strengths (ε ss > 6) assembly yields are lower for two reasons. The first, and more frequent outcome, is long-lived malformed structures with strained interactions. This kinetic trap arises in a wide variety of assembly systems when interactions become strong in comparison to the thermal energy, because strained interactions are unable to anneal before additional subunits lock them in place [39,42,43]. In our simulations, we found that these structures frequently result from an incorrect merger between two partial capsid intermediates; even when each individual interaction is too weak to lock-in the non-optimal structure, multiple erroneous interactions formed by two partial capsids are effectively irreversible on our timescale. The tendency for oligomer binding to lead to malformed structures was seen previously in the context of empty capsid assembly [41,120]. Here, the polymer helps to bring oligomers together, and thus this trap arises when nucleation on the polymer is faster than growth of a nucleus into a complete capsid. This trap resembles Geminivirus particles, which are composed of a union of two nearly complete capsids [9].
The second obstacle to polyelectrolyte encapsulation arises at the highest protein-protein interaction strengths studied (ε ss ≥ 8k B T ), for which subunits not associated with the polyelectrolyte undergo spontaneous assembly. The resulting off-polyelectrolyte assembly depletes the pool of available monomers and small oligomers available for assembly on the polyelectrolyte, leading to a form of the monomer-starvation kinetic trap previously discussed for empty capsid assembly [30,43]. Triggering formation of empty capsids and thus preventing nucleic acid encapsidation by strengthening subunit-subunit interactions has been suggested as a mode of action for a putative antiviral drug for hepatitis B virus [59,60].
At smaller-than-optimal protein-protein interaction strengths (ε ss < 4) assembly is unsuccessful for two reasons, depending upon the ionic strength. At smaller values (C salt ≤ 300 mM), electrostatic interactions are relatively strong, and many proteins adsorb onto to the polymer. However, because of the weak protein-protein interaction strength, these proteins do not form stable capsids, predominantly because nucleation is slow in comparison to t f (the final observation time, t f =2 × 10 8 time steps). In a minority of cases, a nucleus will form, but completion is prevented by the excess number of subunits adsorbed to the polyelectrolyte. We refer to the resulting configurations as disordered, due to the lack of ordered binding between protein subunits. At larger ionic strengths (C salt > 300 mM) electrostatic interactions are relatively weak, and individual subunits rapidly desorb from the polyelectrolyte. In this regime, assembly requires a fluctuation in the number of adsorbed subunits which leads to nucleation of a partial capsid intermediate which has enough subunit-polyelectrolyte interactions to be stable against rapid desorption. The nucleation rate decreases exponentially with subunit-subunit interaction strength [43,62], and thus most simulations at high salt and weak subunit-subunit interactions never undergo nucleation. (We categorize simulations with fewer than 3 subunits adsorbed to the polyelectrolyte and no subunit-subunit interactions (i.e. no progress towards assembly) as 'unnucleated'.) Importantly, we expect that trajectories in this region of parameter space will eventually undergo nucleation. Thus, as the finite assembly time t f is increased the region of successful assembly will expand to lower values of ε ss and higher ionic strength, until eventually reaching values below which capsid assembly is thermodynamically unstable (see [41], Fig. 7). To confirm this possibility, we used the Markov state model (MSM) approach described in Ref. [88] to characterize assembly with ε ss =5, 6 and C salt =500 mM (see Figs. 4C and S4 (in the SI)). In contrast, the malformed capsids encountered under large parameter values typically will not resolve on any relevant timescales, since numerous strong interactions would need to be broken [43]. Thus, the boundaries of successful assembly at large ε ss values are insensitive to t f . Several additional conclusions can be drawn from the variation of packaging efficiency (Fig. 3). Firstly, the yield of complete virus-like particles is relatively high for moderate subunit-subunit interaction strengths ε ss ∈ [4, 6] across a range of ionic strengths (C salt ∈ [100, 300] mM). Even above (C salt > 300 mM) and below (C salt =10 mM) this range we observe moderate yields of complete particles. Secondly, as electrostatic interactions are weakened (moving to the left), the subunit-subunit interaction strength which optimizes the yield increases (i.e. from ε ss = 5k B T at C salt = 100mM to ε ss = 6k B T at C salt : 150 − 400mM to ε ss = 7k B T at C salt = 500mM ). This result suggests that one interaction type can compensate for the other within a limited range. However, though all successful capsids contain the same subunit geometry, the mechanism by which they form depends on the relative interaction strengths, as discussed next.

B. Assembly mechanisms
As noted previously [29,40,43], pathways for assembly around a central core such as a polyelectrolyte can be roughly separated into two classes. In the first class ( Fig. 4a), which we refer to as the 'en masse' mechanism, subunits first adsorb onto the polyelectrolyte in a disordered manner, followed by cooperative rearrangements to form an ordered capsid. In the second class ( Fig. 4b), referred to as the nucleation-and-growth mechanism, a small partial capsid nucleates on the polyelectrolyte followed by the sequential, reversible addition of subunits or small oligomers until assembly completes. In contrast to the earlier models which considered a qualitative subunit-polyelectrolyte interaction, we study here how assembly pathways depend on the ionic strength.
To quantify the degree of order along assembly pathways, we record the total number of subunits adsorbed to the polyelectrolyte n ad and the number of subunits in the largest cluster n. Trajectories that pass through configurations with a large value of n free = n ad − n are disordered, with many adsorbed subunits not participating in ordered assemblies. In Fig. 4, these quantities are shown as a function of time averaged over all simulation trajectories (leading to successful assembly or not), for parameter sets that respectively lead to the en masse mechanism (Fig. 4A) and nucleation-and-growth mechanism (Fig. 4B). In the first case, there are strong subunit-polyelectrolyte interactions (low ionic strength, C salt =100 mM) and weak subunit-subunit interactions (ε ss =3k B T ). Subunits therefore initially adsorb nonspecifically and form only transient subunit-subunit interactions, leading to a rapid rise in n ad with n ≈ 0. Once enough subunits are adsorbed (∼ 12 around an optimal-length polyelectrolyte for this model with a 12-subunit capsid), a cooperative fluctuation in subunit configurations eventually leads to a stable nucleus and then rapid completion of the ordered capsid geometry. Since this nucleation process is stochastic, there is a distribution of waiting times and thus a more gradual increase in the average cluster size n (see Fig. S2). In the nucleation-and-growth case, on the other hand, the subunit-polyelectrolyte interactions are weak (C salt =300 mM) and the subunit-subunit interactions are strong (ε ss =6k B T ). There is limited nonspecific subunit adsorption onto the polyelectrolyte, adsorbed subunits form relatively strong associations, and thus n ad and n increase nearly in lockstep. Snapshots from typical trajectories for each of these parameter sets are shown in Fig. 4A.
To visualize the degree of order as a function of parameter values, we define a trajectory-averaged order parameter n free , which is n free averaged over all configurations with 4 ≤ n ≤ 6 and over all trajectories at a given parameter set [28]. Large values of this parameter (n free 5) indicate the en masse mechanism, while small values (n free 2) indicate the nucleation-and-growth mechanism. As shown in Fig. 5, the degree of order generally increases with ionic strength and subunit-subunit interaction strength, with the most ordered assembly occurring at C salt = 500 mM (where on average fewer than one subunit is adsorbed nonspecifically) and ε ss ≥ 6k B T . However, notice that for ε ss = 3 assembly is always disordered; for such weak subunit-subunit interactions the critical nucleus size is large and a high density of adsorbed subunits is required to achieve nucleation. On the other hand, for moderate subunit-subunit interactions we do not observe the extreme limit of the en masse mechanism even for low ionic strength. Though a low  ionic strength drives strong nonspecific subunit adsorption, absorbed subunits collide frequently due to cooperative polymer motions and subunit sliding along the polymer [29,51]. For ε ss > 3 absorbed subunits therefore achieve nucleation before nonspecific absorption has time to saturate. The nucleation-and-growth trajectories can be further classified based on the relative timescales of nucleation and growth. When nucleation is slow compared to growth (Figs. 4C and S4), the reaction is effectively two-state -each critical nucleus rapidly proceeds to completion, leading to low concentrations of intermediates which are essentially undetectable in a bulk experiment. When nucleation and growth timescales are comparable, multiple capsids within a spatial region can be assembling simultaneously, and thus potentially could be detected in bulk experiments. Below, we consider whether it is possible to experimentally distinguish between the latter case, ordered assembly with rapid nucleation, and en masse assembly pathways characterized by disordered intermediates.
Biological ARM sequences. The assembly simulations described above consider a simplified peptide ARM, which contains only 5 positive ARM charges (Fig. 1). However, our previous work demonstrated that the ARM amino acid sequence (specifically the charges) can significantly affect binding to a polyelectrolyte [89]. Furthermore, results from a more simplified model suggested that the polyelectrolyte binding affinity is a determining factor for the assembly mechanism [29]. To test this hypothesis and to place our results in the context of specific viruses, we calculated the polyelectrolyte binding affinity for the simplified ARM and three ARM sequences from biological viruses (SV40, BMV, PC2), with each amino acid classified as neutral, cationic, or anionic. These calculations were performed by setting the subunit-subunit interaction strength to zero and recording the average number of adsorbed  Figure 5. Dependence of the assembly mechanism on the subunit-subunit attraction strength εss and the ionic strength C salt . The assembly order parametern free (the average number of adsorbed subunits not in the largest partial capsid) is shown as a function of parameter values. As described in the text,n free 2 corresponds to ordered nucleation-and-growth assembly pathways, whereas larger values correspond to disordered pathways. The alternate y-axis shows the subunit-subunit dimerization free energy gss corresponding to εss (see section VII B) and the alternate x-axis shows the linear adsorption density of subunits in the absence of assembly (see Fig. 6 below).
subunits at varying C salt (or Debye length λ D ). The measured equilibrium linear densities of adsorbed subunits, c eq , are shown as a function of C salt in Fig. 6A-C. In all cases the linear density increases monotonically with decreasing C salt ; saturating at a maximum density. The simplified ARM has the largest binding affinity, despite having the lowest net positive charge (+5, 0 neutral segments) of the four species considered. Comparison with SV40 (+6, 14 neutral segments) and BMV (+9, 33 neutral segments) illustrates the expected result that that neutral segments decrease the binding affinity, particularly at high salt. The PC2 subunits, with net charge +22 demonstrate markedly different behavior, with significant subunit absorption at the highest salt concentration simulated (500 mM), but saturating at about 300 mM ionic strength due to subunit-subunit charge repulsions. Variations in adsorption density with subunit concentration are shown in Fig. S1.
In Fig. 6D we connect these measures of adsorption with assembly mechanism by plotting the assembly order parametern free as a function of c eq (controlled by varying C salt ), for several values of the subunit dimerization free energy (determined by ε ss and C salt ) for the 5-ARM and SV40 models. Plotting against c eq rather than C salt (λ D ) allows us to overlay data from these two models while accounting for the differences in affinity due to ARM sequence described above. In support of the proposed link between binding affinity and assembly mechanism, we find rough agreement in the measured assembly order parameters between the two models. The results indicate that ARM sequence can significantly influence the assembly mechanism. For example, if we definen free ≤ 2 as the nucleation-growth mechanism, Fig. 6D indicates that nucleation-growth occurs for c eq ≤ c * eq with the threshold value c * eq = 0.0375 for ε ss ≥ 2.5k B T . From Fig. 6A-C, we can then identify the threshold values of ionic strength C * salt , above which the nucleation-growth mechanism will occur: C * salt ≈ 300 mM for the 5-ARM and C * salt ≈ 175 mM for the BMV and SV40 models; while PC2 is below the threshold value for all salt concentrations considered. This allows us to predict, for example, that recent experiments on SV40 assembly (at C salt =250mM and observed strong subunit-subunit attraction) would have a very lown free (∼ 1), which is consistent with SAXS observations [63].

C. Experimental Observables
We now seek to provide experimental signatures through which the capsid assembly pathways discussed above can be distinguished. We focus on two experimental techniques which have recently been applied to respectively probe the formation of individual capsids and bulk assembly kinetics.
smFCS measurements on individual capsids can distinguish assembly mechanisms. Borodavka et al. nucleocapsid complexes during assembly of MS2 and STNV capsid proteins around cognate RNA or non-cognate RNA. Assembly around cognate RNA was characterized by either constant R H or, in some trajectories, a collapsed complex followed by gradual increase until reaching the size of an assembled capsid. In contrast, assembly around non-cognate RNA led to an increase in R H before eventually decreasing to the size of the capsid. The difference between these two behaviors can be attributed to sequence-specific 'packaging signals' on cognate RNA that interact with the capsid proteins. In this article we do not consider the effect of packaging signals (these will be considered in a subsequent article); instead, we consider whether the pathways described in the previous section can be distinguished by this experimental technique.
We estimated the hydrodynamic radii R H for polymer-subunit intermediates using the program HYDROPRO, which has been shown to accurately predict R H for large protein complexes [85]. The resulting R H are shown during representative en masse and ordered assembly trajectories in Fig. 7E. We see that the complex R H first increases as subunits adsorb onto the complex and then decreases as subunits assemble (Fig. 7E). However, the en masse mechanism leads to a much larger and longer duration increase in R H , due to the extensive and long-lived disordered adsorption of unassembled subunits. The difference in R H between en masse and ordered trajectories is conserved across their respective parameter ranges (see SI Fig. S3 for other trajectories), and also occurs for assembly trajectories around the model nucleic acid with intramolecular base-pairing developed in Ref. [89] (see SI Fig. S3C,F). These results suggest that smFCS can distinguish among the classes of assembly pathways observed in our simulations. They are consistent with an interpretation of the Borodavka et al. [8] results in which assembly around the non-cognate RNA proceeds via the disordered mechanism while packaging signals lead to an ordered mechanism.
To further characterize the differences in polymer conformations between disordered and ordered assembly pathways, we show the polymer radius of gyration, R g , during assembly trajectories in Fig. 7. In contrast to R H , contributions of the capsid proteins are not included in R g . While the results in Fig. 7 are averaged over multiple trajectories, example individual trajectories are shown in SI Fig. S2. In all cases of successful assembly, the polymer is gradually compacted from its free size into its encapsidated size. However, at fixed C salt =100mM, the average rate of compaction increases with ε ss , with a dramatic increase in rate for ε ss > 3k B T (Fig. 7A). Similarly, decreasing C salt increases the rate of compaction (Fig. 7B). Notice that the rate of polymer compaction is not determined by the assembly mechanismincreased order correlates with faster compaction in Fig. 7A but with slower compaction in Figure 7B. When the R g is plotted as a function of number of adsorbed subunits (n ad ), the en masse and ordered pathways clearly split into two groups (Fig. 7C). However, this distinction disappears when R g is plotted as a function of the number of subunits in the largest cluster (n, Fig. 7D). Taken together, these data demonstrate that polymer compaction is driven by adsorbed subunits forming ordered intermediates, with the rate of compaction consequently mirroring the rate of capsid assembly.
trSAXS measurments of bulk assembly kinetics can distinguish assembly mechanisms. While smFCS can detect individual assembly intermediates, Kler et al. [63,64] recently used time resolved small angle x-ray scattering (trSAXS) to elucidate structures of assembling capsids in bulk. They found that the SAXS profiles at all time points could be decomposed into scattering contributions from the unassembled and complete components, suggesting that assembly proceeded by an effectively two-state reaction with undetectable levels of intermediates (i.e. the nucleation-and-growth pathway with relatively slow nucleation). While it is evident that profiles from profiles from a two-state reaction can be decomposed in this way, we investigated the extent to which SAXS profiles from the other pathway classes (en masse or nucleatation-and-growth with rapid nucleation) can be distinguished from the two-state case.
First, Fig. 8A,B shows SAXS profiles calculated (using CRYSOL [109]) from simulation snapshots along ensembles of en masse and ordered assembly trajectories. For each parameter set, SAXS profiles are averaged over 6 time windows. In both cases, the first profile is dominated by scattering from free subunits (Fig. 8C) and the final time segment shows clear minima and maxima corresponding to the complete capsid (Fig. 8C). For comparison, Fig. 8C presents the SAXS profiles for ordered subunit clusters ranging in size from a single subunit (black) to a complete  capsid (yellow). As the capsid grows distinct minima and maxima appear and become more pronounced. We note that the positions of the minima and maxima in the complete model capsid are similar to those observed experimentally for SV40 [63].
To test the extent to which these trajectories can be distinguished from the two-state case, we attempted to fit SAXS profiles using a linear combination of the profiles for unassembled polymer and subunits and the complete capsid. The resulting fits for the second fifth of the trajectory (where intermediates are most plentiful) are shown in Figs. 8D-F. At this stage, the ordered systems contain mostly incomplete subunit clusters (6-11 subunits), while the disordered simulations contain mostly adsorbed but unassembled subunits. For the ordered simulations, we find that the fit reproduces all of the features of the intermediates, except at low q. In contrast, the intermediates in the disordered trajectory display a shoulder at q ∼ 0.3nm −1 that is not observed in any of the unassembled or assembled components. This shoulder is a distinct feature of the disordered intermediate, and thus could be used to identify this class of pathways. Finally, as expected, SAXS profiles from trajectories at parameter sets (Fig. 4C) which lead to two-state kinetics are very well fit by a linear combination of polymer/subunit and complete capsid (Fig. 8F).
A significant distinction between the SAXS profiles is that the ordered pathways lead to an isosbestic point at q ∼ 0.3nm −1 (as seen in SV40 experiments [63]), whereas the disordered pathways do not. The presence of an isosbestic point is frequently assumed to indicate two-state behaviour; however, it occurs in our data for ordered trajectories even when the reaction is far from two-state due to rapid nucleation. In these cases the isosbestic point appears due to due to the similarity in scattering from the ordered intermediates and the complete capsid. This suggests that an isosbestic point may distinguish ordered from disordered assemblies, but is less sensitive to the extent to which the reaction kinetics can be approximated as two-state (i.e., how undetectable the intermediate concentrations are).

IV. DISCUSSION
Our simulations demonstrate that capsid assembly around a NA or other cargo can proceed through two mechanisms, and that which mechanism occurs depends on the relative strengths of protein-protein and protein-cargo interactions. The assembly mechanism can be rationally designed by tuning these interactions through solution conditions or mutagenesis of capsid protein-NA binding domains. However, because productive assembly requires weak interactions, the parameters must be tuned within relatively narrow ranges, and alterations which increase (decrease) the strength of one type of interaction must be compensated by a decrease (increase) in the other. Our results suggest that the subunit-cargo dissociation constant is an important parameter, whose value might be used to map specific viruses onto our phase diagrams (Figs. 2, 3 ,5), although experimental tests of this capability are required. Finally, we have shown that the nature of assembly pathways can be inferred from the results of recently developed techniques to monitor the assembly of individual capsids or bulk assembly kinetics.
Our simulations predict that a single viral species can be driven to assemble via different mechanisms in vitro by changing solution conditions. In particular, under a constant subunit-subunit binding energy (ε ss ∼ 5 − 6k B T ) robust assembly occurs for a range of solvent conditions, with highly ordered or disordered assemblies occurring depending on salt concentration (C salt ∈ [50, 400] mM). To our knowledge, this prediction has not yet been realized experimentally, although the signatures of the two classes of assembly pathways have been seen in experiments on different viruses and/or different cargoes [8,13,36,63,76].
One recent experimental study sought to test the role of of C salt and subunit-subunit attractions (controlled by solution pH) on in vitro assembly of CCMV around RNA [36]. In some regards, these experiments mirror our observations, with malformed capsids observed for strong subunit-subunit attractions, disordered products for weak subunit-subunit attractions, and well-formed capsids at intermediate conditions. However, robust assembly was only observed for a two-step process: first C salt was reduced (from 1M ) to drive RNA-subunit interactions, and then secondly pH was reduced to activate subunit-subunit attractions. The resulting assembly pathways resemble the en masse mechanism described here. On the other hand, one-step assembly led to malformed particles even at moderate pH, suggesting an inability to assemble through a nucleation-and-growth mechanism in these experiments. In our simulations, this outcome would only occur at high salt concentrations (e.g. C salt ∼ 400 − 500mM, see Figs. 2, 3), where the narrow range of ε ss leading to successful assembly indicates that parameters must be finely tuned. Reproducing such a lack of successful assembly at moderate salt concentrations would require a reduction in the orientation-dependence of the subunit-subunit interaction potential (see Methods), or introduction of additional factors as discussed below. Experiments in which solution conditions are changed for other viruses which do undergo one step assembly (e.g. SV40 [63,64]), may elucidate which of these possibilities to pursue and would directly test our prediction that the assembly mechanism can be controlled by solution conditions. Understanding capsid assembly mechanisms and their location within the assembly phase diagram has important implications for the design of antiviral agents. As one example, we consider the recently developed class of HBV inhibitors based on phenylpropenamides [59,60], which act by increasing the strength of subunit-subunit interactions, driving subunits to assemble in the absence of their genome and thus increasing the generation of empty, non-infective capsids [59,60]. Comparing Figs. 3 and 5 shows that a virus that undergoes ordered assembly (e.g. ε ss ∼ 6 − 7, C salt ∼ 300 mM), sits close to parameters that support empty capsid assembly, which are demarcated by dashed lines in Fig. 3. Thus only a small increase in subunit-subunit interactions is required to trigger unproductive RNA-free assembly. In contrast, a much larger perturbation would be required to achieve empty capsid assembly for a virus that assembles via the en masse mechanism.
Outlook. We have described two classes of assembly pathways and several distinct failure modes (Fig. 2) that arise when assembly is driven by nonspecific electrostatic subunit-cargo interactions. Our phase diagrams can serve as a starting point to understand how virus-specific features, such as packaging signals [107], allosteric NA-induced [64,106] or 'autosteric' protein-induced conformational changes [14], base pairing-induced NA-structure [8,89,124], or subcellular localization [5] can robustly avoid failure modes amidst the crowded intracellular milieu while enabling selective assembly around the viral genome [32,93,97]. For example, allosteric or autosteric conformational changes may allow for strong subunit-subunit interactions on the NA while avoiding the off-cargo assembly we observe at large ε ss . Systematically studying how these additional factors expand (or contract) regions of parameter space that lead to successful assembly will ultimately reveal how viruses have optimized their structures and interactions for robust assembly in vivo and how their assembly in vivo or in vitro can be manipulated for biomedical or materials science applications.

Model.
We have recently a presented a complete description of our model system, which we summarize here briefly [89] and in the SI. Our model subunits are based upon that previous used to simulate assembly of empty capsids [31,58,119] which we extended previously to model assembly around cargo [89]. The pseudoatoms in the capsid subunit model are illustrated in Fig. 1. Subunit assembly is mediated through an attractive Morse potential between Attractor ('A') pseudoatoms located at each subunit vertex. The Top ('T') pseudoatoms interact with other 'T' psuedoatoms through a potential consisting of the repulsive term of the Lennard-Jones (LJ) potential, the radius of which is chosen to favor a subunit-subunit angle consistent with a dodecahedron (116 degrees). The Bottom ('B') pseudoatom has a repulsive LJ interaction with 'T' pseudoatoms, intended to prevent 'upside-down' assembly. The 'T', 'B', and 'A' pseudoatoms form a rigid body [31,58,119]. See Refs. [29, 40-42, 46, 58, 70, 75, 79-81, 94-96, 100, 121-123] for related models.
To model electrostatic interaction with a negatively charged NA or polyelectrolyte we extend the model as follows. Firstly, to better represent the capsid shell we add a layer of 'Excluder' pseudoatoms which have a repulsive LJ interaction with the polyelectrolyte and the ARMs. Each ARM is modeled as a bead-spring polymer, with one bead per amino acid. The 'Excluders' and first ARM segment are part of the subunit rigid body. ARM beads interact through repulsive LJ interactions and, if charged, electrostatic interactions modelled by a Debye-Huckel potential. We note that repulsion between subunits due to the positive charges does affect the magnitude of the subunit-subunit interaction. Previously, we estimated this repulsion contributes 1k B T to the dimerization free energy [89]. This contribution (and entropic terms) are not included in ε ss , which is the magnitude of the Morse potential depth. See SI section VII B) for a discussion of binding free energies.
Simulations and units. Simulations were performed with the Brownian Dynamics algorithm of HOOMD, which uses the Langevin equation to evolve positions and rigid body orientations in time [1,68,82]. Simulations were run using a set of fundamental units. The fundamental energy unit is selected to be E u ≡ 1k B T . The unit of length D u is set to the circumradius of a pentagonal subunit, which is taken to be 1D u ≡ 5 nm so that the dodecahedron inradius of 1.46D u = 7.3 nm gives an interior volume consistent with that of the smallest T =1 capsids. Assembly simulations were run at least 10 times for each set of parameters, each of which were concluded at either completion, persistent malformation, or 2×10 8 time steps. For all dynamics simulations there were 60 subunits with box size=200×200×200 nm, resulting in a concentration of 12µM.
SAXS profile and hydrodynamic radius estimations. Small Angle X-ray Scattering (SAXS) analysis was performed using CRYSOL [109]. For this analysis the all-atom structure of an SV40 pentameric subunit [105] was aligned into the position of each coarse-grained subunit and the polymer was replaced with a carbon chain. We note that this entails significant simplification: segments of the protein which were not resolved in the crystal structure were not reconstructed and there is no optimization of structure at the all-atom resolution. We believe this approximation is suitable, given that our analysis is limited to the X-ray scattering profile in the small angle regime, which reflects ∼nm scale structural features. Fitting of the scattering profiles was performed using least squares fitting. Hydrodynamic radius analysis was performed using HYDROPRO [85]. This program is capable of calculating the hydrodynamic radius of large protein complexes. In order to perform this analysis, we treated the synthetic polymer as an amino acid chain. Though this is a gross approximation, it has a negligible effect: removing the polymer entirely does not change the R H trend and only alters the magnitude by ∼ 1%. Simulations were visualized using VMD [54].

VI. ACKNOWLEDGMENTS
We gratefully acknowledge Adam Zlotnick for insightful discussion and critical reading of the manuscript. This work was supported by Award Number R01GM108021 from the National Institute Of General Medical Sciences. Computational resources were provided by the NSF through XSEDE computing resources (Longhorn, Maverick, Condor, and Keeneland) and the Brandeis HPCC which is partially supported by NSF-MRSEC-0820492.

A. Model potentials and parameters
The model details are described in Ref. [89]; we summarize them here for the convenience of the reader. In our model, all potentials can be decomposed into pairwise interactions. Potentials involving capsomer subunits further decompose into pairwise interactions between their constituent building blocks -the excluders, attractors, 'Top' and 'Bottom', and ARM pseudoatoms. It is convenient to write the total energy of the system as the sum of 6 terms: a capsomer-capsomer U cc part (which does not include interactions between ARM pseudoatoms), capsomer-polymer U cp , capsomer-ARM U ca , polymer-polymer U pp , polymer-ARM U pa , and ARM-ARM U aa parts. Each is summed over all pairs of the appropriate type: where cap i cap j<i is the sum over all distinct pairs of capsomers in the system, cap i poly j is the sum over all capsomer-polymer pairs, etc. The capsomer-capsomer potential U cc is the sum of the attractive interactions between complementary attractors, and geometry guiding repulsive interactions between 'Top' -'Top' pairs and 'Top' -'Bottom' pairs. There are no interactions between members of the same rigid body, but ARMs are not rigid and thus there are intra-subunit ARM-ARM interactions. Thus, for notational clarity, we index rigid bodies and non-rigid pseudoatoms in Roman, while the pseudoatoms comprising a particular rigid body are indexed in Greek. E.g., for capsomer i we denote its attractor positions as {a iα } with the set comprising all attractors α, its 'Top' positions {t iα }, and its 'Bottom' positions {b iα }. The capsomer-capsomer interaction potential between two capsomers i and j is then defined as: where ε is an adjustable parameter which both sets the strength of the capsomer-capsomer attraction at each attractor site and scales the repulsive interactions which enforce the dodecahedral geometry. N t , N b , and N a are the number of 'Top', 'Bottom', and attractor pseudoatoms respectively in one capsomer, σ t and σ b are the effective diameters of the 'Top' -'Top' interaction and 'Bottom' -'Top' interactions, which are set to 10.5 nm and 9.0 nm respectively throughout this work, r 0 is the minimum energy attractor distance, set to 1 nm, ̺ is a parameter determining the width of the attractive interaction, set to 2.5, and r cut is the cutoff distance for the attractor potential, set to 10.0 nm. The function L is defined as the repulsive component of the Lennard-Jones potential shifted to zero at the interaction diameter: The function M is a Morse potential: The capsomer-polymer interaction is a short-range repulsion that accounts for excluded-volume. For capsomer i with excluder positions {x iα } and polymer subunit j with position R j , the potential is: where N x is the number of excluders on a capsomer and σ xp = 0.5(σ x + σ p ) is the effective diameter of the excluder -polymer repulsion. The diameter of the polymer bead is σ p = 0.5 nm and the diameter for the excluder beads is σ x = 3.0 nm for the T =1 model and σ x = 5.25 nm for the T =3 model. The capsomer-ARM interaction is a short-range repulsion that accounts for excluded-volume. For capsomer i with excluder positions {x iα } and ARM subunit j with position R j , the potential is: with σ xA = 0.5(σ x + σ A ) as the effective diameter of the excluder -ARM repulsion with σ A = 0.5 nm the diameter of an ARM bead. The polymer-polymer non-bonded interaction is composed of electrostatic repulsions and short-ranged excludedvolume interactions. These polymers also contain bonded interactions which are only evaluated for segments occupying adjacent positions along the polymer chain and angular interactions which are only evaluated for three sequential polymer segments. As noted in the main text, electrostatics are represented by Debye Huckel interactions.
where R ij ≡ |R i − R j | is the center-to-center distance between the polymer subunits, q p = −1 is the valence of charge on each polymer segment, and U DH is a Debye-Huckel potential smoothly shifted to zero at the cutoff: : otherwise (9) λ D is the Debye length, l b is the Bjerrum length, and q 1 and q 2 are the valences of the interacting charges. Bonds are represented by a harmonic potential: Angles are also represented by a harmonic potential: where ϑ ijk is the angle formed by the sequential polymer units i, j, k.
The ARM-ARM interaction is similar to the polymer-polymer interaction, consisting of non-bonded interactions composed of electrostatic repulsions and short-ranged excluded-volume interactions. These ARMs also contain bonded interactions which are only evaluated for segments occupying adjacent positions along the polymer chain: : {i, j} nonbonded (12) where R ij ≡ |R i − R j | is the center-to-center distance between the ARM subunits and q i is the valence of charge on ARM segment i. Finally, the ARM-Polymer interaction is the sum of short-ranged excluded-volume interactions and electrostatic interactions: MSM Calculations. The transition probability matrix T(τ ) is calculated by column-normalizing the count matrix C(τ ), in which each element C ji gives the total number of transitions from state i to state j measured at a lag time τ . The count matrix is calculated from the many, relatively short, trajectories run in parallel using the ratcheting procedure described below. After diagonalizing T(τ ), the time-dependent state probabilities can be written as where ω i is the i th eigenvalue of T(τ ) and i| and |i are the corresponding left/right eigenvectors, which are assumed to be normalized. Since T(τ ) is generally not Hermitian the left and right eigenvectors are not equivalent. There is only one unit eigenvalue, whose associated right eigenvector corresponds to the equilibrium distribution, while all other eigenvalues are positive and real [110]. The average value of an order parameter, Q, can then be calculated as a function of time fromQ with Q i as the order parameter value for state i and P i (t; τ ) as the time dependent probability of state i (from Eq. 14). Eq. 15 was used to calculate the mean cluster size n and mean number adsorbed n ad in SI Fig. 4C. In order to calculate the probability of a given cluster size as a function of time in SI Fig. S4, the same procedure was followed except that all states for a given n were lumped together to give a single curve for each n.
Ratcheting Procedure. Because of the Markov property, starting coordinates for simulations can be chosen to efficiently generate good statistics for all of the relevant transition elements. Many simulations are run in parallel for a time t s , which must be longer than the lag time τ , but can be much shorter than the longest relaxation timescale. For the slow nucleation calculation in this paper, we continuously ran 100 simulations in parallel starting from an unassembled bath of subunits. New simulations were seeded from the states with the fewest number of simulation starts, biasing sampling to unexplored regions of state space. Our ratcheting procedure initially grouped states only by the number of subunits adsorbed to polymer n ad . However, after discovering states for each cluster size on the pathway to assembly, we ran a large number of simulations in parallel seeded from a random selection of existing states. This crude, but effective approach could be improved by using adaptive sampling [11], which explicitly considers the error in the MSM when seeding new simulations. A total of 6.3 × 10 9 simulation steps were used to build the MSM.  Figure S1. Linear adsorption density of subunits on the polyelectrolyte (per nm), as a function of subunit concentration for varying salt concentrations, for simulations with no subunit-subunit attraction (εss=0kBT ).    Figure S3. Hydrodynamic radius, RH, as a function of time for linear polyelectrolytes and the model NA described in Ref. [89], which accounts for some effects of base pairing. En masse assembly leads to a larger increase in RH for assembly around both linear (A,B) and base paired (C) polymers, compared to nucleation-and-growth trajectories around linear (D,E) and base-paired (F) polymers.