INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSIONS APPENDIX REFERENCES

Virus capsid assembly is an example of protein polymerization. Most spherical virus capsids exhibit icosahedral quaternary structure. They are constructed from multiples of 60 monomeric or oligomeric structural subunits, arranged in equivalent or quasi-equivalent environments (Caspar and Klug, 1962). To be viable, virus capsids must assemble efficiently and with high fidelity. The specific mechanisms of capsid assembly are poorly understood for spherical viruses. There is little experimental information describing intersubunit binding energies, rates, and orders for assembly reactions, and formation of nucleating structures.

In vitro capsid assembly has been observed for many viruses. Cowpea chlorotic mottle virus (CCMV) provided the first example of in vitro assembly of a spherical virus (Bancroft et al., 1967). Both empty CCMV capsids and infectious RNA-filled virions have been assembled (Bancroft and Hiebert, 1967; Bancroft et al., 1968; 1969; Adolph and Butler, 1976; Zhao et al., 1995; Fox et al., 1998). Characteristics common to capsid assembly were first observed with CCMV, including extreme sensitivity to solution conditions, concentration dependence, and an apparent critical concentration (Adolph and Butler, 1976). Since 1967 many other systems have been studied. Polyoma virus assembly was pleiomorphic, resulting in the production of icosahedral T = 1, octahedral, and T = 7 particles, as well as sheets and other polymers (Salunke et al., 1989). Poliovirus assembly from oligomeric intermediates demonstrated the fragility of the provirion form of picornaviruses (Rombaut et al., 1990; Basavappa et al., 1994; Rueckert, 1996). Assembly of related bacteriophages MS2 and R17 required a specific nucleic acid oligomer (Beckett et al., 1988; van den Worm et al., 1998). In vitro assembly of bacteriophage P22 demonstrated a relationship between scaffold and capsid proteins; these experiments also demonstrated sigmoidal capsid assembly kinetics for the first time (Prevelige et al., 1993). We present a framework linking experimental data to physically meaningful assembly parameters.

Previously, we have developed models based on thermodynamics and kinetics (Zlotnick, 1994; Zlotnick et al., 1999). Capsid accumulation exhibits the sigmoidal kinetics typical of multistep reactions (c.f. discussion of consecutive reactions in Fersht (1999)). During the initial lag phase, successive intermediates transiently accumulate and are consumed in turn by subsequent reactions. However, without regulation, the assembly reaction is susceptible to kinetic trapping (see Caspar (1980); Zlotnick (1994); Zlotnick et al. (1999, 2000)). Several regulating mechanisms have been proposed that decrease susceptibility to kinetic trapping. Autostery, in which free subunits equilibrate between assembly-competent and -incompetent states, has been proposed (Caspar, 1980); a form of it was used in the local rules model of quasi-equivalent assembly to decrease susceptibility to trapping (Berger et al., 1994; Schwartz et al., 1998). We have focused on regulation through nucleation by incorporating a "kinetically limiting" early step in the reaction, described as KL assembly (Zlotnick et al., 1999), wherein the first few intermediates either assemble slowly or are unstable. Nucleation is common with other biopolymers and has been observed in virus capsid assembly (Sorger et al., 1986; Beckett et al., 1988; Bartenschlager and Schaller, 1992), where a nucleating structure might be a completed vertex or a closed ring of subunits, much as the completion of the first turn nucleates a helical filament (Oosawa and Asakura, 1975). These and other regulatory mechanisms can be incorporated into the general models discussed here.

In this paper we further develop techniques for analyzing assembly (Zlotnick, 1994; Zlotnick et al., 1999) using two specific model systems to illustrate the kinetic analysis derived more generally in the Appendix. The example systems show two specific cases of the differential rate equations (A2) for subunit, intermediate species, and capsid concentrations. In its complete generality the model system allows the assignment of different rates and/or orders for each of the cascade of assembly reactions and is in some respects the spherical analog of a helical filament model (Flyvbjerg et al., 1996). Fundamental geometric differences between filament and shell assemblies require distinct mathematical treatments, as the former is open-ended and the latter is closed.

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSIONS APPENDIX REFERENCES

See Table 1 for naming conventions.

Modeling capsid assembly at equilibrium

Two model capsids were used for calculating simulations (Fig. 1). The first, initially developed without nucleation as the equilibrium (EQ) model in Zlotnick (1994), consists of 12 pentavalent pentagonal subunits assembling as faces of a dodecahedron by making 30 intersubunit contacts (Fig. 1 A). This geometry resembles picornavirus quaternary structure (Rueckert, 1996). The second system is built from 30 tetravalent rectangular subunits (making 60 intersubunit contacts), which can be thought of as a T = 1 icosahedron assembled from "dimeric" subunits (Fig. 1 B). The tetravalent subunit resembles the geometry used in schematic representations of hepatitis B virus capsid protein dimer (Zlotnick et al., 1996).

View larger version (31K): [in this window] [in a new window]   FIGURE 1   Capsid geometries. (A) The dodecahedral model capsid has 12 pentagonal subunits. (B) The icosahedral model capsid has 30 twofold symmetrical subunits. The subunits for the icosahedral model are tetravalent, with binding sites at each corner; subunits for the dodecahedral model are pentavalent, with a binding site at each edge.

To describe capsid polymerization we develop the system of rate equations for the concentrations of species formed in the assembly cascade detailed in the Appendix (Eq. A2). The equations are determined by capsid geometry, assembly path, and intersubunit binding energies. At equilibrium, the differential rate equations reduce to a series of algebraic equations from which it is possible to determine equilibrium values for subunit, capsid, and intermediate species. For a simplest model, at equilibrium, only monomer and capsid are present at appreciable concentrations (Zlotnick, 1994). The equilibrium relationship between the reaction cascade endpoints is:

(1)

where N is the number of subunits in a complete capsid, e.g., N = 12 for the dodecahedral and N = 30 for the icosahedral model.

Equation 1 results in the unwieldy units of M^1N for the association constant K_Acap. Two expressions that are more convenient are K_Dapp, the apparent dissociation constant defined where equilibrium concentrations of subunit and capsid are equal (Eq. 2) and K_Acon, the pairwise (per-contact) association constant between subunits (see Eq. 3) (Zlotnick, 1994). Because N is typically large, it is more practical to use the logarithmic form of equations involving K_Acap when manipulating numerical data.

(2)

Successive substitutions through the sequence of model equations at equilibrium give K_Acap as the product of 1) a statistical coefficient accounting for the overall degeneracy of the reaction, and 2) a power of K_Acon accounting for each contact formed in capsid assembly. This can be written generally for the assembly of a polyhedral structure of N subunits with c contacts per subunit where each subunit has j-fold degeneracy. Such an N-hedral structure will have a total of C = cN/2 pairwise contacts, leading to the general form where j^N1/N is the product of the statistical factors for the degeneracy of the component reactions (Eq. 3).

(3)

Individual statistical factors can be determined explicitly for each assembly step as previously described (Zlotnick, 1994). Applying Eq. 3 to our dodecahedral and icosahedral models gives Eqs. 4.

(4)

In general, equilibrium is described for a system of rate equations by setting the derivative equal to zero. For models like those discussed here (Eqn A2), having a single representative conformation for each size of intermediate, this algebraic system is easily solved by back-substitution, giving the equilibrium concentrations of all species in terms of subunit

(5)

where c_m gives the total number of intersubunit contacts formed in assembling the mth species (composed of m subunits). For models incorporating several paths, the expression (Eq. 5) exhibits a complicated dependence on K_Acon and degeneracy.

In principle, K_Acon need not be the same for all interactions. For example, the K_Acon used for the nucleation interactions might be less than that used for contacts established during elongation. If this were a simple entropy penalty, one would expect the penalty to be recouped once the nucleus is assembled (Oosawa and Asakura, 1975). For simplicity, we use the same K_Acon for all reactions.

Modeling assembly kinetics

The above model of assembly allows calculation of assembly kinetics with few additional assumptions. We assume for simplicity that every contact scores the same per-contact association energy, G_c, from which we define K_Acon via the Arrhenius relation G_c = RT ln K_Acon. Second, we limit our simulations to only two microscopic forward rate constants, for nucleation (f_nuc) and elongation (f_elong), which are modified by appropriate statistical factors. Third, we have restricted assembly to proceed only by addition of individual subunits, without high-order or pseudo-high-order elongation steps. Finally, some consideration must be given to the choice of assembly paths that are used in kinetic simulations: a typical capsid may have hundreds or thousands of possible paths through tens or hundreds of intermediates. This is analogous to an energy landscape model of protein folding (Dinner et al., 2000; Dill, 1999). However, the first and last few reactions in capsid assembly have little degeneracy; the greatest degeneracy is encountered when assembly is half complete.

There are two obvious strategies for limiting the number of paths used in the model: 1) use the "best" paths, and 2) use the "average" path. In the first strategy, we define a measure of "significance" or the likelihood of a path. In the second, we specify only one intermediate of each size, the "average" intermediate of that size. Both of these paths follow steep energy descents with no intrinsic energy minima other than capsid; these paths are akin to a smooth, steep trajectory on an energy landscape.

One realization of the best-path strategy is to use paths through maximally stable intermediates, an msi-measure. For example, there are typically two configurations for the trimer: open (linear) and closed (triangular). Given that subunit geometry admits both, closed trimers are a lower-energy intermediate than open trimers because they gain one additional unit of contact energy from closure. The msi-measure assigns a measure of one to paths rooted in the most stable intermediates and zero to all others, deleting them. We have used the msi-measure to choose intermediates for our dodecahedral model. Statistical coefficients are still required to reflect degeneracy of the msi-path(s).

In the "average-path" strategy, energy and degeneracy are distributed evenly among the steps along a single assembly path, which may or may not be an actual path. Although we developed a version of our icosahedral model using the msi-measure, the behavior of the msi model was not significantly different from that of the average-path version with regard to the concentrations of capsid ([N]) and subunit ([u]). In the average-path construction for the icosahedral model there are 29 reactions over which we distribute the energy from the 60 contacts. The dimerization step yielded only one energy unit because there is only one contact. The next 26 steps were assigned two energy units each (two intersubunit contacts per subunit association) because actual paths have from one to three new contacts formed per step. The final two subunits make three and four contacts on association and were assigned three and four energy units, respectively. Reaction degeneracy was distributed in a similar manner.

Assembly kinetics are simulated as species concentrations [u], [2], ... , [N] evolve according to the differential rate equations (Eq. A2) from which we extract a typical example for convenient reference as Eq. 6.

(6)

This equation takes account of both subunit-intermediate association (forward "f-terms") and dissociation (backward "b-terms") contributing to the kinetics of the mth intermediate. The statistical factors s_m and _m adjust the rates to account for the degeneracy of the reactions; examples of their calculation for the dodecahedral model are given in Zlotnick (1994). The reverse rate constant for the mth reaction

(7)

is calculated from the forward rate (f_m), the per-contact energy (from K_Acon), the number (r) of contacts broken in dissociating the mth subunit, and the reverse reaction's degeneracy _m.

The system of rate equations (A2) developed from the above formalism is not amenable to exact solution, unlike filament assembly (Flyvbjerg et al., 1996). Solutions were calculated using BERKELEY MADONNA (Berkeley Software, Berkeley, CA), with either RK4 or Rosenbrock numerical integration methods, and Waterloo MAPLE (Waterloo, ON, Canada) implementations of the RKF45 integrator. MAPLE simulations were calculated using 14 of the 15 digits available in standard floating-point precision, consuming no more than a few minutes of CPU-time per initial concentration. In no case was high precision required to obtain accurate values for parameter estimators. Different models were fit to synthetic data sets by minimizing the RMSD, as implemented in BERKELEY MADONNA.

	RESULTS AND DISCUSSION

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSIONS APPENDIX REFERENCES

An overview of simulated assembly kinetics

Sigmoidal kinetics such as those exhibited by capsid assembly (Fig. 2) are typical for the synthesis of product by any multistep reaction (Fersht, 1999). For capsid assembly, the three stages of kinetics (lag, rapid assembly, and asymptotic) can be associated with well-defined characteristics of the reaction cascade.

View larger version (40K): [in this window] [in a new window]   FIGURE 2   Assembly kinetics calculated for the icosahedral model. The cascade of reactions exhibits the sigmoidal kinetics experimentally observed for in vitro capsid assembly. Subunit concentration declines. Intermediates successively rise during lag phase and then slowly decrease to extremely low equilibrium concentrations. After all intermediates are synthesized, rapid production of capsid begins.

The initial lag phase is characterized by the sequential building-up of assembly intermediates. Each intermediate reaches a peak concentration in its turn, then slowly relaxes toward its eventual equilibrium concentration. Until the final intermediate begins to approach its peak concentration, capsid production is very slow. The lag phase can be characterized as the initial interval during which a near steady-state "assembly line" of intermediates is constructed. In reactions with nucleation and elongation components, the elongation rate has its clearest effect on the lag duration (see assertion 4).

During the rapid stage of kinetics there is a burst of capsid synthesis as free subunits react with intermediates, effecting a nearly steady conversion of mass in free subunit into mass in capsid state. The rapid rate of capsid synthesis is supported by the relatively high concentrations of reactants, both free subunit and intermediate. During this stage, intermediate concentrations remain almost constant. Both nucleation and elongation rates make contributions to the slope of this part of the kinetic trajectory.

It turns out that assembly path degeneracy is not an important issue for simulations where intermediate concentrations are low. Different assembly paths lead to small changes in the concentrations of different intermediates, but not subunit or capsid. We found that the back reactions are insignificant in the elongation steps until reactions approach equilibrium. Because all assembly reactions proceed one subunit at a time at nearly identical forward rates, all possible paths will take roughly the same time to complete, although some paths are more likely than others.

The asymptotic phase occurs when the depletion of free subunit becomes significant. Low subunit concentration causes the slow rate of nucleus formation to be the major regulating factor in capsid synthesis (see assertion 2).

Conceptually, nucleation occurs when the first few intermediates in a series of reactions are less stable than subsequent products, as might occur in a helical polymer prior to the first interactions between turns (Oosawa and Asakura, 1975; Flyvbjerg et al., 1996). Nucleation of a spherical polymer might depend on formation of a primitive closed structure, such as the first face or vertex of a polyhedron (e.g., an HBV trimer (Zlotnick et al., 1999)). Energetically, nucleation reactions can be differentiated from subsequent elongation reactions by a slower forward rate and/or weaker association energy. For simplicity, we have chosen to focus on a rate-delimited nucleation step, observing that our results do not significantly depend on the mechanism by which nucleation differs from elongation.

Reaction simulations are less susceptible to kinetic traps when assembly is regulated by a nucleation step (Zlotnick et al., 1999). Nucleated reactions can yield near-equilibrium concentrations of capsid without trapping despite relatively high association energies, high initial concentrations of subunit, and fast rates (see assertion 1). Resistance to kinetic trapping results from prevention of over-initiation of assembly. Because the forward rate is proportional to subunit concentration, when subunit concentration is too high, virtually all subunits are assembled into intermediates in the absence of a regulatory step. In nucleated reactions, or in reactions "autosterically" regulated by the persistence of assembly-competent and -incompetent forms of the subunit (Caspar, 1980; Schwartz et al., 1998), unbound subunit is held in reserve. The introduction of an assembly-incompetent subspecies amounts to no more than the formal addition of a slow, initial, first-order reaction step. Mathematically speaking, it is not surprising that either mechanism resists trapping because both introduce slow coordinates in the early equations of the assembly cascade.

Nucleation and the steady-state assembly line

In the absence of a regulating step, successful assembly requires a balance between protein concentration and association energy (Zlotnick, 1994). As we increase either contact affinity, K_Acon, or initial subunit concentration, u₀, the amount of subunit in the assembly line increases. During the rapid-production phase, subunit can be thought of as being converted, via the assembly line, "directly" from free to bound-in-capsid state: the concentrations of intermediates approximate a steady state. When K_Acon or u₀ is too large, assembly is over-initiated; too much mass has entered the assembly line, leaving the reactions starved for available subunits: a kinetic trap. Conversely, too little mass in the assembly line will result in a slow reaction that is starved for intermediates. Because capsid production is slow for either extreme, there is an optimal steady-state concentration of intermediates for K_Acon at some value of u₀ (cf. assertion 1 and Fig. 4).

In nucleated assembly the concentration of free subunit falls much more gradually than during assembly lacking a regulating step (Zlotnick et al., 1999). Nucleated assembly is resistant to kinetic trapping, but not immune to it. The degree of resistance depends on the distinction between nucleation and elongation phases. For example, when the ratio of the forward rate constants, f_nuc/f_elong, is <~0.001, assembly may be extremely resistant to trapping; when f_nuc/f_elong > 0.05, the reaction resembles EQ assembly and may be easily trapped. A weaker contact energy for nucleation reactions enhances the resistance of nucleated assembly to kinetic trapping.

Apparent K_D and pseudo-critical concentrations

As u₀ is increased, the equilibrium capsid concentration increases rapidly because of the high algebraic degree of Eq. 1. As Fig. 3 shows, the point of greatest sensitivity to change in u₀ is well approximated by K_Dapp (Eq. 2), the subunit concentration for which equilibrium subunit and capsid concentrations will be equal. Because K_Dapp is calculated on a molar basis, the initial concentration required to achieve it is (N + 1)K_Dapp. In these models K_Dapp is a pseudo-critical concentration (Fig. 3). From Eq. 1 it is clear that free subunit concentration [u] will vary, albeit slowly, for u₀ (N + 1)K_Dapp. Because [u] increases so slowly as a function of u₀, there is the appearance of a constant [u], approximated by K_Dapp. This same effect was predicted and observed with micelle formation (Tanford, 1980). It is not surprising that critical concentrations have been indicated for CCMV (Adolph and Butler, 1976) and HBV (Seifer et al., 1993), whose capsids have 90 and 120 subunits, respectively, as CCMV and HBV would be expected to have much steeper isotherms than the 30-subunit model (Fig. 3).

View larger version (17K): [in this window] [in a new window]   FIGURE 3   An assembly isotherm for the icosahedral (30-mer) model. The mass fraction capsid (solid line, right axis) as a function of the initial subunit concentration, [u0]. Capsid concentration increases as a 30th power of free subunit (Eq. 3). Free subunit concentration (dashed line) varies little when it approaches KDapp, but it is not constant; i.e., there is no true critical concentration. Subunit concentration [u] will reach KDapp, the horizontal line at 3.78 µM [u], when [u] = [N]; on a mass fraction basis, this will occur at 117 µM total subunit with ~97% of the mass as capsid. The isotherm was calculated for log(KAcon) = 2.5, corresponding to Gc = 3.4 kcal/mol.

Scaling assembly simulations

We searched for intrinsic units of concentration and time for kinetic simulations (see Eq. A5) similar to those observed in the case of filament assembly (Flyvbjerg et al., 1996). Approximately the same trajectories are observed for reactions where concentration is expressed in units of u₀K_Acon and time in units of t/K_Acon. This scale invariance is intrinsic to the model equations A2, simplifying comparison of different kinetic simulations. The rescaling is most nearly exact during the lag and rapid-assembly stages when the forward-reaction terms dominate the kinetics and all cascade reactions are essentially unidirectional.

Assertion 1: K_Acon can be estimated from kinetic trajectories at long times

A reaction must be at equilibrium for precise determination of association energy, yet kinetic simulations of virus assembly reactions (Fig. 2) approach equilibrium asymptotically. We examined simulations over long times to observe how closely they approached equilibrium and to test our estimation of the pairwise association energy, G_c. Reactions were simulated for a range of K_Acon values. For any fixed concentration, reactions at modest K_Acon approached their equilibria rapidly; at higher values reactions were susceptible to stalling. Stalling occurs when the concentration of reactants is so low that bimolecular reactions do not proceed at an appreciable rate.

As assembly kinetics approach equilibrium, capsid concentration changes very slowly. Approximate values for G_c were calculated using the concentrations of capsid [N] and subunit [u] measured at long times in place of their equilibrium values in Eqs. 1, 3, and 4. Intermediates were not included in these calculations to make our results more "realistic" because intermediate concentrations are expected to be so low that they are experimentally undetectable. In long-time simulations with modest K_Acon, the estimated per-contact energy G_c was well within 10% of the exact value. Estimates for reactions with stronger association energies yielded worse approximations because these reactions stall when they became starved for free subunits and intermediates at long times. This is because a large K_Acon corresponds to a very low equilibrium concentration of subunit, [u] (Eqs. 1 and 3) so reactions reached a starvation condition and stalled before approaching equilibrium. Estimation of G_c was slightly more accurate at low initial subunit, even though higher initial subunit concentrations yielded a larger mole fraction of capsid. This seems counterintuitive in light of stalling due to starvation; however, this is simply an effect of u₀ being closer to [u], when u₀ is small. The dodecahedral model proved to be more stringent than the icosahedral model for this estimation (Fig. 4) because subunit geometry results in a lower K_Dapp for a given G_c, as can be seen by combining Eqs. 2 and 4 to generate the approximations of Eq. 8.

(8)

In either case, a K_Dapp much lower than ~0.2 µM will result in a stalled reaction.

View larger version (14K): [in this window] [in a new window]   FIGURE 4   Energy per contact, Gc, calculated from the KAcap observed at the end of kinetic simulations of a 24-h assembly experiment for the dodecahedral (open symbols) and icosahedral (closed symbols) models. The simulations all use fnuc = 100 M1 s1, felong = 106 M1 s1. Association constants were calculated using Eqs. 1 and 2, and converted into energy units using Gc = RTlnKAcon, where T is set to 298 K and R is the gas constant, 1.987 cal (deg mol)1. Simulations for the 12-mer model correspond to 2.5 (circles), 2.7 (squares), and 3.4 kcal/(mol contact) (triangles). Simulations for the 30-mer were calculated for Gc of 3.4 (circles), 4.1 (squares), and 4.8 kcal/(mol contact) (triangles). Though the ranges of Gc are different, the Gapp values are similar due to capsid and subunit geometry. Ignoring the relatively small statistical component, for the 12-mer, Gapp is ~30/11(Gc); for the 30-mer, Gapp is ~60/29(Gc) (see Eq. 7).

Assertion 2: Nucleation rate can be determined from the rate of capsid formation

The regulatory step has a profound effect on assembly kinetics. For nucleated assembly, both the size of the nucleus and the rate of nucleation must be considered. We consider the case where nucleation proceeds from slow association of subunits to form a metastable nucleus.

When a reaction enters the asymptotic phase, the assembly line of intermediates is still in steady state, though the rate of capsid production is attenuated. Because every new nucleus enters the assembly line and proceeds to form a capsid, we assert that the rate of capsid production is proportional to the rate of nucleus formation. More accurately, the average forward rate constant for the nucleation steps is proportional to the rate of capsid formation divided by [u]^nuc. This is expressed in

(9)

where K_subnuc is the equilibrium constant for the formation of [nuc  1]. The value for K_subnuc is the coefficient in Eq. 5. The derivation of Eq. 9 is provided in the Appendix (see A7-A12). We are able to arrive at Eq. 9 because in nucleated assembly 1) elongation reactions are essentially unidirectional, and 2) from relatively early times forward, subnuclear intermediates are already near their equilibrium concentrations. Equation 9 holds at early times in the asymptotic phase, before intermediates in the assembly line have begun to disassemble. Under these conditions, almost every nucleus produced leads to the formation of a complete capsid, i.e., the rate of capsid production is approximately the same as the forward-rate of nucleus production.

Given the size of the nucleus (see assertion 3), nucleus geometry and statistical factors can be deduced. The value for K_subnuc can be calculated from K_Acon (see assertion 1) and the statistical factors. Equation 9 is then an equation in only one unknown, f_nuc. Fig. 5 shows plots of the ratio

(10)

calculated for simulations using three values of f_nuc reflecting the impact that decreasing the forward nucleation rate has on the ratio in Eq. 10. The value of the ratio is concentration-independent after the reaction has reached steady state. The plateau height in each case is equal to the forward rate constant, f_nuc, multiplied by K_subnuc. We emphasize that the estimator depends only on the nucleation rate and not on the elongation rate. In simulations where the nucleation rate was fixed and the elongation rate varied, the plateau height remained fixed.

View larger version (37K): [in this window] [in a new window]   FIGURE 5   The rate of capsid formation can be related the nucleation rate, fnuc. Time and concentration series for the ratio (Eq. 10) of capsid production rate, d[N]/dt, to free subunit concentration raised to the nucleus size power, [u]nuc. The concentration-independent plateau height is directly proportional to the nucleation rate. The value of fnuc was set to 104 M1 s1 (A), 103 M1 s1 (B), and 102 M1 s1 (C). For these three families of reactions, felong = 106 M1 s1, and Gc = 4.1 kcal/mol.

In general, concentration independence of a quantity calculated from concentration-dependent data indicates that the calculated quantity is a parameter of the system. In this case, the concentration independence in the plots illustrates the dependence of the ratio (Eq. 10) on the coefficients in the governing equations (Eq. A2) and not on the initial conditions.

Assertion 3: Nucleus size can be determined from the concentration dependence of the extent of assembly

Simulations show that nucleus size cannot be determined accurately by searching for a best fit to the exponent in Eq. 9. In the analysis given in the Appendix (Eqs. A13-A18) we develop a sensitive estimator (Eq. 12) for nucleus size. This analysis rests on two assumptions: 1) the ratio in Eq. 10 is independent of u₀, and 2) the reactions considered are in a transient steady state where formation of capsid is equal to the depletion of free subunit (Eq. 11).

(11)

Our results (see Appendix) show that the nucleus size appears as the temporal minimum value of the function

(12)

The function L is the slope of a plot of ln[N] versus ln[u] at a single time across a range of initial protein concentrations, u₀ (Zlotnick et al., 1999). The minimum will only be found for a family of reactions that are simultaneously in steady state. This means that there is a range of times and concentrations for which the function L, near its minimum, approximates nucleus size.

Fig. 6 A shows the graph of the function L over a range of times and concentrations for our icosahedral model with nucleus size nuc = 3. The minimum value of L varies little from the nucleus size. Fig. 6 B shows the series of simulations used to calculate the surface in Fig. 6 A, with the shaded region delineating the domain in which nucleus size is approximated to within 0.33. Nucleus size is unambiguously approximated over a fairly wide range of initial concentrations at early time. For lower concentrations, the approximation holds at later times. Fig. 7 shows the surfaces z = L(u₀, t) for nucleus sizes 3, 4, and 5 in the dodecahedral model.

View larger version (40K): [in this window] [in a new window]   FIGURE 6   Determination of nucleus size for the icosahedral model. The minimum of the function L (Eq. 16) gives the nucleus size estimate. (A) The surface of L (light gray) is shown for the nuc = 3,  = 7 × 105, icosahedral model. The minimum of the surface is partially obscured by the error plane z = 3.33 and lies above the error plane z = 2.67, which are dark and medium gray, respectively. (B) The corresponding kinetic simulations (black) show where (shaded region) nucleus size is estimated to within ±0.3.

View larger version (48K): [in this window] [in a new window]   FIGURE 7   Determination of nucleus size for three different families of dodecahedral simulations. Nucleus size estimation surfaces z = L(u0, t) for the dodecahedral model with nucleus sizes of 3 (A), 4 (B), and 5 (C). The figure follows the use of grayscales as in Fig. 6.

Assertion 4: Elongation kinetics are most apparent in the early phases of assembly reactions

Examination of kinetic simulations (Fig. 2) shows that the lag phase is due to the time required to assemble intermediates. The duration of this phase is affected by nucleus size, nucleation rate, the stability of subnuclear intermediates (i.e., K_Acon), and elongation rate. Other parameters remaining constant, simulations indicate a localized dependence of the lag duration on the elongation rate, f_elong (Fig. 8). This is clearest when comparing a series of kinetic trajectories where only f_elong was varied (Fig. 8 A). Note that the different simulations coalesce at longer times, where the nucleation rate regulates the rate of capsid formation (assertion 2). When examining a range of subunit concentrations, the effect of f_elong on the lag duration is most evident at lower [u₀] (Fig. 8, B-E). The variation in lag duration makes it clear that it should be possible to approximate f_elong by curve-fitting.

View larger version (49K): [in this window] [in a new window]   FIGURE 8   The effect of elongation rate is most obvious in the early stages of assembly kinetics. (A) Elongation rate is varied for a constant initial subunit concentration. The elongation rates are 105, 4.6 × 104, 2.2 × 104, and 104 M1 s1, from left to right; the nucleation rate is 100 M1 s1. At long times, where the nucleation governs the rate of capsid formation, the four curves are coincident, as predicted for a common nucleation rate. (B-E) Families of kinetic simulations where [u0] is 20, 10, 6, and 2 µM, corresponding to [u0]/KD of 0.063, 0.032, 0.019, and 0.0063. For these reactions, felong is 105 (B), 4.6 × 104 (C), 2.2 × 104 (D), and 104 M1 s1 (E). All simulations in this figure were calculated using a 30-mer model with a trimeric nucleus, Gc = 4.1, and fnuc = 100 M1 s1.

For small nuclei in large capsids f_elong becomes more important because of the many post-nuclear steps. The elongation rate principally regulates the lag phase while the nucleation rate regulates all phases. Because the nucleation rate affects all phases (Eq. A2), it approximates a global time rescaling factor.

Because the rapid phase includes contributions from f_nuc and f_elong, we examined the concentration dependence of ln(d[N]/dt) against ln u₀, a classical kinetic approach for estimating reaction order in single-step reactions. We observe the concentration dependence of maximal rate is ~1. This is consistent (though not conclusive) with a pseudo-first-order reaction, as previously suggested (Zlotnick et al., 1999). There is no simple relationship between the dependence of ln(d[N]/dt) on ln u₀ and nucleus size.

Sensitivity of kinetic trajectories to model and parameters

How reliable are parameters estimated from assembly kinetics? How well can we discriminate between the correct model of an assembly reaction and a model that mimics the data but does not represent a physically reasonable description of the reaction?

We find that the models are sensitive to variations in parameters. A simulation produced using the wrong parameter values, though it may mimic data at one concentration, could not fit a concentration series. We demonstrate this by examining two reactions (Fig. 9). In the first case we fit K_Acon, f_nuc, and f_elong for an icosahedral (30-mer) model to data from a dodecahedral simulation. In the second case, the correct model was used to fit simulated data for an icosahedron, but we held K_Acon to an incorrect value. In both experiments residuals showed a systematic error that was more pronounced when calculated for a concentration series, though the initial fit may have looked attractive.

View larger version (26K): [in this window] [in a new window]   FIGURE 9   Fitting models to synthetic data. Data calculated for a dodecahedron (solid line, log KAcon = 3.3; fnuc = 30 s1; felong = 1500 s1) with a 10-µM subunit were fit with an icosahedral model (dashed line, log KAcon = 3.84; fnuc = 80 s1; felong = 21,760 s1) (A, B); data calculated for an icosahedron with a 10-µM subunit (solid line, log KAcon = 3.5; fnuc = 100 s1; felong = 30,000 s1) were fit with an icosahedral model with the log(KAcon) restrained to 4 (dashed line, log KAcon = 4.0; fnuc = 44 s1; felong = 440,000 s1) (C, D). Kinetic trajectories of the starting model and fit model were also compared for 5 µM and 2.5 µM subunits. Residuals are expressed as the difference between the correct mass fraction capsid and the mass fraction calculated for the wrong model.

Calculated kinetic trajectories are sensitive to model geometry and kinetic-thermodynamic parameters. Fitting experimental data requires an appropriate model, e.g., a 120-mer msi-path model for hepatitis B assembly. A close fit across a concentration series indicates that the assembly is well represented by the model and that the parameters extracted from the experimental data (as outlined in our assertions 1-3) give a fair representation of the true parameters. We stress, however, that there is no need for such a comparison to extract parameter values from experimental data. A poor fit calculated with parameters extracted from the experimental data indicates that mechanisms not included in the initial model need to be included in further refinements. In this way the extraction of parameters followed by model fitting serves as a basic probe of the assembly process.

	CONCLUSIONS

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSIONS APPENDIX REFERENCES

In conjunction with studies of the assembly of CCMV and HBV, we have developed a simple model of capsid assembly that makes limited assumptions. Based on this mathematically expressed model, we have derived techniques to determine reaction parameters from kinetic data.

The basic assembly model consists of defined capsid geometry and four parameters: the microscopic, per-contact equilibrium constant (K_Acon), the nucleus size (nuc), the nucleation on-rate (f_nuc), and the elongation on-rate (f_elong). We have shown how good approximations of three of these parameters can be extracted from the kinetics of capsid formation. In our experience, these are the most accessible data from assembly experiments. We have illustrated these methods by extracting parameters from simulated data calculated from 12- and 30-subunit models. Practically, the larger model gave more robust results, as mathematically expected. The mathematical analysis was carried out for an arbitrary capsid size and spherical capsid geometry consistent with equivalent subunits. Estimation of other assembly parameters will most likely require fitting observed kinetics with a more detailed model of the assembly reaction.

In the first few minutes of a typical in vitro assembly reaction, ~10¹² nuclei form and go on to yield capsids (Zlotnick et al., 1999, 2000); this is unlike polymerization of a crystal, where a single nucleus leads to formation of a polymer of up to 10¹² subunits. The general model of capsid assembly described in this paper is a logical outgrowth of our previous work (Zlotnick, 1994; Zlotnick et al., 1999, 2000). To our knowledge, it is the only rigorous description of the thermodynamics and kinetics specific to spherical polymerization. Using quasi-molecular dynamics simulations to examine assembly of individual capsids, Berger and co-workers (Schwartz et al., 1998) were able to reproduce the sigmoidal assembly curves expected for assembly of large populations, essentially a statistical mechanical approach to capsid assembly, though their system was not suitable for developing analyses of experimental data. However, our results are distinctly different from analyses of capsid assembly based on formation of open-ended polymers such as crystals, sheets, or filaments (Oosawa and Asakura, 1975; Prevelige et al., 1993).

The model developed here is completely general. For clarity, we have not incorporated many possible features into the model. Such mechanisms include autostery (Schwartz et al., 1998; Caspar, 1980), energy penalties for nucleation, cooperativity of association, cooperativity of dissociation, scaffolding proteins, protein-nucleic acid interaction, and the role of post-assembly conformation changes (c.f. Conway et al., 2001). These modifications to the model represent directions for future investigation. Divergence between predicted assembly kinetics and experimental observation is expected if one fails to account for biologically relevant reactions. It is this divergence that will allow investigators to progressively develop and refine accurate descriptions of virus assembly.

These models of assembly improve our ability to understand the behavior of free capsid protein during an infection and the effect of nucleating factors. They are readily applied to analysis of in vitro assembly reactions for virions and other spherical nanoparticles. We can speculate that some day there will be antiviral drugs directed at interfering with normal virus assembly; evaluating such antivirals will require accurate models of the assembly reactions that they inhibit.