A theory for viral capsid assembly around electrostatic cores


We develop equilibrium and kinetic theories that describe the assembly of viral capsid proteins on a charged central core, as seen in recent experiments in which brome mosaic virus (BMV) capsids assemble around nanoparticles functionalized with polyelectrolyte. We model interactions between capsid proteins and nanoparticle surfaces as the interaction of polyelectrolyte brushes with opposite charge, using the nonlinear Poisson Boltzmann equation. The models predict that there is a threshold density of functionalized charge, above which capsids efficiently assemble around nanoparticles, and that light scatter intensity increases rapidly at early times, without the lag phase characteristic of empty capsid assembly. These predictions are consistent with, and enable interpretation of, preliminary experimental data. However, the models predict a stronger dependence of nanoparticle incorporation efficiency on functionalized charge density than measured in experiments, and do not completely capture a logarithmic growth phase seen in experimental light scatter. These discrepancies may suggest the presence of metastable disordered states in the experimental system. In addition to discussing future experiments for nanoparticle-capsid systems, we discuss broader implications for understanding assembly around charged cores such as nucleic acids.

I Introduction

The cooperative assembly of heterogenius building blocks into ordered structures is crucial for many biological processes, such as the assembly of protein subunits around a viral nucleic acid to form a capsid. Understanding how nucleic acid-protein interactions drive cooperative assembly could enable antiviral drugs that block this essential step in viral replication. At the same time, engineered structures in which viral capsids assemble around synthetic cargoes show great promise as delivery vehicles for drugs Gupta et al. (2005); Garcea and Gissmann (2004); Dietz and Bahr (2004) or imaging agents Soto et al. (2006); Sapsford et al. (2006); Boldogkoi et al. (2004); Dragnea et al. (2003) and as subunits or templates for the synthesis of advanced multicomponent nanomaterials Chatterji et al. (2005); Falkner et al. (2005); Flynn et al. (2003); Douglas and Young (1998, 2006). Realizing these goals, however, requires understanding at a fundamental level how properties of cargoes and capsid proteins control assembly rates and mechanisms. With that goal in mind, this article presents coarse-grained thermodynamic and kinetic models that describe experiments in which brome mosaic virus (BMV) capsid proteins assemble around spherical charge-functionalized nanoparticle cores.

Identifying mechanisms for simultaneous assembly and cargo-encapsidation. Elegant experiments have studied the assembly of capsids around central cores consisting of nucleic acids (e.g. Johnson et al. (2004a); Fox et al. (1994); Valegard et al. (1997); Johnson et al. (2004b); Tihova et al. (2004); Krol et al. (1999); Stockley et al. (2007); Toropova et al. (2008), inorganic polyelectrolytes Bancroft et al. (1969); Hu et al. (2008); Sikkema et al. (2007), and charge-functionalized cores Sun et al. (2007); Dixit et al. (2006); Chen et al. (2005); Dragnea et al. (2003); Chen et al. (2006); Chang et al. (2008). However, assembly mechanisms of viruses and virus-like particles remain incompletely understood because cargo-protein interactions and the structures of transient assembly intermediates are not accessible to experiments. For the case of empty-capsid assembly, Zlotnick and coworkers achieved great insight about assembly mechanisms with simplified theories (e.g. Zlotnick (2007, 1994); Endres and Zlotnick (2002); Zlotnick et al. (2000); Ceres and Zlotnick (2002); Singh and Zlotnick (2003); Kegel and van der Schoot (2004)), for which a small number of parameters are fit by comparison with experimental assembly data over a range of protein concentrations and pH values. The theories developed in this work enable similar comparisons for the assembly of capsids around electrostatic cores, and could thereby elucidate mechanisms for simultaneous assembly and cargo encapsidation.

Prior theoretical works have led to important insights about assembly around polymeric cores, but were equilibrium studies Belyi and Muthukumar (2006); Angelescu et al. (2006); van der Schoot and Bruinsma (2005); Zhang and Glotzer (2004); ?; Hu et al. (2007) or considered specific RNA-mediated assembly pathways with phenomenological descriptions of protein-nucleic acid interactions Rudnick and Bruinsma (2005); Hu and Shklovskii (2007). Computer simulations of a model in which subunits assemble around rigid spherical cores Elrad and Hagan (2008); Hagan (2008) suggest that a core with a geometry commensurate with the lowest free energy empty-capsid morphology can increase assembly rates and the efficiency of assembly (compared to empty capsid assemblyHagan and Chandler (2006); Rapaport et al. (1999); Rapaport (2004, 2008); Nguyen et al. (2007); Zhang and Schwartz (2006); Sweeney et al. (2008); Schwartz et al. (1998); Wilber et al. (2007); Nguyen and Brooks (2008a)) through the effects of heterogeneous nucleation and templating; however, assembly is frustrated for core-subunit interaction strengths that are too strongElrad and Hagan (2008); Nguyen and Brooks (2008b). Cores that are not size-matched with the lowest free energy capsid morphology can direct assembly into alternative capsid morphologies for optimal core-subunit interaction strengths Elrad and Hagan (2008). Important limitations of the simulations are that extensive comparison with experimental data is computationally intractable, and the simulated models do not explicitly represent electrostatics or the polymeric nature of disordered N-terminal tails in capsid proteins.

The present work aims to overcome these limitations by calculating the electrostatic interaction between charge-functionalized nanoparticles and capsid proteins, and then constructing simplified thermodynamic and kinetic theories that describe the simultaneous assembly and cargo encapsidation process. We explore the effect of cargo-capsid protein interactions on assembly mechanisms by presenting predicted assembly rates and nanoparticle incorporation efficiencies as functions of the functionalized charge density, nanoparticle-protein stoichiometric ratio, and the capsid protein concentration. Because each of these parameters can be experimentally controlled, the theoretical predictions suggest a series of experiments that can elucidate the relationship between cargo-capsid protein interactions and assembly mechanisms.

Although there is currently not sufficient experimental data to estimate values for several unknown parameters (adsorption and assembly rate constants, and subunit-subunit binding energies), we qualitatively compare the theoretical predictions to preliminary experimental time-resolved light scattering data (Fig. 1a) and measurements of the efficiency of nanoparticle incorporation in well-formed capsids (Fig. 1b). With physically reasonable values for the unknown parameters, the predicted light scatter increases rapidly at short times without a lag phase, consistent with the experimental measurements, but the theory does not completely capture the functional form of light scatter at long times. The theory predicts that nanoparticles are incorporated into capsids only above a threshold functionalized charge density, which is similar to the lowest charge density for which significant incorporation is observed in experiments. However, predicted incorporation efficiencies rise to 100% over a narrow range of charge densities, whereas the experimental data shows a gradual increase in the incorporation efficiency with increasing charge density. We discuss experiments and potential improvements to the theory which might elucidate these discrepancies.

This paper is organized as follows. In Section II we calculate the interaction between capsid proteins and functionalized charge on nanoparticle surfaces. In Section III, we present a model for the thermodynamics of capsid assembly in the presence of rigid spherical cores, followed in Section IV by a kinetic theory that describes the simultaneous assembly of capsid subunits on nanoparticle surfaces and in bulk solution. We analyze kinetics and the assembly pathways predicted by the model in section V, and discuss the connection to current and future experiments in Section VI. We also compare theory predictions to simulation results in Appendix A.

Figure 1: (a) Time-resolved light scatter for capsid assembly on nanoparticles with a functionalized surface density of acid groups/nm. (b) Packaging efficiency, or the fraction of nanoparticles incorporated into =3 capsids (measured from TEM microgaphs), varies with . The capsid protein-nanoparticle stoichiometric ratio is , core diameters are 12 nm (which promote =3 capsid formation), the capsid protein concentration is 0.4 mg/ml ( M dimer subunit), and there is 100 mM 1:1 salt and 5 mM 2:1 salt. Packaging efficiencies were measured at seconds. The surface density of acid groups is estimated from the fraction of carboxylated TEG molecules, assuming a total surface density of 3 TEG molecules/nm on the nanoparticle surface. Data thanks to Bogdan Dragnea Dragnea (2008a).

Experimental motivation. This work is motivated by experiments in which BMV capsid proteins assemble around gold nanoparticle cores Chen et al. (2006); Sun et al. (2007) functionalized with thiolalkylated tetraethylene glycol (TEG) molecules, a fraction of which are terminated with carboxylate groups Chen et al. (2006). With 12 nm nanoparticles the TEG-nanoparticle complexes have a diameter of roughly 16 nm, which is nearly commensurate with the interior of a T=3 BMV capsid (17-18 nm Lucas et al. (2002)). Ionized carboxyl groups provide a driving force for capsid proteins to adsorb and assemble on nanoparticle surfaces. Time-resolved light scatter measurements (Fig. 1a), which monitor assembly kinetics, are strikingly different than light scatter of empty capsid assembly at short times (e.g. Chen et al. (2008); Zlotnick et al. (1999); Casini et al. (2004)); instead of a lag phase ( seconds – minutes) there is a rapid ( second) rise in signal. As described in Section V, the theoretical predictions in this work begin to explain this observation. Assembly effectiveness is monitored by using TEM to measure packaging efficiencies, or the fraction of nanoparticles incorporated within well-formed capsids. As shown in Fig. 1b, packaging efficiencies increase from 0% to nearly 85% as the fraction of TEG molecules that are carboxylate-terminated increases from 20% to 100%.

While the light scattering data was obtained from an assembly reaction entirely at = 5, the packaging efficiency data shown in Fig. 1a were obtained with a protocol in which capsid subunits were first dialyzed against a buffer at , followed by buffer at . In this work we consider only , since at disordered protein-nanoparticle complexes form but well-formed capsids are not observed without subsequently lowering . Assembly may not be favorable at because binding interactions between capsid proteins are less favorable than at lower Ceres and Zlotnick (2002); also strong TEG-protein interactions might render disordered states kinetically or thermodynamically favorable.

Figure 2: The model system. Surface functionalization molecules (TEG) are end-grafted to an impenetrable gold sphere. The cationic capsid protein N-terminal arms are modeled as polyelectrolytes grafted to the inner surface of a sphere (the capsid), which is impermeable to TEG and polyelectrolyte but permeable to solvent and ions.
Figure 3: The volume fractions of (a) polyelectrolyte segments and (b) carboxylate groups in the grafted TEG layer are shown as a function of layer () above the nanoparticle surface for several surface acid group densities, with the nanoparticle radius nm and the layer dimension nm.

Ii Subunit adsorption on nanoparticle surfaces

For many capsid proteins, there is a high density of basic residues on flexible N-terminal tails, with numbers of charged amino acids that range from six to 30 Belyi and Muthukumar (2006). These positive charges help drive capsid assembly in vivo by interacting with negative charges on nucleic acids, and the interaction of polymeric tails and nucleic acids is explored in Refs. Belyi and Muthukumar (2006); Hu and Shklovskii (2007); Siber and Podgornik (2007). In this section, we analyze the interactions between capsid protein tails and surface functionalized TEG molecules. For simplicity, we neglect other effects, such as hydrophobic interactions, that might drive subunit adsorption, and we assume that all charges on capsid proteins are located in the polymeric tails. For the remainder of this work, we will reserve the word ‘polyelectrolyte’ to apply to the charged polymeric tails on capsid proteins; we will not use this word to refer to the functionalized TEG molecules.

Model system. We first consider a dilute solution of capsid protein subunits with density , and a single nanoparticle; we consider a finite concentration of cores in Section III. At high surface coverage or in an assembled capsid the folded portion of capsid proteins acts as an impenetrable barrier to the polymeric tails. We therefore represent the nanoparticle-capsid complex as a concave polyelectrolyte brush (the interior surface of the capsid) interacting with a convex brush of the opposite charge (the TEG layer); the model system geometry is illustrated in Fig. 2. We consider surface functionalization molecules grafted to the outer surface of an impenetrable sphere with radius 6 nm. The TEG surface functionalization molecules used in Ref. Sun et al. (2007) consist of a 10 carbon alkanethiol chemically linked to a polyethylene glycol chain of 4 monomer units. The surface density of functionalized charge groups, , is experimentally controlled by terminating a fraction TEG molecules with carboxylates, with the remainder terminated with neutral groups Dragnea (2008a). We represent TEG molecules as freely jointed chains with 9 neutral segments end-grafted to the surface and one charged segment (with variable valence) at the free end.

The N-terminal tails, or ’arms’ of BMV capsid proteins contain nine positive charges in the first 20 amino acids and the first 25-40 residues are disordered in capsid crystal structures Lucas et al. (2002). Since BMV capsids assemble from protein dimer subunits Pfeiffer and Hirth (1974); Cuillel et al. (1983); Berthetcolominas et al. (1987), in our model each subunit has a net charge of . We represent the two N-terminal arms from each dimer subunit as a single polymer of 36 segments with Kuhn length 3.5 , and valence 0.5. These polyelectrolytes are end-grafted to the inner surface of a sphere with radius 9 nm, the inner diameter of a T=3 BMV capsid Lucas et al. (2002). The model capsid is permeable to solvent and ions, but impenetrable to TEG and polyelectrolytes segments. The latter condition is reasonable for high surface coverages of capsid subunits and in the strong coupling limit, in which case the polyelectrolyte configurational entropy is negligible compared to electrostatic effects. For low charge densities, the results do not change qualitatively if the impermeable capsid is eliminated, but the equilibrium concentration of adsorbed polyelectrolytes increases by up to 50%.

Numerical calculations. To our knowledge, interacting polyelectrolyte brushes with this geometry have not been studied, although the interaction of capsid arms with a nucleic acid is considered in Refs. Belyi and Muthukumar (2006); Siber and Podgornik (2008). For many cases we consider, the electrostatic potential in the region of the capsid and nanoparticle surface is large; therefore, we will solve the nonlinear Poisson Boltzmann (PB) equation. To do so, we employ the method of Scheutjens and Fleer (SF) Scheutjens and Fleer (1979); Bohmer et al. (1990, 1991), in which the spatial distributions of polyelectrolyte, TEG, ions, and water, as well as the dissociation equilibrium for carboxylate and water groups, are solved numerically with a self consistent field approximation on a lattice. The calculation accounts for the finite size of all species and the calculated free energies explicitly consider the entropy of ions and solvent molecules. The method is thoroughly described in Refs. Bohmer et al. (1990); Israels et al. (1994); our implementation and additional terms that must be included in the free energy to account for dissociation of weak acid groups are described in Appendix B. We neglect spatial variations of segment densities in the directions lateral to the surface (and thus neglect the effect of ion-ion correlations), and determine the variation of densities in the direction normal to the surface.

All calculations consider the conditions used for the experimental data shown in Fig. 1a, with pH=5, and 100 mM 1:1 salt with an additional 5mM 2:1 salt. The total surface density of functionalized molecules (neutral and charged) is nm at the nanoparticle surfaceDragnea (2008b). Rather than explicitly modeling two species of surface molecules (neutral and acid-terminated), we vary the surface charge density by changing the valency of the ionic end group. To explore a wide range of surface charge we consider end groups with valencies ; the experimental data in Fig. 1 corresponds to the range . For the number of functionalized charge groups corresponds to of the charge on a BMV capsid. The lattice size and statistical segment length for polyelectrolyte and TEG molecules Kawakami et al. (2006); Oesterhelt et al. (1999); Kienberger et al. (2000) are set to nm (which roughly corresponds to the size of a single amino acid and the correlation length of water), and there are 0.5 charges per polyelectrolyte segment, which corresponds to the linear charge density predicted by counterion condensation Manning (1969); Muthukumar (1996); Patra and Yethiraj (1999); Kundagrami and Muthukumar (2008). Except for charge, all species are treated as chemically identical, with the Flory interaction parameter and dielectric permittivity ). Free energies are insensitive to the introduction of unfavorable interactions () between the alkane segments and hydrophilic species or species-dependent dielectric constants. In the strong coupling limit, the equilibrium density of positive charge due to adsorbed subunits is not sensitive to statistical segment length, lattice size, polymer length or charge per segments. The dissociation constant () for carboxylate groups is not known in the vicinity of the surface functionalized gold particle; we use , so that approximately 25% of the acid groups are dissociated for at full coverage (see below), consistent with Ref. Gershevitz and Sukenik (2004). We do not explicity consider complexation of acid groups, a possibility suggested in that reference. We will also consider experiments in which the terminal group is a strong acid, with .

Figure 4: The free energy per area at the capsid surface for the adsorption of a concave polymer brush onto a spherical brush with opposite charge. For varying surface densities of polyelectrolyte (normalized by the density in a complete capsid), , the free energy relative to is shown for nanoparticle functionalization with weak acid groups () at acid groups/nm (), () and strong acid groups at (+), (o). There is 100 mM 1:1 and 5 mM 2:1 salt and .
Figure 5: The equilibrium surface density of adsorbed polyelectrolytes, , normalized by the surface density of polyelectrolyte on a core with a complete capsid. Predictions are shown for surface functionalization with weak acid () and strong acid (+) groups. The bulk subunit concentration is 10 M, with other parameters as in Fig. 4.

Free energies and equilibrium concentrations of adsorbed polyelectrolytes. The equilibrium surface density of adsorbed polyelectrolyte can be calculated by considering a system for which there is free exchange between surface adsorbed polyelectrolytes and a bath containing polyelectrolyte, salt, and solvent. The kinetic theory developed in Section IV, however, requires the free energy for non-equilibrium adsorption densities. We therefore consider a restricted number of grafted molecules (which cannot exchange with the bath), as a function of surface charge and density of grafted molecules. For each grafting density , the spatial distributions of polyelectrolyte, TEG, ions, and solvent are determined by minimizing the total free energy. The free energy, , is calculated from Eqs. 31 - 37. Note that the density of polyelectrolytes (Fig. 3) segments is nonmonotonic with distance from the capsid, consistent with studies of capsid assembly around nucleic acids Belyi and Muthukumar (2006); Siber and Podgornik (2008).

Free energies for two functionalized charge densities with weak and strong acid groups are shown in Fig. 4 and the equilibrium adsorption densities, , which correspond to the minimum excess free energy, are shown in Fig. 5. The adsorption densities (, ) are normalized with respect to the density of dimer subunits in a complete capsid (0.09 subunits/nm on the interior surface of the capsid with inner radius 8.9 nm). The predicted adsorption densities for weak acids () are significantly lower than for the case of strong acids because the high local density of negative charges shifts the carboxylate dissociation equilibrium; for a carboxylate surface density the average dissociated fraction shifts from 75% for an isolated molecule to 35%. However, dissociation is enhanced as polyelectrolytes with positive charges adsorb and the average dissociation fraction on a core with a complete capsid is 78%. As discussed in Section V, this effect results in significantly different dependencies of assembly kinetics on for weak and strong acids.

There are several reasons why our calculations may underestimate at low and for weak acid groups. Nanoparticle-protein interactions in addition to those involving the N-terminal arms could contribute to subunit adsorption. The calculated interactions involving N-terminal arms assume uniform charge densities in directions parallel to the surface, which will overestimate polyelectrolyte interactions at low surface densities, when the average distance between adsorbed subunits is greater than the polyelectrolyte radius of gyration; the approach described for neutral subunits in Appendix A might be appropriate in this regime. As discussed above, reducing restrictions on polyelectrolyte conformations due to the folded portion of capsid proteins results in somewhat higher adsorption densities for low surface charges. In addition, the of carboxylate groups and the effect of interactions between the charged species and the gold surface are poorly understood. For this reason, we examine assembly over a wide range of surface charges for both and to understand the interplay between surface charge, dissociation, and assembly. It would be desirable to experimentally measure the surface charge as a function of the fraction of charged TEG molecules on the surface; while it will be difficult to quantitatively verify the surface charge or corresponding magnitude of the electrostatic potential, the variation of the potential at the core surface with (see Fig. 6) is roughly consistent with preliminary experimental measurements of zeta potentials for functionalized nanoparticles Dragnea (2008b). While our calculations may underestimate , these approximations are less important at high coverages and hence should not significantly affect equilibrium packaging efficiency calculations.

Figure 6: The electrostatic potential at the nanoparticle surface for varying surface functionalization () with strong (+) and weak () acid groups. There is 100 mM 1:1 and 5 mM 2:1 salt.
Figure 7: The free energy per area at the capsid surface of an empty capsid due to electrostatic repulsions and ion and solvent entropy is shown as a function of the charge density on capsid subunits, as predicted by SF calculations (solid lines with symbols) and the linearized PB equation (dashed line). The calculations use 100 mM 1:1 salt and 5 mM 2:1 salt.

ii.1 Empty-capsid free energies.

We use the same approach, but without the core and TEG molecules, to estimate the electrostatic contribution to the free energy of forming an empty capsid. The free energy, , as a function of polyelectrolyte density () is shown in Fig. 7, where it is compared to the free energy calculated with the linearized Poisson Boltzmann equation for charge smeared on a spherical surface Kegel and van der Schoot (2004); Siber and Podgornik (2007) with radius nm, , with subunits, =18 charges/subunit, nm the Bjerrum length of water, and the Debye length. Although the linearized Poisson Boltzmann equation was shown to agree closely with the full nonlinear calculation in Ref. Siber and Podgornik (2007) for this ionic strength, there is a large discrepancy between the PB and SF calculations; the SF calculations predict much lower free energies because the flexible tails adopt stretched conformations which reduce electrostatic repulsions (see also Ref. Siber and Podgornik (2008)). We note, though, that we have not considered interactions between charges located in the structured region of capsid proteins.

Zlotnick and coworkers Ceres and Zlotnick (2002); Ceres et al. (2004); Zlotnick (2007) have shown that capsid formation free energies can be fit using the law of mass action to experimental data for capsid and free subunit concentrations as the total capsid protein concentration is varied. The resulting free energies depend on H as well as salt concentration Ceres and Zlotnick (2002) and likely include effects arising from hydrophobic as well as electrostatic interactions  Kegel and van der Schoot (2004). We can subtract the polyelectrolyte free energy from the total capsid formation free energy derived from experimental data to estimate the free energy due to attractive subunit-subunit interactions, , which may arise from hydrophobic contacts. Following Ceres and coworkers Ceres and Zlotnick (2002), we assume that each subunit in a capsid makes favorable contacts with neighboring subunits, so the attractive energy per subunit is


For protein dimer subunits in a =3 BMV capsid with an inner capsid radius of nm Lucas et al. (2002), we obtain . With a typical subunit-subunit contact free energy of (-3 kcal/mol) Ceres and Zlotnick (2002); Ceres et al. (2004); Zlotnick (2007), this gives (compare to the estimate from experimental data for hepatitis B virus assembly by Kegel and van der Schoot Kegel and van der Schoot (2004) of an attractive interaction strength of per subunit).

Iii Equilibrium theory for core-controlled assembly

In this section we consider the thermodynamics for a system of subunits that can assemble into capsids around cores and assemble in bulk solution to form empty capsids. The free energies due to electrostatic interactions for subunits adsorbed on surfaces or assembled into empty capsids developed in section II are used to determine the relative free energies of core-associated capsid intermediates.

iii.1 The thermodynamics of core-controlled assembly

We consider a dilute solution of capsid protein subunits with density , and solid cores with density . Subunits can associate to form capsid intermediates in bulk solution or on core surfaces. A complete capsid is comprised of subunits. Zandi and van der Schoot ? develop the free energy for core-controlled assembly in the context of capsids assembling around polymers, but neglect the contribution of partial capsid intermediates, which will be important for describing the kinetics of core-controlled assembly in the next section. Here, we consider the free energy for a system of subunits, capsids, and partial capsid intermediates. To simplify the calculation, we neglect the possibility of malformed capsids, but note that disordered configurations could be kinetically arrested or even the favored equilibrium state in the case of core-subunit interactions that are much stronger than the thermal energy Elrad and Hagan (2008); Nguyen and Brooks (2008b).

The total free energy density is given by


with the free energy density of empty capsid intermediates


where is the density of intermediates with subunits and is a standard state volume. For simplicity we follow Zlotnick and coworkers Zlotnick (1994); Endres and Zlotnick (2002) and assume that there is only one intermediate for each size , with a free energy defined below in Eq. 9.

The free energy density of core-associated intermediates is given by


where is the fraction of cores with an adsorbed intermediate of size and with adsorbed but unassembled subunits. The free energy for such a core is given by and is defined below. For simplicity, we assume that a core can have at most one assembled intermediate and a total of subunits can fit on a single core surface, where is the size of a complete capsid, so states with subunits have zero probability. The prime on the second sum indicates that because, by definition, a state with has no assembled intermediates.

To obtain the equilibrium concentration of intermediates we minimize the total free energy subject to the constraints




where Eq. 5 enforces that the total density of subunits is given by .

Minimization yields




where is the inverse of the thermal energy, and are chemical potentials for subunits and cores, respectively.

iii.2 Free energies for empty-capsid assembly intermediates

The free energy of an empty capsid intermediate with subunits is written as


where is the number of new subunit-subunit contacts formed by the binding of subunit to the intermediate, is the contact free energy, accounts for degeneracy in the number of ways subunits can bind to or unbind from an intermediate (see the factors in Refs. Zlotnick (1994); Endres and Zlotnick (2002)), and is a rotational binding entropy penalty. For simplicity and to ease comparison with earlier works, we take except when comparing theoretical predictions to simulation results (see Appendix A). As discussed at the end of Section II, the subunit-subunit contact free energy has been estimated from experiments for several viruses Ceres and Zlotnick (2002); Ceres et al. (2004); Zlotnick et al. (1999); Zlotnick (2007) and depends on salt concentration and H Ceres and Zlotnick (2002); Kegel and van der Schoot (2004).

iii.3 Free energies for assembly intermediates on core surfaces

The free energy for a core without a partial capsid is given by with the core area, the functionalized acid group density, the normalized density of adsorbed subunits, and calculated in Section II. Likewise, the free energy for a complete capsid is , where the polyelectrolyte contribution to the capsid free energy () is subtracted because it is included in .

In the case of a core with a partially assembled capsid and adsorbed but unassembled subunits, there is an inhomogeneous distribution of charge on the surface, with higher charge density in the region of the assembled cluster. We determine the free energy for such a core, with adsorbed subunits and an intermediate of size , in the following schematic way. We assume the unassembled subunits spread on the surface of the core not occupied by the intermediate. There are thus two ‘pools’ of charge on the surface – the partial capsid with a high subunit density and the remaining area of the core surface with a charge density spread uniformly over the area not covered by the intermediate , . The total electrostatic free energy is then given by


In general there should be a term in to account for strain energy if the core curvature is not commensurate with that of the intermediate; we assume that capsid and core curvatures match.

Core encapsidation. Above a threshold or critical subunit concentration (CSC), the majority of subunits will be in capsidsZlotnick (1994); Bruinsma et al. (2003); Hagan and Chandler (2006). Because core-subunit interactions can lead to high localized surface charge densities, the threshold concentration for assembly in the nanoparticle system can be below the CSC for spontaneous assembly. Above the CSC, both empty and core-filled capsids can assemble. Due to the entropy cost of assembling a capsid on a core, Eqs. 8 and 7 show that there is a threshold surface free energy ( for a stoichiometric amount of capsid protein and nanoparticles), below which empty capsids and free cores are favored over full capsids. Above the threshold surface free energy and above the CSC, nearly all cores will be encapsidated () at equilibrium. However, we will see that there can be a higher threshold surface free energy for efficient encapsidation kinetics.

Iv Kinetic theory for core-controlled assembly

In this section we use the free energies developed in section III for capsids and capsid intermediates as the basis for a kinetic theory of core-controlled assembly.

Zlotnick and coworkers Zlotnick (1994); Endres and Zlotnick (2002) have developed a system of rate equations that describe the time evolution of concentrations of empty capsid intermediates


where and are respectively the forward and reverse binding rates and is a statistical factor that describes the number of different ways to transition from intermediate to Zlotnick (1994); Endres and Zlotnick (2002); corrects for double counting the number of pairwise combinations of free subunits. Following Ref. Endres and Zlotnick (2002), we simplify the calculation by requiring that transitions between intermediates are only allowed through binding or unbinding of a single subunit and there is only one intermediate for each size . We assume that the assembly rate constant is independent of intermediate size , but the number of contacts per subunit is not (see Table 1).

We extend this approach to describe simultaneous capsid assembly and adsorption to spherical cores, giving the following expression for the time evolution of core states


The first two lines in each of Eqs. 12 describe adsorption and desorption from the core surface with respective rate constants and , while the following lines describes binding and unbinding to a capsid intermediate. Rate constants for binding and unbinding of adsorbed subunits to an adsorbed intermediate are denoted by and , respectively. Finally, and are respectively the rate constants for the association (dissociation) of non-adsorbed subunits to (from) an adsorbed intermediate. The latter process becomes important when an adsorbed capsid nears completion, and electrostatic repulsions render non-specific adsorption of subunits unfavorable, causing the surface capture and diffusion mechanism to be slow in comparison to the empty-capsid assembly mechanism.

Adsorption is calculated from , with the adsorption rate constant and the probability that an adsorbing subunit is not blocked by previously adsorbed subunits. For simplicity, we assume Langmuir kinetics, , but note that this is a poor approximation at high surface coverages (see Appendix A).

The desorbtion rate is related to the adsorption rate by detailed balance


Assembly and disassembly of adsorbed subunits is described in a manner analogous to that used for empty capsids, with assembly rates given by


where the effective surface density is .

Eq. IV states that subunit-subunit binding rates are proportional to the frequency of subunit-subunit collisions, which can be calculated from the Smoluchowski equation for the diffusion limited rate in three dimensions, with the density of subunits () within a layer above the surface with thickness of the subunit size, , or from the diffusion limited rate in two dimensions Berg and Purcell (1977) with a surface density given by . The result of the two calculations differs only by a logarithmic factor that is of order 1 in this case (see appendix B of Ref. Berg and Purcell (1977)). Surface assembly rate constants () can be smaller than those for empty capsid assembly (), if interactions between subunits and the TEG layer impede diffusion or reorientation of adsorbed subunits. Association rates could also change in comparison to bulk solution if subunits take different conformations on the core surface; however, imaging experiments show that protein structures in nanoparticle-filled capsids are nearly identical to those in empty capsids Sun et al. (2007), which suggests that subunits do not denature on the surface.

The disassembly rate is given by detailed balance


Similarly the association and dissociation rates for solubilized subunits are related by detailed balance:


The expression for the time evolution of empty capsid intermediates is now


We note that although Eqs. 12 and 17 have the form of a Master equation, the finite concentration of subunits introduces a nonlinear dependence of the time evolution on state probabilities.

V Results

Figure 8: Packaging efficiencies, or the fraction of nanoparticles incorporated in well-formed capsids, at varying surface charge densities () of weak acid groups predicted by (a) the equilibrium theory and (b) the kinetic theory at seconds. Parameters are as listed in Fig. 1 with , , and subunit binding free energies are indicated on the plots.

In this section we report the predictions of the equilibrium and kinetic theories for two experimentally observable quantities: packaging efficiencies, or the fraction of nanoparticles encapsulated by well-formed capsids (estimated from TEM experiments), and time-resolved light scattering intensity. Preliminary experimental data for these quantities is shown in Fig. 1. Since there is not sufficient data to estimate the values of all unknown parameters, we make qualitative predictions about the effect of changing the experimental control parameters: the functionalized surface charge density , the stoichiometric ratio of capsid protein to nanoparticles , and the capsid protein concentration . The important unknown parameters are the subunit-subunit binding energy, , the dissociation constant of the functionalized charge groups, the rate constant for adsorption of subunits onto the nanoparticle surface, , and the rate constant for the assembly of surface-adsorbed subunits, . The definitions and values of model parameters are given in Table 1.

v.1 Packaging efficiencies

We begin by discussing the relationship between the density of functionalized charge groups on nanoparticle surfaces, , and packaging efficiencies. As discussed in Section I, the experimental measurements of packaging efficiencies were obtained from an experimental protocol in which the solution varied over time; because we do not know the dependence of the subunit-subunit binding energy on , we only present predictions of the equilibrium and kinetic theories for packaging efficiencies at for several binding energies. We will compare these predictions to the currently existing data, but note that they should be compared to packaging efficiencies from experiments carried out entirely at , as was the case for the light scattering data shown in Fig. 1b.

As shown in Fig. 8a, the equilibrium theory predicts effective packaging only above a threshold functionalized charge density, which depends on the binding free energy. Effective packaging even occurs at a binding free energy , showing that assembly on cores can occur under conditions for which spontaneous empty-capsid assembly is not favorable, as seen in experiments Sun et al. (2007) and simulations Hagan (2008); Elrad and Hagan (2008). On core surfaces nanoparticle-subunit electrostatic interactions and attractive subunit-subunit interactions overcome subunit-subunit electrostatic repulsions to drive assembly ( kcal/mol per subunit for , Eq. 1).

As shown in Fig. 8b, for weak subunit-subunit interactions ( kcal/mol) the kinetic theory predictions at seconds are essentially identical to the equilibrium results. As subunit-subunit interactions increase in magnitude, however, the predicted threshold charge density becomes larger than the equilibrium value, and even begins to increase for . With stronger subunit-subunit binding interactions, the assembly of empty capsids competes with assembly on core surfaces by depleting the concentration of free subunits. Assembly on cores remains favorable at high surface charge densities because initial subunit adsorption leads to rapid nucleation, and subsequent nonspecific adsorption enhances capsid growth rates (see below). The predicted packaging efficiencies at long times are relatively insensitive to variations of the kinetic parameters and .

Interestingly, the threshold surface charge density predicted by the kinetic theory, , is nearly independent of the subunit binding energy for kcal/mol, and comparable to the lowest experimental charge densities with measurable packaging efficiencies (Fig. 1). However, the experimental results show a gradual increase in packaging efficiency as the functionalized charge density increases. We do not observe such a gradual increase in packaging efficiencies except for shorter observation times or an extremely low surface assembly rate constant. This discrepancy could arise because we have assumed monodisperse cores that are perfectly matched to the capsid size, or it may indicate the presence of kinetically frustrated states (see below).

Figure 9: The predicted time dependence of light scattering for functionalization with strong acid groups. (a) and (b) Indicated functionalization densities with rate constants (a) and (b) . (c) Indicated surface assembly rate constants, with nm. Other parameters for (a)-(c) are ,, M, and kcal/mol.
Figure 10: The predicted time dependence of light scattering for surface charge groups with and , for (a) indicated functionalization densities , (b) indicated functionalization densities with , and (c) indicated surface assembly rate constants with nm. Other parameters are as in Fig. 9.

v.2 Assembly kinetics on core surfaces

In this section we explore the effect cargo-capsid protein interaction strength on assembly kinetics. As discussed in Section I, assembly kinetics can be monitored in experiments with time-resolved dynamic light scattering; preliminary experimental data for the highest functionalized charge density is shown in Fig. 1. Because it is unlikely that the light scattering signal can distinguish between disordered and assembled subunits on core surfaces, we estimate light scattering from the kinetic theory by calculating the mass-averaged amount of adsorbed and assembled () subunits on core surfaces. We simplify the presentation by varying only the charge density () and the surface assembly rate constant (), with physically reasonable values for the subunit binding free energy, kcal/mol, and the the subunit adsorption rate constant, kad (). As discussed in section Section IV, the effect of reduction in dimensionality due to adsorption onto a surface is accounted for in the kinetic theory, but assembly rate constants could still be dramatically different from those in bulk because of impeded diffusion of adsorbed subunits. We therefore explore a wide range of surface assembly rate constants .

Predicted light scatter as a function of time is shown for various functionalized charge densities and surface assembly rate constants in Figs. 9 and  10 for functionalization with strong and weak acidsstr (), respectively. The analysis is easiest for the case presented in Fig. 9a (), for which assembly on the core surface is slow compared to subunit adsorption. For each value of , there is a rapid initial rise in light scatter, followed by a charge-dependent plateau, and finally a slow increase to saturation. The initial fast kinetics correspond to nonspecific subunit adsorption (i.e. with little or no assembly) and the plateau occurs when adsorption saturates at what would be the equilibrium surface density if there were no assembly (, see Section II). If adsorbed subunits assemble into a relatively stable intermediate, the chemical potential of remaining adsorbed subunits decreases, resulting in further adsorption and increased light scatter. This process continues as additional subunits bind to the growing intermediate, until assembly stops.

At high functionalization densities, nonspecific subunit adsorption approaches the density of a complete capsid; the high local or surface concentration of subunits enables rapid nucleation and assembly and thus the light scatter reaches its final saturation value at relatively short times. At the lower functionalization densities, nonspecific adsorption saturates at smaller surface concentrations, with corresponding smaller initial plateau heights for the light scatter. Nucleation times increase dramatically as the surface concentration is reduced (see the discussion in Ref. Hagan (2008)), resulting in the long plateau before assembly completes.

For a larger assembly rate constant , the nonspecific adsorption and assembly phases can only be clearly distinguished at low functionalized charge densities, as shown in Fig. 9b; at larger adsorption and assembly are concomitant. For weak acid groups (Fig. 10), equilibrium surface densities increase only slowly with (, see Fig. 5) because the local concentration of negative charge disfavors carboxylate dissociation. There is still a strong dependence of overall assembly rates on , however, because the positive charge from capsid proteins drives additional carboxylate dissociation (see Section II). Thus, assembly effectiveness cannot be predicted from alone, as was found for simulations of neutral subunits in Ref. Elrad and Hagan (2008).

Cores can enhance assembly rates. From the preceding discussion, it is clear that functionalized charge leads to an increased concentration of capsid proteins nanoparticle surfaces relative to that in bulk solution, which can increase capsid nucleation rates and even induce assembly below the threshold protein concentration for formation of empty capsids. Because nonspecific adsorption onto core surfaces can be much faster than subunit-subunit binding rates, the core surface can also enhance assembly rates during the growth phase by increasing the effective subunit-subunit interaction cross-section. This effect was hypothesized for RNA by Hu et al. Hu and Shklovskii (2007) and is reported for simulation results in Ref. Hagan (2008).

Assembly rates decrease as capsids near completion. Net assembly rates decrease as capsids near completion because adsorption of subunits onto cores with large partial capsids is impeded by two effects. First, the excluded volume of the partial capsid decreases adsorption rates; a similar effect has been observed in simulations representing empty capsid assembly, but the magnitude of rate reduction seems to be model dependent (see Refs. Nguyen et al. (2007); Hagan (2008)). Electrostatic repulsions with adsorbed and/or assembled subunits on the core surface can have a larger effect on subunit adsorption. When the subunit surface density is large compared to , the surface capture and diffusion assembly mechanism that enables rapid core-controlled assembly becomes inefficient. Even though direct association of free subunits onto core-bound intermediates is included in the kinetic theory, there is a reduction in assembly rates as capsids near completion; this effect is most noticeable for surface functionalization with weak acid groups (Fig. 10a), or low fractions of strong acid groups ( nm). The decreased assembly rates are reflected in slow growth of light scatter as assembly nears completion.

Figure 11: Estimated light scatter from simulations in which subunits assemble around model nanoparticle with a size mismatch. The preferred empty-capsid morphology is a T=3 capsid, but the lowest free energy configuration is a T=1 capsid assembled around the T=1-size nanoparticle. Estimated light scatter is shown as a function of time for different values of the parameter (in units of ), which gives the free energy cost for adopting a subunit conformation consistent with a T=1 morphology. As described in Ref. Elrad and Hagan (2008), the larger values of lead to the formation of frustrated nanoparticle-associated partial capsids with a T=3 morphology, which cannot close around the nanoparticle and block assembly into the lowest free energy configuration. The data is replotted from simulations reported in Fig. 6 of Ref. Elrad and Hagan (2008), with the nanoparticle-subunit interaction energy .

Comparison with experimental data. As discussed above, there is not sufficient experimental data to estimate values for the unknown parameters; however, the predicted light scatter plots share some features with the experimental light scatter data shown in Fig. 1b. In particular, for some parameter sets light scatter increases rapidly at short times, without the lag time seen for empty capsid assembly Zlotnick et al. (1999); Casini et al. (2004); Chen et al. (2008). It is tempting to attribute the initial rise to nonspecific adsorption that does not saturate until nearly complete coverage, as seen for strong acid surface groups (Fig. 9a), but the calculation with weak acid groups predicts values of that are well below complete coverage even for nm. Interestingly, the predicted light scatter curves in this case with fast assembly kinetics (, Fig. 9a) still demonstrate a rapid initial rise followed by slow saturation; the initial rise corresponds to rapid nucleation and assembly concomitant with adsorption. Thus, as discussed in the next section, comparison of theoretical predictions and experimental light scatter data at varying functionalized charge density will be necessary to unequivocally determine the origin of the initial light scatter phase.

The experimental light scattering data (Fig. 1b) appears to demonstrate a logarithmic increase in light scattering between 20 and 200 seconds. While the theoretical light scatter curves also demonstrate slow growth near completion (see above) due to subunit-subunit electrostatic repulsions, this effect is limited to the last five or 10 subunits. We did not find any parameter sets for which predicted light scatter grew logarithmically over several decades in time, although our parameter search was not exhaustive. Intriguingly, the predicted light scatter from simulations in which frustrated off-pathway intermediates block assembly Elrad and Hagan (2008) do show logarithmic growth over several decades (see Fig. 11).

v.3 Concentration effects

Varying nanoparticle concentrations. We first examine assembly effectiveness as the nanoparticle concentration is varied at fixed subunit concentrations. We define the capsid protein-nanoparticle stoichiometric ratio as so that when there is exactly enough capsid protein to assemble a complete capsid on every core. The variation of packaging efficiencies with the stoichiometric ratio is shown in Fig. 12; the predicted curves are sigmoidal in shape, much like the experimentally observed packaging efficiencies shown in Fig. 1A of Ref. Sun et al. (2007). Packaging efficiencies for are typically lower than the equilibrium values obtained by solving Eqs. 7 and 8, because initially capsid proteins undergo nonspecific adsorption onto every core, which depletes the concentration of free subunits. The formation of complete capsids then requires desorption and subsequent re-adsorption onto cores with large intermediates. The timescale for subunit desorption increases with functionalized surface charge density and therefore the kinetic packaging efficiency is nonmonotonic with respect to for . This scenario is consistent with the slow assembly kinetics experimentally observed for capsid assembly in the presence of RNA Johnson et al. (2004a).

Varying subunit concentrations. As shown in Fig. 13, at acid groups/nm the kinetic theory predicts a much weaker dependence of assembly effectiveness on subunit concentration than has been seen for empty-capsid assembly. This trend arises because the chemical potential of adsorbed subunits is dominated by electrostatics rather than subunit translational entropy.

Figure 12: Packaging efficiencies depend on the capsid protein/core stoichiometric ratio. The predictions of the kinetic theory for the variation of packaging efficiencies with the capsid protein/core stoichiometric ratio, at a fixed subunit concentration of M, are shown for functionalization with weak acid groups, with other parameters as in Fig. 8.
Figure 13: Packaging efficiencies depend only weakly on the initial subunit concentration, , for fixed . The predictions of the kinetic theory for the variation of packaging efficiencies with subunit concentration are shown for functionalization with weak acid groups, with other parameters as in Fig. 12.

Vi Discussion

The feasibility of estimating the unknown parameters. The results in the last section are consistent with some features of the preliminary experimental data, and the predicted dependencies of packaging efficiency and light scattering on the functionalized charge density and nanoparticle-protein stoichiometric ratio can be qualitatively compared with additional experimental data. Quantitative comparisons of experimental data, however, will require estimates for at least three currently unknown parameters, the subunit-subunit binding energy , the assembly rate constant on core surfaces , the subunit adsorption rate . In Appendix A, we show that, with good estimates for the unknown parameters, the theoretical predictions agree reasonably with simulation results for a model of core-controlled assembly. While it is not certain that a simplified theory with a few parameters will exactly match experimental data from a system with thousands of distinct intermediates, each of which could have different rate constants, rate equations with a similar level of coarse graining have successfully explained many features of empty-capsid assembly (e.g. Zlotnick (2007, 1994); Endres and Zlotnick (2002); Zlotnick et al. (2000); Ceres and Zlotnick (2002); Singh and Zlotnick (2003)). Given the number of experimental control parameters in the nanoparticle-capsid assembly system, it should be possible to generate sufficient experimental data to further test the qualitative trends predicted here, and to make quantitative estimates for the unknown parameters through data fitting. With estimated values for and , it would be possible to predict the assembly state of subunits on core surfaces during the initial rise in light scatter, which would be challenging to determine with experiments alone.

In addition to measuring packaging efficiencies and light scattering data over a range of functionalized charge densities (with both weak and strong acid groups) and capsid protein-nanoparticle stoichiometric ratios, the theoretical predictions in the last section suggest additional experiments that could be useful to provide independent estimates of some unknown parameters. In particular, the predicted light scatter curves in Figs. 9 and 10 show that adsorption and assembly happen simultaneously for many parameter sets. Measuring light scatter of assembly incompetent proteins in the presence of functionalized nanoparticles would enable independent estimation of the subunit adsorption rate and the equilibrium subunit surface charge density ( Fig. 5). In addition to varying the functionalized charge density, the subunit-nanoparticle interaction could modified by mutations to the N-terminal arm on capsid proteinsAniagyei et al. (2009); Tang et al. Tang et al. (2006) have shown that cowpea chlorotic mottle virus capsid proteins with 34 residues deleted from the N-terminal arms are assembly competent (although they lose selectivity for the native capsid structure).

Identifying metastable disordered states. We make several important approximations in this work. Most significantly, the theory considers only only intermediates consistent with the native capsid geometry, while experiments Dragnea (2008b) and simulations Elrad and Hagan (2008); Nguyen and Brooks (2008b) indicate that other capsid morphologies and/or asymmetric, malformed structures can be kinetically or even thermodynamically favored when core-subunit interaction strengths are large compared to the thermal energy . These effects could be accounted for in the theory by extending the state space to include structures other than complete and partial well-formed capsids, and by introducing subunit diffusion constants that depend on the strength of core-subunit interactions. If the present theory proves capable of describing experiments with parameters for which predominately well-formed capsids assemble, however, inconsistencies between theoretical predictions and experimental measurements for other system parameter values could identify kinetic traps. For example, the calculated packaging efficiencies show a less gradual increase with functionalized charge density than measured in experiments and predicted light scatter does not fully capture the a logarithmic increase in the experimental light scatter curve. Both of these discrepancies could be explained by frustrated off-pathway of assembly intermediates or slow subunit diffusion rates, as evidenced by predicted light scatter from simulations which do demonstrate frustrated assembly (Fig. 11). Confirming this possibility will require additional experimental data.

Parameter Value Definition
90 Number of subunits in a complete =3 capsid
995.4 nm Inner surface area of a complete capsid (Section III.3
4.2 nm Subunit diameter
kcal/mol Binding free energy per subunit-subunit contact (Eq. 9)
kcal/mol Free energy per subunit-subunit contact due to attractive interactions (Eq. 1)
Subunit adsorption rate constant
Empty capsid assembly rate constant Johnson et al. (2005) (Eq. 11)
Surface assembly rate constant (Eq. IV)
charge/(nm gold surface) Surface-density of functionalized charge
18 Number of positive charges per protein-dimer subunit
10 M Total subunit concentration
Capsid protein–nanoparticle stoichiometric ratio,
M Nanoparticle concentration
Table 1: Parameter values used for calculations in this work. The parameters used for capsid assembly are as follows. The number of subunit-subunit contacts are for , for , for , and . This choice obviates the need for different nucleation and elongation rate constants Endres and Zlotnick (2002). The nucleus size is . The forward and reverse reaction degeneracies are , for (the average value calculated from simulations in Ref. Hagan (2008)) and , . The binding degeneracy entropy (Eq. 9) is .

Vii Conclusions

In summary, we develop simplified thermodynamic and kinetic theories that describe the simultaneous assembly of viral proteins and encapsidation of cargo. When applied to the encapsidation of charge-functionalized nanoparticles, the theories predicts a transition from inefficient core encapsidation to nearly 100% incorporation efficiency as the functionalized charge density is increased beyond a threshold value, and the estimated light scatter signal shows an initial rapid increase followed by slow rise to saturation, as seen in experiments. The predicted increase in incorporation efficiency with increasing charge density is sharper than seen in experiments, which could indicate the presence of kinetic traps that are not accounted for in the present theory; comparison of theory predictions with experimental data collected with a wide range of control parameters will be important to assess this possibility. The trends predicted here for varying surface charge, subunit concentration, and capsid-nanoparticle stoichiometric ratio could guide the design of experiments that identify fundamental principles and/or additional complexities for simultaneous assembly and cargo encapsulation.

Acknowledgments I am indebted to Bogdan Dragnea and the Dragnea group for insightful conversations and for providing me with experimental data, and to Oren Elrad for Fig 11. Funding was provided by Award Number R01AI080791 from the National Institute Of Allergy And Infectious Diseases and by the National Science Foundation through the Brandeis Materials Research Science and Engineering Center (MRSEC).

Appendix A

Figure 14: The kinetic theory predictions (smooth lines) and simulation results (noisy lines) for the time dependence of adsorbed and assembled subunits, , are shown for the neutral subunit model (Eq. 18). The simulations consider T=1 capsids and commensurate cores, so a complete capsid is comprised of dimer subunits. Curves at increasing height correspond to reduced subunit densities of 2.04, 4.07, 8.14, 20.4, 40.7, with a surface free energy of and subunit-subunit binding energies of (a) and (b)

In this appendix we compare predictions of the kinetic theory for the time dependence of the number of adsorbed particles, , to results presented in Ref. Hagan (2008) from a computational model that represents the assembly of T=1 capsid-like objects around a commensurate nanoparticle. This comparison elucidates the effect of approximations used in the kinetic theory. In particular, in rate equation approaches that assume a single reaction pathway as described above, or multiple reaction pathways Endres et al. (2005); Sweeney et al. (2008); Zhang and Schwartz (2006), and even binding between large intermediates Sweeney et al. (2008); Zhang and Schwartz (2006), only intermediates that are consistent with partially assembled capsids are considered. The simulations, on the other hand, explicitly track the spatial coordinates of subunits and thereby make no a priori assumptions about assembly pathways. Because we can either calculate parameters or estimate them from independent simulations, this comparison is a stringent test of the validity of these approximations. It also tests the utility of using the kinetic theory to describe experimental adsorption results. The procedure by which we determine parameters for the kinetic theory that correspond to particular sets of simulation parameters is described below, as well as modifications to the theory in order to represent neutral subunits (or high salt concentrations). We do not consider here the simulations reported in Fig. 11 and Ref. Elrad and Hagan (2008), for which the nanoparticle is not commensurate with the preferred capsid morphology.

Adsorption kinetics predicted by the kinetic theory and observed in simulations are shown in Fig. 14 for several subunit concentrations and subunit-subunit binding energies for an adsorption energy of . The agreement between the kinetic theory and simulation results is surprisingly good, considering that no parameters were adjusted to fit this data and, as discussed below, relatively crude estimates are used for the subunit-subunit binding rate constants and binding entropy .

The computational model considered in Ref. Hagan (2008) consists of rigid subunits, for which excluded volume interactions are modeled by spherically symmetric repulsive forces, and complementary subunit-subunit interactions that drive assembly are modeled by directional attractions. The lowest energy states in the model correspond to ”capsids”, which consist of multiples of 60 monomers in a shell with icosahedral symmetry. The parameters of the model are the energy associated with the attractive potential, , which is measured in units of the thermal energy , and the specificity of the directional attractions, which is controlled by the angular parameters and . Capsid subunits also experience short ranged isotropic interactions with a rigid sphere placed at the center of the simulation box; these interactions are minimized when capsid subunits are adjacent to the surface of the sphere. Subunit positions and orientations are propagated according to overdamped Brownian dynamics, with the unit of time , where D is the subunit diffusion coefficient. Full details of the model are given in Ref. Hagan (2008)

Free energies for neutral subunits with no internal structure or high salt concentrations. The free energies for partial capsids on core surfaces are modified from the expressions given in Section III.3 for the case of neutral subunits or high salt concentrations, in which case interactions are limited to excluded volume, directional attractions that drive assembly, and favorable core-subunit interactions. The free energy of a core with adsorbed but unassembled subunits and an adsorbed intermediate of size is written as


where is the core-subunit surface free energy strength.

The second term, accounts for the entropy for two-dimensional motions of adsorbed subunits on the core surface. A simple approach would assume Langmuir adsorption, but we find that much more accurate results are obtained by integrating an empirical formula for the chemical potential of a fluid of hard disks Siepmann et al. (1992) because the Langmuir model overestimates subunit mixing entropy at high surface coverages:


with the packing fraction .

The last term in Eq. 18, , is the entropy penalty for subunits to be localized a smaller distance from the core surface than the subunit diameter, which we obtained by numerically integrating the partition function for an adsorbed subunit when comparing the kinetic theory to simulations.

Calculating the time dependence of adsorbed and assembled subunits requires choosing values for several theoretical parameters; here we describe the parameter values and how they were chosen.

  1. Subunit binding entropy penalty. Eq. 9 includes a term for the entropy loss (in addition to the subunit mixing entropy included in ) for a subunit to bind to a capsid that accounts for translational restrictions on scales smaller than the subunit diameter and rotational restrictions Hagan and Chandler (2006); the entropy penalty should increase only slightly when a subunit changes from one to more than one bond. The subunit-subunit binding entropy penalty is calculated as described in Ref. Hagan and Chandler (2006). There is a typo in Eq. (15) of that reference; the correct formula is

  2. Subunit binding rate constant . The assembly rate constant was chosen as (in the dimensionless units of Ref. Hagan (2008)) by comparing kinetic theory predictions to simulation results for the assembly of empty capsids. The resulting fit (not shown) shows less agreement between theory and simulation results than we find for core simulations, perhaps because approximations such as an average assembly pathway without multimer binding are less significant when subunits are confined to the surface of a core.

  3. The effective surface concentration of adsorbed subunits is increased if subunits are confined within less than a molecular diameter by the core-subunit interaction, so we take . This factor ensures that entropy losses due to subunit localization are not counted twice.

  4. Adsorption rate, . The rate of subunit adsorption to core surfaces in our simulations can be calculated from the Smoluchowski equation, with forces due to the adsorption potential explicitly included. Instead, we use the diffusion limited rate for the zero force case, but take , where is the subunit-core center of mass distance at the point when the interaction force becomes nonzero. This choice was validated by comparing theory and simulation results for the time dependence of the number of adsorbed subunits in the case of no assembly, i.e. . Since the excluded volume of adsorbed subunits is accounted for in Eq. 19 we set the adsorption blocking factor (see Eqs. 12). This approach slightly over estimates the adsorption rate at moderate coverage, but yields the proper equilibrium surface density if there is no assembly. Because of the definition of the unit time , the subunit diffusion constant is .

We find that the agreement in most cases is surprisingly good, considering the degree of error in estimating and , and the approximation that there is only one reaction pathway.

Appendix B

In this appendix we summarize the application of the method of Scheutjens and Fleer Scheutjens and Fleer (1979) to polyelectrolytes. Our presentation closely follows that given in Bohmer et al. Bohmer et al. (1990), except there is a different geometry, the Flory-Huggins parameter for interactions between polymer segments is set to , and there are additional terms in the free energy equations.

We consider a lattice around an impenetrable sphere with radius , contacting bulk solution in layer , with large enough that the presence of the surface is negligible in layer . All quantities are assumed uniform within a layer . Fixed numbers of polyelectrolytes and TEG molecules are grafted in layers and , respectively, with the lattice dimension nm.

In the calculation, a segment of species is subjected to a potential field , which is given by (relative to the bulk solution)


where the term represents the hard-core interaction, which is the same for every species and ensures that the total volume fraction for all species sums to unity in layer . The second term on the right-hand side accounts for electrostatic interactions, with the valency for a segment of species , the charge on an electron, and the electrostatic potential in layer . The statistical weight of a segment for species in layer , relative to the bulk solution, is given by the segment weighting factor


so that the volume fraction of a free monomer of type is given by .

The volume fractions for segments in chain molecules are calculated with a segment distribution function , which gives the statistical weight of a chain conformation starting with a segment 1, located anywhere in the lattice, and ending in layer after steps. The segment distribution function is calculated from a recurrence relation:


with , , and the fraction of contacts for a lattice site in layer with respective layers ,, and (see Ref. van Lent and Scheutjens (1989)). The statistical weight of a chain conformation that starts at segment and ends with segment in layer is calculated in the same way. The segment distribution functions for grafted molecules are modified to constrain segment to begin in the grafted layer as described in Ref. Cosgrove et al. (1987). To account for impenetrable surfaces, for for all species, and for for polyelectrolyte and TEG molecules.

The volume fraction of segment in a chain of segments is determined from the joint probability that a chain conformation starts at segment 1 and ends at segment in layer , and a chain conformation starts at segment and ends at segment in layer :


where the probabilities are divided by to avoid double counting of segment and is a normalization constant. If molecules of species are free to exchange with the bulk solution, is determined by the fact that all segment weighting factors in bulk so that


If the total amount of molecules of species is fixed (i.e. for a grafted brush) is determined from


with the total amount of molecules of species calculated from


and the chain weighting factor , which gives the statistical weight for a chain of type to be found anywhere in the lattice. The number of lattice sites depends on the layer because of curvature and is calculated in Ref. van Lent and Scheutjens (1989).

The electrostatic potential. Electrostatics are accounted for with a multi-Stern-layer model in which all charges within a lattice layer are located on the plane at the center of the layer, with no charge outside of that plane. The plane charge density in layer , is calculated from a sum over all species


where we set the cross-sectional area per lattice site with the distance between layers. For weak acid groups the degree of dissociation is given byBohmer et al. (1991)


with the dimensionless dissociation parameter related to the dissociation constant by with the molarity of pure segments for species .

The electrostatic potential with respect to bulk solution at layer is calculated through electroneutrality Bohmer et al. (1990) and the potential in each layer is calculated from Gauss’ law, for spherical layers it is Bohmer et al. (1991)


with the volume fraction averaged dielectric permittivity for layer .

Free energies. The free energy is obtained from the canonical partition function for the lattice system Scheutjens and Fleer (1979); Evers et al. (1990); Bohmer et al. (1990)


with the electrostatic energy for the case of a fixed surface charge given by


The final term in Eq. 31 was omitted in Ref. Bohmer et al. (1990), but is necessary to describe the free energy due to dissociation of acid groups