A theory for viral capsid assembly around electrostatic cores
aa r X i v : . [ q - b i o . B M ] M a y A theory for viral capsid assembly around electrostatic cores
Michael F. Hagan
Department of Physics, Brandeis University, Waltham, MA, 02454 (Dated: November 11, 2018)We develop equilibrium and kinetic theories that describe the assembly of viral capsid proteins on a chargedcentral core, as seen in recent experiments in which brome mosaic virus (BMV) capsids assemble aroundnanoparticles functionalized with polyelectrolyte. We model interactions between capsid proteins and nanopar-ticle surfaces as the interaction of polyelectrolyte brushes with opposite charge, using the nonlinear PoissonBoltzmann equation. The models predict that there is a threshold density of functionalized charge, above whichcapsids 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 en-able interpretation of, preliminary experimental data. However, the models predict a stronger dependence ofnanoparticle incorporation efficiency on functionalized charge density than measured in experiments, and donot completely capture a logarithmic growth phase seen in experimental light scatter. These discrepancies maysuggest the presence of metastable disordered states in the experimental system. In addition to discussing fu-ture experiments for nanoparticle-capsid systems, we discuss broader implications for understanding assemblyaround charged cores such as nucleic acids.
I. INTRODUCTION
The cooperative assembly of heterogenius building blocksinto ordered structures is crucial for many biological pro-cesses, such as the assembly of protein subunits around a vi-ral nucleic acid to form a capsid. Understanding how nucleicacid-protein interactions drive cooperative assembly could en-able antiviral drugs that block this essential step in viral repli-cation. At the same time, engineered structures in whichviral capsids assemble around synthetic cargoes show greatpromise as delivery vehicles for drugs [1, 2, 3] or imagingagents [4, 5, 6, 7] and as subunits or templates for the synthesisof advanced multicomponent nanomaterials [8, 9, 10, 11, 12].Realizing these goals, however, requires understanding at afundamental level how properties of cargoes and capsid pro-teins control assembly rates and mechanisms. With that goalin mind, this article presents coarse-grained thermodynamicand kinetic models that describe experiments in which bromemosaic virus (BMV) capsid proteins assemble around spheri-cal charge-functionalized nanoparticle cores.
Identifying mechanisms for simultaneous assembly andcargo-encapsidation.
Elegant experiments have studied theassembly of capsids around central cores consisting of nu-cleic acids (e.g. [13, 14, 15, 16, 17, 18, 19, 20], inorganicpolyelectrolytes [21, 22, 23], and charge-functionalized cores[7, 24, 25, 26, 27, 28]. However, assembly mechanisms ofviruses and virus-like particles remain incompletely under-stood because cargo-protein interactions and the structures oftransient assembly intermediates are not accessible to experi-ments. For the case of empty-capsid assembly, Zlotnick and coworkers achieved great insight about assembly mechanismswith simplified theories (e.g. [29, 30, 31, 32, 33, 34, 35]),for which a small number of parameters are fit by compari-son with experimental assembly data over a range of proteinconcentrations and pH values. The theories developed in thiswork enable similar comparisons for the assembly of capsidsaround electrostatic cores, and could thereby elucidate mech-anisms for simultaneous assembly and cargo encapsidation.Prior theoretical works have led to important insights aboutassembly around polymeric cores, but were equilibrium stud-ies [36, 37, 38, 39, 40 ? ] or considered specific RNA-mediated assembly pathways with phenomenological descrip-tions of protein-nucleic acid interactions [41, 42]. Computersimulations of a model in which subunits assemble aroundrigid spherical cores [43, 44] suggest that a core with a geom-etry commensurate with the lowest free energy empty-capsidmorphology can increase assembly rates and the efficiency ofassembly (compared to empty capsid assembly[45, 46, 47, 48,49, 50, 51, 52, 53, 54]) through the effects of heterogeneousnucleation and templating; however, assembly is frustrated forcore-subunit interaction strengths that are too strong[43, 55].Cores that are not size-matched with the lowest free energycapsid morphology can direct assembly into alternative capsidmorphologies for optimal core-subunit interaction strengths[43]. Important limitations of the simulations are that ex-tensive comparison with experimental data is computationallyintractable, and the simulated models do not explicitly rep-resent electrostatics or the polymeric nature of disordered N-terminal tails in capsid proteins.The present work aims to overcome these limitationsby calculating the electrostatic interaction between charge-functionalized nanoparticles and capsid proteins, and thenconstructing simplified thermodynamic and kinetic theoriesthat describe the simultaneous assembly and cargo encapsi-dation process. We explore the effect of cargo-capsid proteininteractions on assembly mechanisms by presenting predictedassembly rates and nanoparticle incorporation efficiencies asfunctions of the functionalized charge density, nanoparticle-protein stoichiometric ratio, and the capsid protein concentra-tion. Because each of these parameters can be experimentallycontrolled, the theoretical predictions suggest a series of ex-periments that can elucidate the relationship between cargo-capsid protein interactions and assembly mechanisms.Although there is currently not sufficient experimental datato estimate values for several unknown parameters (adsorptionand assembly rate constants, and subunit-subunit binding en-ergies), we qualitatively compare the theoretical predictions topreliminary experimental time-resolved light scattering data(Fig. 1a) and measurements of the efficiency of nanoparticleincorporation in well-formed capsids (Fig. 1b). With physi-cally reasonable values for the unknown parameters, the pre-dicted light scatter increases rapidly at short times without alag phase, consistent with the experimental measurements, butthe theory does not completely capture the functional form oflight scatter at long times. The theory predicts that nanoparti-cles are incorporated into capsids only above a threshold func-tionalized charge density, which is similar to the lowest chargedensity for which significant incorporation is observed in ex-periments. However, predicted incorporation efficiencies riseto 100% over a narrow range of charge densities, whereas theexperimental data shows a gradual increase in the incorpo-ration efficiency with increasing charge density. We discussexperiments and potential improvements to the theory whichmight elucidate these discrepancies.This paper is organized as follows. In Section II we calcu-late the interaction between capsid proteins and functionalizedcharge on nanoparticle surfaces. In Section III, we presenta model for the thermodynamics of capsid assembly in thepresence of rigid spherical cores, followed in Section IV bya kinetic theory that describes the simultaneous assembly ofcapsid subunits on nanoparticle surfaces and in bulk solution.We analyze kinetics and the assembly pathways predicted bythe model in section V, and discuss the connection to currentand future experiments in Section VI. We also compare theorypredictions to simulation results in Appendix A. L i gh t sc a tt e r [ A . U .] t [sec] (a) P a ck ag i ng E ff i c i en cy σ C [acid groups / nm ] (b) FIG. 1: (a)
Time-resolved light scatter for capsid assembly onnanoparticles with a functionalized surface density of σ C = 3 acidgroups/nm . (b) Packaging efficiency, or the fraction of nanopar-ticles incorporated into T =3 capsids (measured from TEM micro-gaphs), varies with σ C . The capsid protein-nanoparticle stoichiomet-ric ratio is r CP:NP = 1 . , core diameters are 12 nm (which pro-mote T =3 capsid formation), the capsid protein concentration is 0.4mg/ml ( C p = 10 µ M dimer subunit), and there is 100 mM 1:1 saltand 5 mM 2:1 salt. Packaging efficiencies were measured at t = 10 seconds. The surface density of acid groups is estimated from thefraction of carboxylated TEG molecules, assuming a total surfacedensity of 3 TEG molecules/nm on the nanoparticle surface. Datathanks to Bogdan Dragnea [56]. Experimental motivation.
This work is motivated by ex-periments in which BMV capsid proteins assemble aroundgold nanoparticle cores [24, 27] functionalized with thio-lalkylated tetraethylene glycol (TEG) molecules, a fraction ofwhich are terminated with carboxylate groups [27]. With 12nm nanoparticles the TEG-nanoparticle complexes have a di-ameter of roughly 16 nm, which is nearly commensurate withthe interior of a T=3 BMV capsid (17-18 nm [57]). Ionizedcarboxyl groups provide a driving force for capsid proteins to
FIG. 2: The model system. Surface functionalization molecules(TEG) are end-grafted to an impenetrable gold sphere. The cationiccapsid protein N-terminal arms are modeled as polyelectrolytesgrafted to the inner surface of a sphere (the capsid), which is im-permeable to TEG and polyelectrolyte but permeable to solvent andions. adsorb and assemble on nanoparticle surfaces. Time-resolvedlight scatter measurements (Fig. 1a), which monitor assem-bly kinetics, are strikingly different than light scatter of emptycapsid assembly at short times (e.g. [58, 59, 60]); instead ofa lag phase ( ∼ seconds – minutes) there is a rapid ( ∼ second) rise in signal. As described in Section V, the theoret-ical predictions in this work begin to explain this observation.Assembly effectiveness is monitored by using TEM to mea-sure packaging efficiencies, or the fraction of nanoparticles in-corporated within well-formed capsids. As shown in Fig. 1b,packaging efficiencies increase from 0% to nearly 85% as thefraction of TEG molecules that are carboxylate-terminated in-creases from 20% to 100%.While the light scattering data was obtained from an as-sembly reaction entirely at p H = 5, the packaging efficiencydata shown in Fig. 1a were obtained with a protocol in whichcapsid subunits were first dialyzed against a buffer at p H =7 . , followed by buffer at p H = 5 . In this work we con-sider only p H = 5 , since at p H = 7 . disordered protein-nanoparticle complexes form but well-formed capsids are notobserved without subsequently lowering p H. Assembly maynot be favorable at p H = 7 . because binding interactions be-tween capsid proteins are less favorable than at lower p H [33];also strong TEG-protein interactions might render disorderedstates kinetically or thermodynamically favorable.
II. SUBUNIT ADSORPTION ON NANOPARTICLESURFACES
For many capsid proteins, there is a high density of basicresidues on flexible N-terminal tails, with numbers of chargedamino acids that range from six to 30 [36]. These positivecharges help drive capsid assembly in vivo by interacting with v o l u m e f r a c t i on ∆ z (a) σ C =-3 σ C =-2 σ C =-1 0 0.01 0.02 0.03 0.04 0.05 0.06 1 2 3 4 5 6 7 8 9 v o l u m e f r a c t i on ∆ z (b) σ C =-3 σ C =-2 σ C =-1 FIG. 3: The volume fractions of (a) polyelectrolyte segments and (b) carboxylate groups in the grafted TEG layer are shown as a func-tion of layer ( ∆ z = z − R C /l + 1 ) above the nanoparticle surfacefor several surface acid group densities, with the nanoparticle radius R C = 6 nm and the layer dimension l = 0 . nm. negative charges on nucleic acids, and the interaction of poly-meric tails and nucleic acids is explored in Refs. [36, 42, 61].In this section, we analyze the interactions between capsidprotein tails and surface functionalized TEG molecules. Forsimplicity, we neglect other effects, such as hydrophobic in-teractions, that might drive subunit adsorption, and we assumethat all charges on capsid proteins are located in the polymerictails. For the remainder of this work, we will reserve the word‘polyelectrolyte’ to apply to the charged polymeric tails oncapsid proteins; we will not use this word to refer to the func-tionalized TEG molecules. Model system.
We first consider a dilute solution of cap-sid protein subunits with density ρ p , and a single nanoparticle;we consider a finite concentration of cores in Section III. Athigh surface coverage or in an assembled capsid the foldedportion of capsid proteins acts as an impenetrable barrier tothe polymeric tails. We therefore represent the nanoparticle-capsid complex as a concave polyelectrolyte brush (the inte-rior surface of the capsid) interacting with a convex brush ofthe opposite charge (the TEG layer); the model system geom-etry is illustrated in Fig. 2. We consider surface functional-ization molecules grafted to the outer surface of an impene-trable sphere with radius 6 nm. The TEG surface function-alization molecules used in Ref. [24] consist of a 10 carbonalkanethiol chemically linked to a polyethylene glycol chainof 4 monomer units. The surface density of functionalizedcharge groups, σ C , is experimentally controlled by terminat-ing a fraction TEG molecules with carboxylates, with the re-mainder terminated with neutral groups [56]. We representTEG molecules as freely jointed chains with 9 neutral seg-ments 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 proteinscontain nine positive charges in the first 20 amino acids andthe first 25-40 residues are disordered in capsid crystal struc-tures [57]. Since BMV capsids assemble from protein dimersubunits [62, 63, 64], in our model each subunit has a netcharge of q p = 18 . We represent the two N-terminal armsfrom each dimer subunit as a single polymer of 36 segmentswith Kuhn length 3.5 ˚ A , and valence 0.5. These polyelec-trolytes are end-grafted to the inner surface of a sphere withradius 9 nm, the inner diameter of a T=3 BMV capsid [57].The model capsid is permeable to solvent and ions, but im-penetrable to TEG and polyelectrolytes segments. The lat-ter condition is reasonable for high surface coverages of cap-sid subunits and in the strong coupling limit, in which casethe polyelectrolyte configurational entropy is negligible com-pared to electrostatic effects. For low charge densities, theresults do not change qualitatively if the impermeable capsidis eliminated, but the equilibrium concentration of adsorbedpolyelectrolytes increases by up to 50%. Numerical calculations.
To our knowledge, interactingpolyelectrolyte brushes with this geometry have not been stud-ied, although the interaction of capsid arms with a nucleicacid is considered in Refs. [36, 65]. For many cases we con-sider, the electrostatic potential in the region of the capsid andnanoparticle surface is large; therefore, we will solve the non-linear Poisson Boltzmann (PB) equation. To do so, we employthe method of Scheutjens and Fleer (SF) [66, 67, 68], in whichthe spatial distributions of polyelectrolyte, TEG, ions, and wa- ter, as well as the dissociation equilibrium for carboxylate andwater groups, are solved numerically with a self consistentfield approximation on a lattice. The calculation accounts forthe finite size of all species and the calculated free energiesexplicitly consider the entropy of ions and solvent molecules.The method is thoroughly described in Refs. [67, 69]; our im-plementation and additional terms that must be included in thefree energy to account for dissociation of weak acid groups aredescribed in Appendix B. We neglect spatial variations of seg-ment densities in the directions lateral to the surface (and thusneglect the effect of ion-ion correlations), and determine thevariation of densities in the direction normal to the surface.All calculations consider the conditions used for the exper-imental data shown in Fig. 1a, with pH=5, and 100 mM 1:1salt with an additional 5mM 2:1 salt. The total surface den-sity of functionalized molecules (neutral and charged) is ∼ nm − at the nanoparticle surface[70]. Rather than explicitlymodeling two species of surface molecules (neutral and acid-terminated), we vary the surface charge density σ C by chang-ing the valency of the ionic end group. To explore a widerange of surface charge we consider end groups with valen-cies v ∈ [0 , ; the experimental data in Fig. 1 correspondsto the range v ∈ [0 , . For v = 1 the number of function-alized charge groups corresponds to ∼ of the charge ona BMV capsid. The lattice size and statistical segment lengthfor polyelectrolyte and TEG molecules [71, 72, 73] are set to l = 0 . nm (which roughly corresponds to the size of a singleamino acid and the correlation length of water), and there are0.5 charges per polyelectrolyte segment, which correspondsto the linear charge density predicted by counterion condensa-tion [74, 75, 76, 77]. Except for charge, all species are treatedas chemically identical, with the Flory interaction parameter χ = 0 and dielectric permittivity ǫ = 80 ). Free energiesare insensitive to the introduction of unfavorable interactions( χ > ) between the alkane segments and hydrophilic speciesor species-dependent dielectric constants. In the strong cou-pling limit, the equilibrium density of positive charge due toadsorbed subunits is not sensitive to statistical segment length,lattice size, polymer length or charge per segments. The disso-ciation constant ( pK a ) for carboxylate groups is not known inthe vicinity of the surface functionalized gold particle; we use pK a = 4 . , so that approximately 25% of the acid groups aredissociated for p H = 5 at full coverage (see below), consistentwith Ref. [78]. We do not explicity consider complexation ofacid groups, a possibility suggested in that reference. We will -1.5-1-0.5 0 0.5 1 1.5 2 2.5 0 0.2 0.4 0.6 0.8 1 f S F [ k B T / n m ] ρ surf pKa 4.5, σ C =3pKa 4.5, σ C =1pKa 0.0, σ C =3pKa 0.0, σ C =1 FIG. 4: The free energy per area at the capsid surface for the adsorp-tion of a concave polymer brush onto a spherical brush with oppositecharge. For varying surface densities of polyelectrolyte (normalizedby the density in a complete capsid), ρ surf , the free energy relative to ρ surf = 0 is shown for nanoparticle functionalization with weak acidgroups ( pK a = 4 . ) at σ C = 3 acid groups/nm ( (cid:4) ), σ C = 1 ( N )and strong acid groups at σ C = 3 (+), σ C = 1 (o). There is 100 mM1:1 and 5 mM 2:1 salt and p H = 5 . ρ eq s u r f σ C [acid groups / nm ]pKa 4.5pKa 0.0 FIG. 5: The equilibrium surface density of adsorbed polyelectrolytes, ρ eqsurf , normalized by the surface density of polyelectrolyte on a corewith a complete capsid. Predictions are shown for surface function-alization with weak acid ( (cid:4) ) and strong acid (+) groups. The bulksubunit concentration is 10 µ M, with other parameters as in Fig. 4. also consider experiments in which the terminal group is astrong acid, with pK a = 0 . . Free energies and equilibrium concentrations of adsorbedpolyelectrolytes.
The equilibrium surface density of adsorbedpolyelectrolyte can be calculated by considering a system forwhich there is free exchange between surface adsorbed poly-electrolytes and a bath containing polyelectrolyte, salt, andsolvent. The kinetic theory developed in Section IV, however, requires the free energy for non-equilibrium adsorption den-sities. We therefore consider a restricted number of graftedmolecules (which cannot exchange with the bath), as a func-tion of surface charge and density of grafted molecules. Foreach grafting density ρ surf , the spatial distributions of poly-electrolyte, TEG, ions, and solvent are determined by mini-mizing the total free energy. The free energy, f SF , is calculatedfrom Eqs. B11 - B17. Note that the density of polyelectrolytes(Fig. 3) segments is nonmonotonic with distance from the cap-sid, consistent with studies of capsid assembly around nucleicacids [36, 65].Free energies for two functionalized charge densities withweak and strong acid groups are shown in Fig. 4 and theequilibrium adsorption densities, ρ eqsurf , which correspond tothe minimum excess free energy, are shown in Fig. 5. Theadsorption densities ( ρ surf , ρ eqsurf ) are normalized with respectto the density of dimer subunits in a complete capsid (0.09subunits/nm on the interior surface of the capsid with innerradius 8.9 nm). The predicted adsorption densities for weakacids ( pK a = 4 . ) are significantly lower than for the case ofstrong acids because the high local density of negative chargesshifts the carboxylate dissociation equilibrium; for a carboxy-late surface density σ C = 3 / nm the average dissociated frac-tion shifts from 75% for an isolated molecule to 35%. How-ever, dissociation is enhanced as polyelectrolytes with positivecharges adsorb and the average dissociation fraction on a corewith a complete capsid is 78%. As discussed in Section V,this effect results in significantly different dependencies of as-sembly kinetics on ρ eqsurf for weak and strong acids.There are several reasons why our calculations may un-derestimate ρ eqsurf at low σ C and for weak acid groups.Nanoparticle-protein interactions in addition to those involv-ing the N-terminal arms could contribute to subunit adsorp-tion. The calculated interactions involving N-terminal armsassume uniform charge densities in directions parallel to thesurface, which will overestimate polyelectrolyte interactionsat low surface densities, when the average distance betweenadsorbed subunits is greater than the polyelectrolyte radius ofgyration; the approach described for neutral subunits in Ap-pendix A might be appropriate in this regime. As discussedabove, reducing restrictions on polyelectrolyte conformationsdue to the folded portion of capsid proteins results in some-what higher adsorption densities for low surface charges. Inaddition, the pK a of carboxylate groups and the effect of in-teractions between the charged species and the gold surface -70-60-50-40-30-20-10 0 10 0 0.5 1 1.5 2 2.5 3 s u r f a c e po t en t i a l [ m V o l t] σ C [acid groups / nm ]pKa 4.5pKa 0.0 FIG. 6: The electrostatic potential at the nanoparticle surface forvarying surface functionalization ( σ C ) with strong (+) and weak ( (cid:4) )acid groups. There is 100 mM 1:1 and 5 mM 2:1 salt. are poorly understood. For this reason, we examine assemblyover a wide range of surface charges for both pK a = 4 . and pK a = 0 to understand the interplay between surface charge,dissociation, and assembly. It would be desirable to experi-mentally measure the surface charge as a function of the frac-tion of charged TEG molecules on the surface; while it willbe difficult to quantitatively verify the surface charge or cor-responding magnitude of the electrostatic potential, the vari-ation of the potential at the core surface with σ C (see Fig. 6)is roughly consistent with preliminary experimental measure-ments of zeta potentials for functionalized nanoparticles [70].While our calculations may underestimate ρ eqsurf , these approxi-mations are less important at high coverages and hence shouldnot significantly affect equilibrium packaging efficiency cal-culations. A. Empty-capsid free energies.
We use the same approach, but without the core and TEGmolecules, to estimate the electrostatic contribution to the freeenergy of forming an empty capsid. The free energy, f SF,EC , asa function of polyelectrolyte density ( ρ surf ) is shown in Fig. 7,where it is compared to the free energy calculated with thelinearized Poisson Boltzmann equation for charge smearedon a spherical surface [35, 61] with radius R cap = 8 . nm, f PB = λ B λ D q s n s / [2 R cap ( λ D + R cap )] , with n s = 90 subunits, q s =18 charges/subunit, λ B ≈ . nm the Bjerrum length ofwater, and λ D the Debye length. Although the linearized Pois- f E C , f PB [ k B T / n m ] ρ surf FIG. 7: The free energy per area at the capsid surface of an emptycapsid due to electrostatic repulsions and ion and solvent entropyis shown as a function of the charge density on capsid subunits, aspredicted by SF calculations (solid lines with symbols) and the lin-earized PB equation (dashed line). The calculations use 100 mM 1:1salt and 5 mM 2:1 salt. son Boltzmann equation was shown to agree closely with thefull nonlinear calculation in Ref. [61] for this ionic strength,there is a large discrepancy between the PB and SF calcu-lations; the SF calculations predict much lower free ener-gies because the flexible tails adopt stretched conformationswhich reduce electrostatic repulsions (see also Ref. [65]). Wenote, though, that we have not considered interactions be-tween charges located in the structured region of capsid pro-teins.Zlotnick and coworkers [29, 33, 79] have shown that capsidformation free energies can be fit using the law of mass actionto experimental data for capsid and free subunit concentra-tions as the total capsid protein concentration is varied. Theresulting free energies depend on p H as well as salt concentra-tion [33] and likely include effects arising from hydrophobicas well as electrostatic interactions [35]. We can subtractthe polyelectrolyte free energy from the total capsid forma-tion free energy derived from experimental data to estimatethe free energy due to attractive subunit-subunit interactions, g H , which may arise from hydrophobic contacts. FollowingCeres and coworkers [33], we assume that each subunit in acapsid makes favorable contacts with n c = 4 neighboring sub-units, so the attractive energy per subunit is g H = g b n c / − πR cap f SF,EC /N (1)For N = 90 protein dimer subunits in a T =3 BMV cap-sid with an inner capsid radius of R cap = 8 . nm [57], weobtain f SF,EC = 2 . k B T / nm . With a typical subunit-subunit contact free energy of g b = − k B T (-3 kcal/mol)[29, 33, 79], this gives g H ≈ − k B T / subunit (compare tothe estimate from experimental data for hepatitis B virus as-sembly by Kegel and van der Schoot [35] of an attractive in-teraction strength of − k B T per subunit). III. EQUILIBRIUM THEORY FOR CORE-CONTROLLEDASSEMBLY
In this section we consider the thermodynamics for a sys-tem of subunits that can assemble into capsids around coresand assemble in bulk solution to form empty capsids. Thefree energies due to electrostatic interactions for subunits ad-sorbed on surfaces or assembled into empty capsids developedin section II are used to determine the relative free energies ofcore-associated capsid intermediates.
A. The thermodynamics of core-controlled assembly
We consider a dilute solution of capsid protein subunitswith density ρ p , and solid cores with density ρ C . Subunits canassociate to form capsid intermediates in bulk solution or oncore surfaces. A complete capsid is comprised of N subunits.Zandi and van der Schoot [ ? ] develop the free energy forcore-controlled assembly in the context of capsids assemblingaround polymers, but neglect the contribution of partial cap-sid intermediates, which will be important for describing thekinetics 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 thatdisordered configurations could be kinetically arrested or eventhe favored equilibrium state in the case of core-subunit inter-actions that are much stronger than the thermal energy k B T [43, 55].The total free energy density is given by F = F EC + F core (2)with F EC the free energy density of empty capsid intermedi-ates F EC = N X i =1 ρ i log( ρ i a ) − ρ i + ρ i G cap i (3) where ρ i is the density of intermediates with i subunits and a is a standard state volume. For simplicity we follow Zlot-nick and coworkers [30, 31] and assume that there is only oneintermediate for each size i , with a free energy G cap i definedbelow in Eq. 9.The free energy density of core-associated intermediates isgiven by F core = ρ C (cid:20) N X n =0 N X m =0 ′ P n,m log( P n,m ρ C a ) − P n,m + P n,m G core n,m (cid:21) (4)where P n,m is the fraction of cores with an adsorbed inter-mediate of size m and with n adsorbed but unassembled sub-units. The free energy for such a core is given by G core n,m andis defined below. For simplicity, we assume that a core canhave at most one assembled intermediate and a total of N sub-units can fit on a single core surface, where N is the size ofa complete capsid, so states with m + n > N subunits havezero probability. The prime on the second sum indicates that m = 0 , , . . . N because, by definition, a state with m = 1 has no assembled intermediates.To obtain the equilibrium concentration of intermediates weminimize the total free energy subject to the constraints ρ C N X n =1 N X m =0 ′ ( n + m ) P n,m + N X i =1 iρ i = ρ p (5)and N X n =1 N X m =0 ′ P n,m = 1 (6)where Eq. 5 enforces that the total density of subunits is givenby ρ p .Minimization yields ρ i a = exp[ − β ( G cap i − iµ )] (7)and P n,m ρ C a = exp[ − β ( G core n,m − ( n + m ) µ − µ C )] (8)where β = 1 /k B T is the inverse of the thermal energy, µ = k B T log( a ρ ) and µ C = k B T log( a ρ C P , ) are chemi-cal potentials for subunits and cores, respectively. B. Free energies for empty-capsid assembly intermediates
The free energy of an empty capsid intermediate with i sub-units is written as G cap i ( g b ) = i X j =1 ( n c j g b ) − k B T S degen i − ik B T s rot (9)where n c j is the number of new subunit-subunit contactsformed by the binding of subunit j to the intermediate, g b is the contact free energy, S degen accounts for degeneracy inthe number of ways subunits can bind to or unbind from anintermediate (see the s factors in Refs. [30, 31]), and s rot isa rotational binding entropy penalty. For simplicity and toease comparison with earlier works, we take s rot = 0 exceptwhen comparing theoretical predictions to simulation results(see Appendix A). As discussed at the end of Section II, thesubunit-subunit contact free energy has been estimated fromexperiments for several viruses [29, 33, 59, 79] and dependson salt concentration and p H [33, 35].
C. Free energies for assembly intermediates on core surfaces
The free energy for a core without a partial capsid is givenby G core n, = A C f SF ( σ C , ρ surf ) with A C the core area, σ C thefunctionalized acid group density, ρ surf = n/N the normal-ized density of adsorbed subunits, and f SF calculated in Sec-tion II. Likewise, the free energy for a complete capsid is G core = G cap m + A C [ f SF ( σ C , − f SF,EC ] , where the polyelec-trolyte contribution to the capsid free energy ( f SF,EC ) is sub-tracted because it is included in f SF .In the case of a core with a partially assembled capsid andadsorbed but unassembled subunits, there is an inhomoge-neous distribution of charge on the surface, with higher chargedensity in the region of the assembled cluster. We determinethe free energy for such a core, with n adsorbed subunits andan intermediate of size m , in the following schematic way.We assume the unassembled subunits spread on the surfaceof the core not occupied by the intermediate. There are thustwo ‘pools’ of charge on the surface – the partial capsid witha high subunit density ρ surf = 1 and the remaining area of thecore surface with a charge density spread uniformly over thearea not covered by the intermediate , A free m = A C ( N − m ) /N .The total electrostatic free energy is then given by G surf n,m = G cap m + (cid:0) A C − A free m (cid:1) [ f SF ( σ C , − f SF,EC ] + A free m f SF ( σ C , nq p N − m ) for n + m ≤ N ∞ for n + m > N . (10)In general there should be a term in G core to account for strainenergy if the core curvature is not commensurate with that ofthe intermediate; we assume that capsid and core curvaturesmatch. Core encapsidation.
Above a threshold or critical sub-unit concentration (CSC), the majority of subunits will be incapsids[30, 45, 80]. Because core-subunit interactions canlead to high localized surface charge densities, the thresh-old concentration for assembly in the nanoparticle systemcan be below the CSC for spontaneous assembly. Abovethe 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( G core N, − G cap N ≈ − k B T log( ρ C a ) for a stoichiometric amountof capsid protein and nanoparticles), below which empty cap-sids and free cores are favored over full capsids. Above thethreshold surface free energy and above the CSC, nearly allcores will be encapsidated ( P ,N ≈ ) at equilibrium. How-ever, we will see that there can be a higher threshold surfacefree energy for efficient encapsidation kinetics. IV. KINETIC THEORY FOR CORE-CONTROLLEDASSEMBLY
In this section we use the free energies developed in sec-tion III for capsids and capsid intermediates as the basis for akinetic theory of core-controlled assembly.Zlotnick and coworkers [30, 31] have developed a systemof rate equations that describe the time evolution of concen-trations of empty capsid intermediates dρ dt = − k as s ρ − N X i =2 ˆ k as s i +1 ρ i ρ + ˆ k b ρ i dρ i dt = ˆ k as s i ρ ρ i − − ˆ k as s i +1 ρ ρ i i = 2 . . . N − ˆ k b ρ i + ˆ k b ρ i +1 (11)where ˆ k as and ˆ k b are respectively the forward and reversebinding rates and s i is a statistical factor that describes thenumber of different ways to transition from intermediate i to i + 1 [30, 31]; s corrects for double counting the number ofpairwise combinations of free subunits. Following Ref. [31],we simplify the calculation by requiring that transitions be-tween intermediates are only allowed through binding or un-binding of a single subunit and there is only one intermediatefor each size i . We assume that the assembly rate constant ˆ k as is independent of intermediate size i , but the number ofcontacts per subunit is not (see Table 1).We extend this approach to describe simultaneous capsidassembly and adsorption to spherical cores, giving the follow-ing expression for the time evolution of core states dP n, dt = ρ Ω n − , P n − , + Ω n +1 , P n +1 , − (cid:0) ρ Ω n, + Ω n, (cid:1) P n, + W n − , P n − , + W spont P n, − ( W n, + ρ W spont ) P n, dP n, dt = Ω n − , ρ P n − , + Ω n +1 , P n +1 , − (cid:0) ρ Ω n, + Ω n, (cid:1) P n, + W n +2 , P n +2 , + W n − , P n − , + W spont P n, + W spont P n, − (cid:16) W n, + W n, + ρ W spont + W spont (cid:17) P n, dP n,m dt = ρ Ω n − ,m P n − ,m m = 3 . . . N +Ω n +1 ,m P n +1 ,m − (cid:0) ρ Ω n,m + Ω n,m (cid:1) P n,m + W n +1 ,m − P n +1 ,m − + W n − ,m +1 P n − ,m +1 + ρ W spont m − P n,m − + W spont m +1 P n,m +1 − (cid:16) W n,m + W n,m + ρ W spont m + W spont m (cid:17) P n,m (12)The first two lines in each of Eqs. 12 describe adsorptionand desorption from the core surface with respective rate con-stants Ω and Ω , while the following lines describes bindingand unbinding to a capsid intermediate. Rate constants forbinding and unbinding of adsorbed subunits to an adsorbedintermediate are denoted by W and W , respectively. Finally, W spont and W spont are respectively the rate constants for theassociation (dissociation) of non-adsorbed subunits to (from)an adsorbed intermediate. The latter process becomes im-portant when an adsorbed capsid nears completion, and elec-trostatic repulsions render non-specific adsorption of subunitsunfavorable, causing the surface capture and diffusion mech-anism to be slow in comparison to the empty-capsid assemblymechanism. Adsorption is calculated from Ω n,m = k ad Φ n,m , with k ad the adsorption rate constant and Φ n,m the probability that anadsorbing subunit is not blocked by previously adsorbed sub-units. For simplicity, we assume Langmuir kinetics, Φ n,m =( N − m − n ) /N , but note that this is a poor approximation athigh surface coverages (see Appendix A).The desorbtion rate Ω n,m is related to the adsorption rateby detailed balance Ω n,m = Ω n − ,m exp (cid:0) G core n,m − G core n − ,m (cid:1) /a . (13)Assembly and disassembly of adsorbed subunits is de-scribed in a manner analogous to that used for empty capsids,with assembly rates given by W n, = k as s ( n − ρ surf n, (14) W n,m = k as s n +1 ρ surf n,m D C m = 2 . . . N − where the effective surface density is ρ surf n,m = n/ [ a ( N − m )] .Eq. 15 states that subunit-subunit binding rates are propor-tional to the frequency of subunit-subunit collisions, whichcan be calculated from the Smoluchowski equation for thediffusion limited rate in three dimensions, with the density ofsubunits ( ρ surf n,m ) within a layer above the surface with thick-ness of the subunit size, a , or from the diffusion limitedrate in two dimensions [81] with a surface density given by ρ = ρ surf n,m a . The result of the two calculations differs onlyby a logarithmic factor that is of order 1 in this case (seeappendix B of Ref. [81]). Surface assembly rate constants( k as ) can be smaller than those for empty capsid assembly( ˆ k as ), if interactions between subunits and the TEG layer im-pede diffusion or reorientation of adsorbed subunits. Associ-ation rates could also change in comparison to bulk solutionif subunits take different conformations on the core surface;however, imaging experiments show that protein structuresin nanoparticle-filled capsids are nearly identical to those inempty capsids [24], which suggests that subunits do not dena-ture on the surface.The disassembly rate is given by detailed balance W n,m = W n +1 ,m − exp (cid:0) G core n,m − G core n +1 ,m − (cid:1) /a . (15)Similarly the association and dissociation rates for solubi-lized subunits are related by detailed balance: W spont = 2 k as s ρ W spont m = k as s m +1 m = 2 . . . N − W spont m = W spont m exp (cid:0) G core n,m − G core n,m − (cid:1) /a (16)0The expression for the time evolution of empty capsid in-termediates is now dρ dt = − k as s ρ − N X i =2 ˆ k as s i +1 ρ i ρ + ˆ k b ρ i + ρ C N X n =1 N X m =1 (cid:16) Ω n,m − Ω n,m + W spont m − ρ W spont m (cid:17) P n,m dρ i dt = ˆ k as s i ρ ρ i − − ˆ k as s i +1 ρ ρ i i = 2 ..N − ˆ k b u i + ˆ k b ρ i +1 . (17)We note that although Eqs. 12 and 17 have the form of aMaster equation, the finite concentration of subunits intro-duces a nonlinear dependence of the time evolution on stateprobabilities. V. RESULTS
In this section we report the predictions of the equilib-rium and kinetic theories for two experimentally observablequantities: packaging efficiencies, or the fraction of nanopar-ticles encapsulated by well-formed capsids (estimated fromTEM experiments), and time-resolved light scattering inten-sity. Preliminary experimental data for these quantities isshown in Fig. 1. Since there is not sufficient data to esti-mate the values of all unknown parameters, we make qualita-tive predictions about the effect of changing the experimentalcontrol parameters: the functionalized surface charge density σ C , the stoichiometric ratio of capsid protein to nanoparticles r CP:NP , and the capsid protein concentration C C . The impor-tant unknown parameters are the subunit-subunit binding en-ergy, g b , the dissociation constant of the functionalized chargegroups, the rate constant for adsorption of subunits onto thenanoparticle surface, k ad , and the rate constant for the assem-bly of surface-adsorbed subunits, k as . The definitions and val-ues of model parameters are given in Table I. A. Packaging efficiencies
We begin by discussing the relationship between the densityof functionalized charge groups on nanoparticle surfaces, σ C ,and packaging efficiencies. As discussed in Section I, the ex-perimental measurements of packaging efficiencies were ob-tained from an experimental protocol in which the solution p Hvaried over time; because we do not know the dependence of P a ck ag i ng E ff i c i en cy σ C [acid groups / nm ] (a) g b =-3 g b =-2 g b =-1 g b =0 0 0.2 0.4 0.6 0.8 1 0.5 1 1.5 2 2.5 P a ck ag i ng E ff i c i en cy σ C [acid groups / nm ] (b) g b =-5 g b =-3 g b =-2 g b =-1 FIG. 8: Packaging efficiencies, or the fraction of nanoparticles incor-porated in well-formed capsids, at varying surface charge densities( σ C ) of weak acid groups predicted by (a) the equilibrium theoryand (b) the kinetic theory at seconds. Parameters are as listedin Fig. 1 with k ad = 10 ( M · s ) − , k as = 10 M s − , and subunitbinding free energies g b are indicated on the plots. the subunit-subunit binding energy g b on p H, we only presentpredictions of the equilibrium and kinetic theories for packag-ing efficiencies at p H = 5 for several binding energies. Wewill compare these predictions to the currently existing data,but note that they should be compared to packaging efficien-cies from experiments carried out entirely at p H = 5 , as wasthe case for the light scattering data shown in Fig. 1b.As shown in Fig. 8a, the equilibrium theory predicts effec-tive packaging only above a threshold functionalized chargedensity, which depends on the binding free energy. Effec-tive packaging even occurs at a binding free energy g b = 0 ,showing that assembly on cores can occur under conditionsfor which spontaneous empty-capsid assembly is not favor-able, as seen in experiments [24] and simulations [43, 44]. On1core surfaces nanoparticle-subunit electrostatic interactionsand attractive subunit-subunit interactions overcome subunit-subunit electrostatic repulsions to drive assembly ( g H = − . kcal/mol per subunit for g b = 0 , Eq. 1).As shown in Fig. 8b, for weak subunit-subunit interactions( g b ≥ − kcal/mol) the kinetic theory predictions at sec-onds are essentially identical to the equilibrium results. Assubunit-subunit interactions increase in magnitude, however,the predicted threshold charge density becomes larger than theequilibrium value, and even begins to increase for g b < − .With stronger subunit-subunit binding interactions, the assem-bly of empty capsids competes with assembly on core sur-faces by depleting the concentration of free subunits. Assem-bly on cores remains favorable at high surface charge den-sities because initial subunit adsorption leads to rapid nucle-ation, and subsequent nonspecific adsorption enhances capsidgrowth rates (see below). The predicted packaging efficien-cies at long times are relatively insensitive to variations of thekinetic parameters k ad and k as .Interestingly, the threshold surface charge density predictedby the kinetic theory, σ C ≈ . , is nearly independent of thesubunit binding energy for g b ≤ kcal/mol, and comparableto the lowest experimental charge densities with measurablepackaging efficiencies (Fig. 1). However, the experimentalresults show a gradual increase in packaging efficiency as thefunctionalized charge density increases. We do not observesuch a gradual increase in packaging efficiencies except forshorter observation times or an extremely low surface assem-bly rate constant. This discrepancy could arise because wehave assumed monodisperse cores that are perfectly matchedto the capsid size, or it may indicate the presence of kineticallyfrustrated states (see below). B. Assembly kinetics on core surfaces
In this section we explore the effect cargo-capsid proteininteraction strength on assembly kinetics. As discussed inSection I, assembly kinetics can be monitored in experimentswith time-resolved dynamic light scattering; preliminary ex-perimental data for the highest functionalized charge densityis shown in Fig. 1. Because it is unlikely that the light scatter-ing signal can distinguish between disordered and assembledsubunits on core surfaces, we estimate light scattering fromthe kinetic theory by calculating the mass-averaged amountof adsorbed and assembled ( n + m ) subunits on core surfaces. L i gh t sc a tt e r [ A . U .] t [sec] σ c =0.8 σ c =1 σ c =2 σ c =3 σ c =5 (a) L i gh t sc a tt e r [ A . U .] t [sec] σ c =3 σ c =0.9 σ c =0.8 (b)
10 20 30 40 50 60 70 80 90 0.1 1 10 100 1000 L i gh t sc a tt e r [ A . U .] t [sec] k as =10 k as =10 k as =1 (c) FIG. 9: The predicted time dependence of light scattering for func-tionalization with strong acid groups. (a) and (b)
Indicated func-tionalization densities with rate constants (a) k as = 10 M s − and (b) k as = 10 M s − . (c) Indicated surface assembly rate con-stants, with σ C = 3 nm − . Other parameters for (a) - (c) are k ad =10 ( M s ) − , r CP:NP = 1 . , C p = 10 µ M, and g b = − kcal/mol. We simplify the presentation by varying only the charge den-sity ( σ C ) and the surface assembly rate constant ( k as ), withphysically reasonable values for the subunit binding free en-ergy, g b = − kcal/mol, and the the subunit adsorption rate2 L i gh t sc a tt e r [ A . U .] t [sec] (a) σ c =1.5 σ c =1.4 σ c =1.3 σ c =3 0 10 20 30 40 50 60 70 80 90 0.1 1 10 100 1000 L i gh t sc a tt e r [ A . U .] t [sec] σ c =1.5 σ c =1.4 σ c =1.3 σ c =3 (b) L i gh t sc a tt e r [ A . U .] t [sec] k as =10 k as =10 k as =10 (c) FIG. 10: The predicted time dependence of light scattering for sur-face charge groups with pK a = 4 . and pH = 5 , for (a) indicatedfunctionalization densities k as = 10 M s − , (b) indicated function-alization densities with k as = 10 M s − , and (c) indicated surfaceassembly rate constants with σ C = 3 nm − . Other parameters are asin Fig. 9. constant, k ad = 10 ( M · s ) − [82]. As discussed in sectionSection IV, the effect of reduction in dimensionality due toadsorption onto a surface is accounted for in the kinetic the-ory, but assembly rate constants could still be dramatically different from those in bulk because of impeded diffusion ofadsorbed subunits. We therefore explore a wide range of sur-face assembly rate constants k as .Predicted light scatter as a function of time is shown forvarious functionalized charge densities σ C and surface assem-bly rate constants k as in Figs. 9 and 10 for functionalizationwith strong and weak acids[83], respectively. The analysis iseasiest for the case presented in Fig. 9a ( k as = 10 M s − ), forwhich assembly on the core surface is slow compared to sub-unit adsorption. For each value of σ C , there is a rapid initialrise in light scatter, followed by a charge-dependent plateau,and finally a slow increase to saturation. The initial fast ki-netics correspond to nonspecific subunit adsorption (i.e. withlittle or no assembly) and the plateau occurs when adsorptionsaturates at what would be the equilibrium surface density ifthere were no assembly ( ρ eqsurf , see Section II). If adsorbed sub-units assemble into a relatively stable intermediate, the chem-ical potential of remaining adsorbed subunits decreases, re-sulting in further adsorption and increased light scatter. Thisprocess continues as additional subunits bind to the growingintermediate, until assembly stops.At high functionalization densities, nonspecific subunit ad-sorption approaches the density of a complete capsid; the highlocal or surface concentration of subunits enables rapid nu-cleation and assembly and thus the light scatter reaches itsfinal saturation value at relatively short times. At the lowerfunctionalization densities, nonspecific adsorption saturatesat smaller surface concentrations, with corresponding smallerinitial plateau heights for the light scatter. Nucleation timesincrease dramatically as the surface concentration is reduced(see the discussion in Ref. [44]), resulting in the long plateaubefore assembly completes.For a larger assembly rate constant k as = 10 M s − ,the nonspecific adsorption and assembly phases can only beclearly distinguished at low functionalized charge densities,as shown in Fig. 9b; at larger σ C adsorption and assemblyare concomitant. For weak acid groups (Fig. 10), equilibriumsurface densities increase only slowly with σ C ( ρ eqsurf ≤ . ,see Fig. 5) because the local concentration of negative chargedisfavors carboxylate dissociation. There is still a strong de-pendence of overall assembly rates on σ C , however, becausethe positive charge from capsid proteins drives additional car-boxylate dissociation (see Section II). Thus, assembly effec-tiveness cannot be predicted from ρ eqsurf alone, as was found forsimulations of neutral subunits in Ref. [43].3 Cores can enhance assembly rates.
From the precedingdiscussion, it is clear that functionalized charge leads to anincreased concentration of capsid proteins nanoparticle sur-faces relative to that in bulk solution, which can increasecapsid nucleation rates and even induce assembly below thethreshold protein concentration for formation of empty cap-sids. Because nonspecific adsorption onto core surfaces canbe much faster than subunit-subunit binding rates, the coresurface can also enhance assembly rates during the growthphase by increasing the effective subunit-subunit interactioncross-section. This effect was hypothesized for RNA by Hu etal. [42] and is reported for simulation results in Ref. [44].
Assembly rates decrease as capsids near completion.
Netassembly rates decrease as capsids near completion becauseadsorption of subunits onto cores with large partial capsidsis impeded by two effects. First, the excluded volume of thepartial capsid decreases adsorption rates; a similar effect hasbeen observed in simulations representing empty capsid as-sembly, but the magnitude of rate reduction seems to be modeldependent (see Refs. [44, 49]). Electrostatic repulsions withadsorbed and/or assembled subunits on the core surface canhave a larger effect on subunit adsorption. When the subunitsurface density is large compared to ρ eqsurf , the surface captureand diffusion assembly mechanism that enables rapid core-controlled assembly becomes inefficient. Even though directassociation of free subunits onto core-bound intermediates isincluded in the kinetic theory, there is a reduction in assemblyrates as capsids near completion; this effect is most noticeablefor surface functionalization with weak acid groups (Fig. 10a),or low fractions of strong acid groups ( σ C < nm − ). Thedecreased assembly rates are reflected in slow growth of lightscatter as assembly nears completion. Comparison with experimental data.
As discussed above,there is not sufficient experimental data to estimate values forthe unknown parameters; however, the predicted light scatterplots share some features with the experimental light scatterdata shown in Fig. 1b. In particular, for some parameter setslight scatter increases rapidly at short times, without the lagtime seen for empty capsid assembly [58, 59, 60]. It is tempt-ing to attribute the initial rise to nonspecific adsorption thatdoes not saturate until nearly complete coverage, as seen forstrong acid surface groups (Fig. 9a), but the calculation withweak acid groups predicts values of ρ eqsurf that are well belowcomplete coverage even for σ C = 3 nm − . Interestingly, thepredicted light scatter curves in this case with fast assembly FIG. 11: Estimated light scatter from simulations in which subunitsassemble around model nanoparticle with a size mismatch. Thepreferred empty-capsid morphology is a T=3 capsid, but the low-est free energy configuration is a T=1 capsid assembled around theT=1-size nanoparticle. Estimated light scatter is shown as a functionof time for different values of the parameter ǫ AA (in units of k B T ),which gives the free energy cost for adopting a subunit conformationconsistent with a T=1 morphology. As described in Ref. [43], thelarger values of ǫ AA lead to the formation of frustrated nanoparticle-associated partial capsids with a T=3 morphology, which cannotclose around the nanoparticle and block assembly into the lowest freeenergy configuration. The data is replotted from simulations reportedin Fig. 6 of Ref. [43], with the nanoparticle-subunit interaction en-ergy ǫ s = 8 k B T . kinetics ( k as = 10 M s − , Fig. 9a) still demonstrate a rapidinitial rise followed by slow saturation; the initial rise corre-sponds to rapid nucleation and assembly concomitant with ad-sorption. Thus, as discussed in the next section, comparisonof theoretical predictions and experimental light scatter dataat varying functionalized charge density will be necessary tounequivocally determine the origin of the initial light scatterphase.The experimental light scattering data (Fig. 1b) appears todemonstrate a logarithmic increase in light scattering between20 and 200 seconds. While the theoretical light scatter curvesalso demonstrate slow growth near completion (see above)due to subunit-subunit electrostatic repulsions, this effect islimited to the last five or 10 subunits. We did not find anyparameter sets for which predicted light scatter grew logarith-mically over several decades in time, although our parametersearch was not exhaustive. Intriguingly, the predicted lightscatter from simulations in which frustrated off-pathway in-termediates block assembly [43] do show logarithmic growthover several decades (see Fig. 11).4 C. Concentration effects
Varying nanoparticle concentrations.
We first examineassembly effectiveness as the nanoparticle concentration isvaried at fixed subunit concentrations. We define the cap-sid protein-nanoparticle stoichiometric ratio as r CP:NP = C p / ( N C C ) so that r CP:NP = 1 when there is exactly enoughcapsid protein to assemble a complete capsid on every core.The variation of packaging efficiencies with the stoichiomet-ric ratio is shown in Fig. 12; the predicted curves are sigmoidalin shape, much like the experimentally observed packaging ef-ficiencies shown in Fig. 1A of Ref. [24]. Packaging efficien-cies for r CP:NP . are typically lower than the equilibriumvalues obtained by solving Eqs. 7 and 8, because initially cap-sid proteins undergo nonspecific adsorption onto every core,which depletes the concentration of free subunits. The forma-tion of complete capsids then requires desorption and subse-quent re-adsorption onto cores with large intermediates. Thetimescale for subunit desorption increases with functionalizedsurface charge density and therefore the kinetic packaging ef-ficiency is nonmonotonic with respect to σ C for r CP:NP ≤ .This scenario is consistent with the slow assembly kinetics ex-perimentally observed for capsid assembly in the presence ofRNA [13]. Varying subunit concentrations.
As shown in Fig. 13, at σ C ≥ acid groups/nm the kinetic theory predicts a muchweaker dependence of assembly effectiveness on subunit con-centration than has been seen for empty-capsid assembly. Thistrend arises because the chemical potential of adsorbed sub-units is dominated by electrostatics rather than subunit trans-lational entropy. VI. DISCUSSION
The feasibility of estimating the unknown parameters.
The results in the last section are consistent with some featuresof the preliminary experimental data, and the predicted de-pendencies of packaging efficiency and light scattering on thefunctionalized charge density and nanoparticle-protein stoi-chiometric ratio can be qualitatively compared with additionalexperimental data. Quantitative comparisons of experimentaldata, however, will require estimates for at least three cur-rently unknown parameters, the subunit-subunit binding en-ergy g b , the assembly rate constant on core surfaces k as , thesubunit adsorption rate k ad . In Appendix A, we show that, P a ck ag i ng E ff i c i en cy σ C [acid groups / nm ] σ C =2.0 σ C =1.5 σ C =1.3 FIG. 12: Packaging efficiencies depend on the capsid protein/corestoichiometric ratio. The predictions of the kinetic theory for thevariation of packaging efficiencies with the capsid protein/core sto-ichiometric ratio, r CP:NP = C p /NC C at a fixed subunit concentra-tion of C p = 10 µ M, are shown for functionalization with weak acidgroups, with other parameters as in Fig. 8. -9 -7 -5 P a ck ag i ng E ff i c i en cy C p (M) σ C =2.0 σ C =1.6 σ C =1.2 FIG. 13: Packaging efficiencies depend only weakly on the initialsubunit concentration, C p , for fixed r CP:NP . The predictions of thekinetic theory for the variation of packaging efficiencies with sub-unit concentration C p are shown for functionalization with weak acidgroups, with other parameters as in Fig. 12. with good estimates for the unknown parameters, the theoret-ical predictions agree reasonably with simulation results for amodel of core-controlled assembly. While it is not certain thata simplified theory with a few parameters will exactly matchexperimental data from a system with thousands of distinct in-termediates, each of which could have different rate constants,rate equations with a similar level of coarse graining have suc-cessfully explained many features of empty-capsid assembly(e.g. [29, 30, 31, 32, 33, 34]). Given the number of experi-5mental control parameters in the nanoparticle-capsid assem-bly system, it should be possible to generate sufficient experi-mental data to further test the qualitative trends predicted here,and to make quantitative estimates for the unknown parame-ters through data fitting. With estimated values for k as and k ad ,it would be possible to predict the assembly state of subunitson core surfaces during the initial rise in light scatter, whichwould be challenging to determine with experiments alone.In addition to measuring packaging efficiencies and lightscattering data over a range of functionalized charge densities(with both weak and strong acid groups) and capsid protein-nanoparticle stoichiometric ratios, the theoretical predictionsin the last section suggest additional experiments that couldbe useful to provide independent estimates of some unknownparameters. In particular, the predicted light scatter curves inFigs. 9 and 10 show that adsorption and assembly happen si-multaneously for many parameter sets. Measuring light scat-ter of assembly incompetent proteins in the presence of func-tionalized nanoparticles would enable independent estimationof the subunit adsorption rate and the equilibrium subunit sur-face charge density ρ eqsurf ( Fig. 5). In addition to varying thefunctionalized charge density, the subunit-nanoparticle inter-action could modified by mutations to the N-terminal arm oncapsid proteins[84]; Tang et al. [85] have shown that cowpeachlorotic mottle virus capsid proteins with 34 residues deletedfrom the N-terminal arms are assembly competent (althoughthey lose selectivity for the native capsid structure). Identifying metastable disordered states.
We make sev-eral important approximations in this work. Most signifi-cantly, the theory considers only only intermediates consis-tent with the native capsid geometry, while experiments [70]and simulations [43, 55] indicate that other capsid morpholo-gies and/or asymmetric, malformed structures can be kineti-cally or even thermodynamically favored when core-subunitinteraction strengths are large compared to the thermal en-ergy k B T . These effects could be accounted for in the theoryby extending the state space to include structures other thancomplete and partial well-formed capsids, and by introduc-ing subunit diffusion constants that depend on the strength ofcore-subunit interactions. If the present theory proves capableof describing experiments with parameters for which predom-inately well-formed capsids assemble, however, inconsisten-cies between theoretical predictions and experimental mea-surements for other system parameter values could identifykinetic traps. For example, the calculated packaging efficien- cies show a less gradual increase with functionalized chargedensity than measured in experiments and predicted light scat-ter does not fully capture the a logarithmic increase in theexperimental light scatter curve. Both of these discrepanciescould be explained by frustrated off-pathway of assembly in-termediates or slow subunit diffusion rates, as evidenced bypredicted light scatter from simulations which do demonstratefrustrated assembly (Fig. 11). Confirming this possibility willrequire additional experimental data. VII. CONCLUSIONS
In summary, we develop simplified thermodynamic and ki-netic theories that describe the simultaneous assembly of vi-ral proteins and encapsidation of cargo. When applied to theencapsidation of charge-functionalized nanoparticles, the the-ories predicts a transition from inefficient core encapsidationto nearly 100% incorporation efficiency as the functionalizedcharge density is increased beyond a threshold value, and theestimated light scatter signal shows an initial rapid increasefollowed by slow rise to saturation, as seen in experiments.The predicted increase in incorporation efficiency with in-creasing charge density is sharper than seen in experiments,which could indicate the presence of kinetic traps that are notaccounted for in the present theory; comparison of theory pre-dictions with experimental data collected with a wide rangeof control parameters will be important to assess this possi-bility. The trends predicted here for varying surface charge,subunit concentration, and capsid-nanoparticle stoichiometricratio could guide the design of experiments that identify fun-damental principles and/or additional complexities for simul-taneous assembly and cargo encapsulation.
Acknowledgments
I am indebted to Bogdan Dragnea andthe Dragnea group for insightful conversations and for provid-ing me with experimental data, and to Oren Elrad for Fig 11.Funding was provided by Award Number R01AI080791 fromthe National Institute Of Allergy And Infectious Diseases andby the National Science Foundation through the Brandeis Ma-terials Research Science and Engineering Center (MRSEC).
APPENDIX A
In this appendix we compare predictions of the kinetic the-ory for the time dependence of the number of adsorbed parti-cles, n + m , to results presented in Ref. [44] from a computa-6 Parameter Value Definition N
90 Number of subunits in a complete T =3 capsid A C Inner surface area of a complete capsid (Section III C a g b [ − , kcal/mol Binding free energy per subunit-subunit contact (Eq. 9) g H [ − . , − . kcal/mol Free energy per subunit-subunit contact due to attractive interactions (Eq. 1) k ad ( Ms ) − Subunit adsorption rate constant ˆ k as M s − Empty capsid assembly rate constant [86] (Eq. 11) k as [10 , ] M s − Surface assembly rate constant (Eq. 15) σ C [ − , charge/(nm gold surface) Surface-density of functionalized charge q p
18 Number of positive charges per protein-dimer subunit C p µ M Total subunit concentration r CP:NP [0 . , Capsid protein–nanoparticle stoichiometric ratio, r CP:NP = C p / ( NC C ) C C [0 . , . µ M Nanoparticle concentration
TABLE I: Parameter values used for calculations in this work. The parameters used for capsid assembly are as follows. The number ofsubunit-subunit contacts are n c j = 1 for j ≤ n nuc , n c j = 2 for j ∈ ( n nuc , N − n nuc + 1] , n c N − = 3 for j ∈ ( N − n nuc , N ) , and n c N = 4 . Thischoice obviates the need for different nucleation and elongation rate constants [31]. The nucleus size is n nuc = 5 . The forward and reversereaction degeneracies are s = ¯ s = 2 , s j = ¯ s j = 3 for ≤ j < N (the average value calculated from simulations in Ref. [44]) and s N = 1 , ¯ s N = N . The binding degeneracy entropy (Eq. 9) is S degen = k B log( s j / ¯ s j ) . tional model that represents the assembly of T=1 capsid-likeobjects around a commensurate nanoparticle. This compari-son elucidates the effect of approximations used in the kinetictheory. In particular, in rate equation approaches that assumea single reaction pathway as described above, or multiple re-action pathways [50, 51, 87], and even binding between largeintermediates [50, 51], only intermediates that are consistentwith partially assembled capsids are considered. The simu-lations, on the other hand, explicitly track the spatial coordi-nates of subunits and thereby make no a priori assumptionsabout assembly pathways. Because we can either calculateparameters or estimate them from independent simulations,this comparison is a stringent test of the validity of these ap-proximations. It also tests the utility of using the kinetic the-ory to describe experimental adsorption results. The proce-dure by which we determine parameters for the kinetic theorythat correspond to particular sets of simulation parameters isdescribed below, as well as modifications to the theory in or-der to represent neutral subunits (or high salt concentrations).We do not consider here the simulations reported in Fig. 11and Ref. [43], for which the nanoparticle is not commensuratewith the preferred capsid morphology.Adsorption kinetics predicted by the kinetic theory and ob-served in simulations are shown in Fig. 14 for several subunitconcentrations and subunit-subunit binding energies for an ad-sorption energy of ε c = 7 . The agreement between the kinetictheory and simulation results is surprisingly good, consider-ing that no parameters were adjusted to fit this data and, as discussed below, relatively crude estimates are used for thesubunit-subunit binding rate constants f and binding entropy s rot .The computational model considered in Ref. [44] consistsof rigid subunits, for which excluded volume interactions aremodeled by spherically symmetric repulsive forces, and com-plementary subunit-subunit interactions that drive assemblyare modeled by directional attractions. The lowest energystates in the model correspond to ”capsids”, which consist ofmultiples of 60 monomers in a shell with icosahedral sym-metry. The parameters of the model are the energy associ-ated with the attractive potential, ε b , which is measured inunits of the thermal energy k B T , and the specificity of thedirectional attractions, which is controlled by the angular pa-rameters θ m and φ m . Capsid subunits also experience shortranged isotropic interactions with a rigid sphere placed at thecenter of the simulation box; these interactions are minimizedwhen capsid subunits are adjacent to the surface of the sphere.Subunit positions and orientations are propagated accordingto overdamped Brownian dynamics, with the unit of time t = a / D , where D is the subunit diffusion coefficient.Full details of the model are given in Ref. [44] Free energies for neutral subunits with no internal struc-ture or high salt concentrations.
The free energies for partialcapsids on core surfaces are modified from the expressionsgiven in Section III C for the case of neutral subunits or highsalt concentrations, in which case interactions are limited toexcluded volume, directional attractions that drive assembly,7 ad s o r bed s ubun i t s t / 10 t ad s o r bed s ubun i t s t / 10 t FIG. 14: The kinetic theory predictions (smooth lines) and simula-tion results (noisy lines) for the time dependence of adsorbed andassembled subunits, n + m , are shown for the neutral subunit model(Eq. A1). The simulations consider T=1 capsids and commensuratecores, so a complete capsid is comprised of N = 30 dimer subunits.Curves at increasing height correspond to reduced subunit densitiesof ρ p a = ε c = 7 and subunit-subunit binding energies of (a) ε b = 9 . and (b) ε b = 11 . and favorable core-subunit interactions. The free energy ofa core with n adsorbed but unassembled subunits and an ad-sorbed intermediate of size m is written as G surf n,m = ( n + m ) ε c − T S mix m,n − T S ad n, m (A1)where ε b is the core-subunit surface free energy strength.The second term, S mix accounts for the entropy for two-dimensional motions of adsorbed subunits on the core surface.A simple approach would assume Langmuir adsorption, butwe find that much more accurate results are obtained by in-tegrating an empirical formula for the chemical potential of afluid of hard disks [88] because the Langmuir model overesti- mates subunit mixing entropy at high surface coverages: S mix ( η ) /k B = Z η dη (cid:20) −
78 log(1 − η ) + 2 η − η + 9 η − η ) (cid:21) (A2)with the packing fraction η = ( n + m ) a / [16 π ( R C + a/ ] .The last term in Eq. A1, S ad , is the entropy penalty for sub-units to be localized a smaller distance from the core surfacethan the subunit diameter, which we obtained by numericallyintegrating the partition function for an adsorbed subunit whencomparing the kinetic theory to simulations.Calculating the time dependence of adsorbed and assem-bled subunits requires choosing values for several theoreticalparameters; here we describe the parameter values and howthey were chosen.1. Subunit binding entropy penalty. Eq. 9 includes a term s rot for the entropy loss (in addition to the subunit mix-ing entropy included in g b ) for a subunit to bind to a cap-sid that accounts for translational restrictions on scalessmaller than the subunit diameter a and rotational re-strictions [45]; the entropy penalty should increase onlyslightly when a subunit changes from one to more thanone bond. The subunit-subunit binding entropy penaltyis calculated as described in Ref. [45]. There is a typoin Eq. (15) of that reference; the correct formula is s rot ≈ −
32 log β∂ u att ( r ) ∂r (cid:12)(cid:12)(cid:12)(cid:12) r =0 −
12 log ( βε b ) π θ m φ m (A3)2. Subunit binding rate constant f . The assembly rateconstant was chosen as f = 0 . (in the dimension-less units of Ref. [44]) by comparing kinetic theorypredictions to simulation results for the assembly ofempty capsids. The resulting fit (not shown) shows lessagreement between theory and simulation results thanwe find for core simulations, perhaps because approxi-mations such as an average assembly pathway withoutmultimer binding are less significant when subunits areconfined to the surface of a core.3. The effective surface concentration of adsorbed sub-units is increased if subunits are confined within lessthan a molecular diameter a by the core-subunit inter-action, so we take ρ surf n,m = na ( N − m ) exp ( S ad /k B ) . Thisfactor ensures that entropy losses due to subunit local-ization are not counted twice.4. Adsorption rate, k ad . The rate of subunit adsorption tocore surfaces in our simulations can be calculated from8the Smoluchowski equation, with forces due to the ad-sorption potential explicitly included. Instead, we usethe diffusion limited rate for the zero force case, buttake k ad = 4 πDR cut , where R cut = 2 . a is the subunit-core center of mass distance at the point when the in-teraction force becomes nonzero. This choice was val-idated by comparing theory and simulation results forthe time dependence of the number of adsorbed sub-units in the case of no assembly, i.e. ε b = 0 . Sincethe excluded volume of adsorbed subunits is accountedfor in Eq. A2 we set the adsorption blocking factor φ n,m = 1 (see Eqs. 12). This approach slightly overestimates the adsorption rate at moderate coverage, butyields the proper equilibrium surface density if there isno assembly. Because of the definition of the unit time t , the subunit diffusion constant is D = 1 / .We find that the agreement in most cases is surprisinglygood, considering the degree of error in estimating s rot and f ,and the approximation that there is only one reaction pathway. APPENDIX B
In this appendix we summarize the application of themethod of Scheutjens and Fleer [66] to polyelectrolytes. Ourpresentation closely follows that given in Bohmer et al. [67],except there is a different geometry, the Flory-Huggins pa-rameter for interactions between polymer segments is set to χ = 0 , and there are additional terms in the free energy equa-tions.We consider a lattice around an impenetrable sphere withradius R C , contacting bulk solution in layer M + 1 , with M large enough that the presence of the surface is negli-gible in layer z = M . All quantities are assumed uni-form within a layer z . Fixed numbers of polyelectrolytes andTEG molecules are grafted in layers z p = R cap /l = 26 and z TEG = R C /l + 1 = 17 , respectively, with the lattice dimen-sion l = 0 . nm.In the calculation, a segment of species i is subjected to apotential field u i ( z ) , which is given by (relative to the bulksolution) u i ( z ) = u ′ ( z ) + v i e Ψ( z ) (B1)where the term u ′ ( z ) 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 z . Thesecond term on the right-hand side accounts for electrostaticinteractions, with the valency v i for a segment of species i , e the charge on an electron, and Ψ( z ) the electrostatic potentialin layer z . The statistical weight of a segment for species i inlayer z , relative to the bulk solution, is given by the segmentweighting factor G i ( z ) = exp( − u i ( z ) /k B T ) (B2)so that the volume fraction of a free monomer of type i isgiven by φ i ( z ) = φ b i G i ( z ) .The volume fractions for segments in chain molecules arecalculated with a segment distribution function G i ( z, s | ,which gives the statistical weight of a chain conformationstarting with a segment 1, located anywhere in the lattice, andending in layer z after s − steps. The segment distributionfunction is calculated from a recurrence relation: G i ( z, |
1) = G i ( z ) G i ( z, s |
1) = G i ( z )[ λ − G i ( z − , s − | λ ( G i ( z, s − |
1) + λ G i ( z + 1 , s − | (B3)with λ − , λ , and λ the fraction of contacts for a lattice site inlayer z with respective layers z − , z , and z +1 (see Ref. [89]).The statistical weight G i ( z, s | r ) of a chain conformation thatstarts at segment r and ends with segment s in layer z is cal-culated in the same way. The segment distribution functionsfor grafted molecules are modified to constrain segment tobegin in the grafted layer as described in Ref. [90]. To ac-count for impenetrable surfaces, u i ( z ) = ∞ for z < z TEG forall species, and u i ( z ) = ∞ for z > z p for polyelectrolyte andTEG molecules.The volume fraction of segment s in a chain of r segmentsis determined from the joint probability that a chain conforma-tion starts at segment 1 and ends at segment s in layer z , anda chain conformation starts at segment r and ends at segment s in layer z : φ i ( z, s ) = C i G i ( z, s | G i ( z, s | r ) /G i ( z ) (B4)where the probabilities are divided by G i ( z ) to avoid doublecounting of segment s and C i is a normalization constant. Ifmolecules of species i are free to exchange with the bulk so-lution, C i is determined by the fact that all segment weightingfactors G i = 1 in bulk so that C i = φ b i /r i . (B5)9If the total amount of molecules of species i is fixed (i.e. for agrafted brush) C i is determined from C i = θ i / [ r i G i ( r | (B6)with the total amount of molecules of species i calculatedfrom θ i = M X z =1 L ( z ) φ i ( z ) (B7)and the chain weighting factor G i ( r |
1) = P Mz =1 L ( z ) G i ( z, r | , which gives the statistical weightfor a chain of type i to be found anywhere in the lattice. Thenumber of lattice sites L ( z ) depends on the layer z becauseof curvature and is calculated in Ref. [89]. The electrostatic potential.
Electrostatics are accounted forwith a multi-Stern-layer model in which all charges within alattice layer are located on the plane at the center of the layer,with no charge outside of that plane. The plane charge densityin layer z , σ ( z ) is calculated from a sum over all species σ ( z ) = X i α i ( z ) v i eφ i ( z ) /a s (B8)where we set the cross-sectional area per lattice site a s = l with l the distance between layers. For weak acid groups thedegree of dissociation α i ( z ) is given by[68] α i ( z ) = K i h φ H O ( z ) ih φ H O + ( z ) i + K i h φ H O ( z ) i (B9)with the dimensionless dissociation parameter K i related tothe dissociation constant K D by K i = K D /c i with c i the mo-larity of pure segments for species i .The electrostatic potential with respect to bulk solution atlayer M is calculated through electroneutrality [67] and thepotential in each layer is calculated from Gauss’ law, forspherical layers it is [68] Ψ( z + 1) = Ψ( z ) − l πl (cid:20) ǫ ( z )( z − / z +1 ǫ ( z + 1) z ( z + 1 / (cid:21) z X z ′ =1 L ( z ′ ) σ ( z ′ ) (B10)with ǫ ( z ) the volume fraction averaged dielectric permittivityfor layer z . Free energies.