INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION CONCLUSIONS REFERENCES

Actin filaments are key components of the cytoskeleton and play many important roles in both muscle and nonmuscle cells. The filaments are two-stranded helical polymers formed from actin monomers assembled in a polar fashion. Because of this polarity, the two ends of the filament, called the barbed and pointed ends, have different properties, both in terms of structure and dynamics. Actin filament polymerization has been extensively studied for many years and the factors controlling the kinetics have been well characterized (Pollard, 1986, 1990; Carlier, 1991). Actin polymerization follows a nucleation-elongation scheme characterized by unfavorable nucleation followed by more favorable elongation after a stable nucleus is formed. Despite the amount of time and effort that has been devoted to studying the polymerization phase, the process of spontaneous nucleation is still not well understood. Due to the size of the system and time scales involved, it is not possible in experiments to view the intermediates formed in the nucleation process, but measurements made during the polymerization phase can only be extrapolated to make estimates about the nucleation phase (Wegner and Engel, 1975; Tobacman and Korn, 1982; Cooper et al., 1983; Frieden and Goddette 1983; Frieden, 1983; Buzan and Frieden, 1996). All of these previous studies used kinetic modeling to fit polymerization curves, but we now have more advanced simulation techniques that offer us the unique opportunity to look at protein interactions at the level of the proteins involved. Using a combination of different computational methods, we hope to answer several outstanding questions about actin nucleation: what is the size of the critical nucleus, what are the steps taken in forming the critical nucleus, and what are the rate constants for each of the nucleation steps?

The study of protein-protein interactions through computational means has been well established in recent years. Brownian dynamics (BD) simulations have been shown to be very effective at both reproducing and predicting protein association rates (Nambi et al., 1991; Northrup et al., 1993; Kozack et al., 1995; Gabdoulline and Wade, 1997; Elcock et al., 1999, 2001; Sept et al., 1999). Similarly, the calculation of binding free energies allows one to estimate the contributions from many different sources, such as electrostatics, configurational entropy, hydrophobic interactions and desolvation (e.g., Sharp et al., 1991; Horton and Lewis, 1992; Simonson and Brünger, 1994; Brady et al., 1997; Gilson et al., 1997; Hummer et al., 1998; plus many more). These two types of calculations are complementary because we know that the thermodynamics and kinetics are related through the relation

(1)

where G_b is the binding free energy, K_d is the dissociation equilibrium constant, and k₊/k are the association/dissociation rate constants for a two-state binding reaction. This study will involve three separate steps. First is the identification of all protein complexes that could be formed during nucleation. It is important to note that we will make no assumptions about the size of the critical nucleus, but will investigate all possible nucleation pathways. Next, we will perform two independent sets of calculations: Brownian dynamics simulations to get the association rate constants and free energy calculations to estimate the binding free energy. Last, by combining these results with Eq. 1, we will be able to find dissociation rate constants for each of the protein complexes. By combining the nucleation kinetics with the known polymerization rates, and comparing experimental data and the polymerization curves predicted using these rate constants, we should be able to elucidate details about the nucleation process.

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION CONCLUSIONS REFERENCES

Actin structures

To perform our calculations, we first needed to define all of the protein structures that we were interested in. Our interest was in the complexes formed during nucleation, and, as such, we constructed all reasonable dimers, trimers, and tetramers that could be formed along this pathway, as illustrated in Figs. 1 and 2. Because of the helical structure of the filament, each monomer is in contact with four other neighboring monomers. This means that there are two unique dimers that can be formed: a cross-filament dimer (between red and blue monomers in Fig. 1), and a longitudinal dimer (corresponds to red-red or blue-blue dimers in Fig. 1). Although some of the complexes in Fig. 2 appear to be identical, it should be again noted that we are dealing with a polar structure and the two ends of the polymer have different properties. We were working under the assumption that nucleation occurs only through the addition of monomers (no dimer-dimer or higher-order interactions). This assumption was later validated when we saw the short lifetime of the dimer states. We used a monomer-tetramer complex to represent a monomer interacting with a longer filament (i.e., polymerization). The results obtained in the monomer-tetramer simulations were matched to experimental values and used to scale all the results for the smaller structures. As we had done previously (Sept et al., 1999), we used the actin filament structure produced by Holmes et al. (1990), because this allowed us to easily define the contacts between adjacent monomers.

View larger version (62K): [in this window] [in a new window]   FIGURE 1   A spacefilling model of the actin filament of Holmes et al. (1990) showing how the different monomers correspond to our cartoon models. The two strands of the filament are colored blue and red with the individual monomers in different shades. Note that there are two distinct dimers that can be formed: the cross filament dimer on the right (one red monomer and one blue monomer) and the longitudinal dimer on the left (either blue-blue or red-red).

View larger version (18K): [in this window] [in a new window]   FIGURE 2   Depiction of the possible pathways we considered in our model. The paths within the dotted line are the reactions that were included in the nucleation-elongation scheme. The bold arrows indicate what we found to be the preferred nucleation pathway, accounting for 99.7% of the polymer that is formed.

Binding free energy calculations

The binding free energy of forming a protein-protein complex has many different components resulting from electrostatic and van der Waals interactions, changes in internal, external and solvent entropy, hydrophobic interactions, etc. Although most of these interactions can be treated theoretically, the accuracy of these calculations was not sufficient for our needs because relatively small changes in G_b will result in large changes in K_d. Because we were always dealing with the binding of a monomer to different-sized polymers, many of the free energy contributions were the same in each case we examined (conformational changes, translational/rotational entropy), whereas others depended on the amount of surface area that is buried in forming each complex (hydration, solvent and side-chain entropy, etc.). We made the assumption that the binding free energy for each step of the nucleation and polymerization process could be given by the equation

(2)

where we have assumed that all of the contributions to the binding free energy can be grouped into three terms: G_elec is the electrostatic interaction energy, A represents the contributions proportional to the change in surface area, and G₀ is the sum of the energy contributions that are constant for each step. In each case, we calculated the electrostatic interaction energies using the University of Houston Brownian Dynamics program (Madura et al., 1995) to solve the Poisson-Boltzmann equation for both the protein complex and the structures that came together to form the complex. By subtracting the energies of the individual proteins from that of the complex, we found the interaction energy that was due to electrostatic interactions. We used the full nonlinear form of the Poisson-Boltzmann equation and assumed an ionic strength of 50 mM and a solvent dielectric of 78.4. The protein dielectric was chosen to be 12 on the basis of pK_a calculations (e.g., Antosiewicz et al., 1996; García-Morena E. et al., 1997) and should help compensate for the fact that we do not allow for protein flexibility. Using a 1.4-Å-radius probe, we calculated the solvent-accessible surface again for both the complex and its parts to determine the surface area A that was buried in each case. The buried surface is treated uniformly, and we ignore and details such as the curvature of the surface or the hydrophobic or other natures of the residues that make up the interface.

Although Eq. 2 still contained two unknowns,  and G₀, we had the additional constraint that the monomer-tetramer interactions (reaction m in Fig. 2 and Table 1) matched experimental results. For Mg-ATP actin at an ionic strength of 50 mM, the rate constants for association and dissociation are 11.6 µM¹ s¹ and 1.4 s¹, respectively (Pollard, 1986), which tells us that G_b = 9.49 kcal/mol. Taking the values for reaction m in Table 1 and inserting them in Eq. 2, we see that our scaling relation has the form,

(3)

Because we have one equation with two unknowns, we still have the freedom to choose one of our parameters,  or G₀. We will select G₀ as our degree of freedom, but our method of choosing a specific value for this variable will be discussed later.

Brownian dynamics simulations

The Brownian dynamic (BD) simulations were performed using the program SDA (Gabdoulline and Wade, 1997, 1998), as done in previous actin polymerization simulations (Sept et al., 1999). The basis of BD simulations is the solution of the equation (Ermak and McCammon, 1978),

(4)

where R is the position of the protein, D is the diffusion constant (translational or rotational), t is time step, k is Boltzmann's constant, and T is the temperature. The relative position of the proteins is affected by two parameters: F, interaction forces between the proteins, and S, a stochastic term that captures the Brownian motion caused by solvent interactions. The electrostatic calculations for each protein complex were done using the same electrostatic parameters as in the free energy calculations, and the same binding criteria were used for all simulations (three independent contacts formed at 10-Å separation). We needed to perform simulations for the association of two monomers (four possible binding sites), the binding of a dimer and a monomer (two different dimers each with multiple binding sites), and, finally, a trimer with a monomer (again three possible configurations). Apart from different binding contacts, the parameters of each simulation were the same except for the rotational and translational diffusion constants. These will obviously vary with the size and shape of the molecule and were set for each case of a monomer, dimer, or trimer using formulae developed for the diffusion of ellipsoids (Bereolos et al., 1993). The diffusion constants are given in Table 2. For each possible association reaction, we performed about 100,000 trajectories where the proteins were started at a separation of b = 120 Å, and each trajectory was halted if the binding criteria were satisfied or the proteins escaped beyond a separation q = 500 Å. From the fraction of trajectories that satisfied the binding criteria, we calculated the association rate constant k₊ as (Northrup et al., 1984)

(5)

where k_D(b) is the rate at which molecules arrive at a separation b, and  is the fraction of trajectories that form a successful protein complex. Because the interaction potential was negligible beyond 120 Å, the rates k_D(x) could simply be replaced by the Smoluchowski rate 4Dx (Smoluchowski, 1916). Once we had the relative association rates for each possible nucleation step, we scaled them all by the same factor such that the rate for the monomer-tetramer simulation (i.e., polymerization) matched the elongation rate of 11.6 µM¹ s¹ measured in experiments.

Nucleation-elongation equations

Using the association and dissociation rates for the complete nucleation process, we can solve a set of nucleation-elongation equations to get the time course of polymerization. For a given choice of G₀, we will get rate constants for every pathway depicted in Fig. 2, however, because many of the complexes are extremely unfavorable, including all these possibilities needlessly complicates the set of equations we need to solve. From the rates shown in Table 1, we determined that it was only reasonable to include the structures within the dotted line in Fig. 2. This means we have monomers, two possible dimers, one trimer and one tetramer, and our set of equations looks like:

(6)

where A represents actin monomers, A₂ and A^*₂ are the two possible dimers, A₃ and A₄ represent trimers and tetramers, and F represent all filaments longer than 4 monomers. The letters for the various reactions correspond to the reactions in Fig. 2 and Table 1. So we do not need to track the distribution of filaments represented by F, we have assumed the back reaction rate to be zero for the last "nucleation" reaction. This assumption will not affect the time course of polymerization and greatly simplifies solving the resulting set of differential equations.

	RESULTS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION CONCLUSIONS REFERENCES

Table 1 summarizes all of the simulation data and calculated binding energies for each structure along the nucleation pathway. The structures for each complex were derived from the actin filament structure of Holmes et al. (1990) (see Fig. 1). The association rates from the BD simulations were uniformly scaled so that the rate for the monomer-tetramer system (reaction m) was equal to the experimentally measured rate of 11.6 µM¹ s¹. The unscaled rate constant for reaction m was 37.8 µM¹ s¹ and thus the scaling factor was about 0.307.

For each nucleation and polymerization step, G_elec and A were found using University of Houston Brownian dynamics program (Madura et al., 1995). For a given choice of G₀, the corresponding value of  is determined by our scaling relation (Eq. 3), and by inserting the values for G_elec, A, G₀, and  in Eq. 2, we can determine G_b for each reaction in Fig. 2. The dissociation rate constants were then found by inserting the values for G_b and k₊ into Eq. 1. Figure 3 shows the effect of different values of G₀ on the predicted time course of polymerization for the same G-actin concentration. To produce these plots, the rate constants resulting from each choice of G₀ were inserted into the kinetic scheme given in Eq. 6. A smaller value of G₀ results in more nucleation and faster polymerization due to a higher concentration of filament ends, while larger values inhibits the nucleation process. A true test of our model is to compare the predictions of these simulations with experimental results. Figure 4 shows both the experimental and simulated polymerization curves for five different G-actin concentrations between 3 and 10 µM. Through a trial and error procedure, we found a value of G₀ = 20.5 kcal/mol resulted in the best fit of the polymerization data, and this choice of G₀ resulted in a value for  of 10.9 cal/mol/Ų. It should be noted that G₀ was the only free parameter to simultaneously fit all six curves.

View larger version (20K): [in this window] [in a new window]   FIGURE 3   Polymerization plots showing the effect of different values of G0 (in kcal/mol) on the kinetics of nucleation and the time course of polymerization. The corresponding  values (in cal/mol/Å2) were calculated from Eq. 3. All curves were calculated for an G-actin concentration of 5 µM.

View larger version (27K): [in this window] [in a new window]   FIGURE 4   Plots of the predicted time course of polymerization using the rate constants given in Table 1. The actin concentrations in each plot are (bottom to top) 3, 4, 5, 6, 8, and 10 µM. The experimental data is courtesy of Dr. Harry Higgs (Salk Institute).

	DISCUSSION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION CONCLUSIONS REFERENCES

The only measure that we have to validate the predictions of our model is to compare our simulated polymerization results with corresponding experimental findings, as shown in Fig. 4. Even though we have only one free parameter to fit all of the polymerization curves, our results match very well over the complete range of actin concentrations.

Nucleation pathway and the critical nucleus

Our assumption that nucleation proceeds via monomer addition appears to be valid based on the rate constants that we find. Due to the large dissociation rate constant and low concentration of dimers, the chance of these structures coming together to form a tetramer is extremely unlikely. Not surprisingly, the most critical step in the nucleation process is the formation of the dimer. For completeness, we investigated all possible trimers that could be formed by monomer addition to the two dimers (see Fig. 2). The only probable trimer that resulted was through reactions e and f, and, although we present the results for the other possible pathways (c, d, g, h, and i), they will be ignored in the subsequent discussion.

The key point in the predicted nucleation pathway is that the longitudinal dimer (b) is more favorably formed than the cross-filament dimer (a). Although the difference in G_elec is relatively small between the two dimers, the amount of buried surface area differs significantly. Because of the additional 400 Å² buried by the longitudinal dimer (b), its binding energy becomes about 5 kcal/mol more favorable than that of dimer (a), making it the dominant pathway. Figure 5 shows the time course of polymerization and the concentrations of the two dimers and the trimer. We see that the concentration of the cross-filament dimer (a) is about four orders of magnitude less than the longitudinal dimer (b), and even though step (e) is thermodynamically more favorable than step (f), the (a)-(e) pathway only contributes about 0.3% of the trimers that are formed. Although we have included both pathways in our nucleation-elongation scheme, the dominant pathway is (b)-(f)-(k)-(m)-, as indicated by the bold arrows in Fig. 2, and a kinetic scheme using only these reactions is indistinguishable from the results we present here.

View larger version (21K): [in this window] [in a new window]   FIGURE 5   The simulated time course of polymerization for 5 µM actin showing the relative concentrations of the two dimers formed by reactions a and b, and the total trimer concentration formed by reactions e and f. The amount of trimer resulting from the a-e pathway only accounts for about 0.3% of the total trimer concentration.

It may not be immediately obvious, but there are significant differences that arise depending on the dimer that is formed. If we imagine that the preferred dimer would be the cross-filament dimer (a), the next step in forming the trimer (e) is nearly identical to subsequent polymerization steps (k) and (m). There are two reasons for this. First, the surface area buried in each of the steps (e), (k), and (m) is identical in that the only difference between these steps are additional monomers at the end opposite of where the binding is occurring. Second, the only difference in the electrostatic interactions in each step is again the interaction between the new monomer and the monomers opposite the binding end. From Table 1, we see that the values of G_elec are very similar for these three steps. Hence, the kinetics resulting from this pathway would have only one unfavorable step in the nucleation process, and growth beyond the dimer (a) would basically follow polymerization kinetics, implying that the dimer (a) would be the critical nucleus. Despite the many attempts in the past by us and other researchers, it is not possible to fit polymerization data using a kinetic model with the dimer as the critical nucleus (data not shown). The main problem appears to be that, with only one nucleation step, the variation in the rate of nucleation with concentration is not enough to give an ample spread in the polymerization plots.

The nucleation pathway predicted through our modeling is fundamentally different because the dimer (b) is formed between monomers within the same protofilament, and the formation of the trimer (f) is also different from further polymerization steps. This results in two nucleation steps, the formation of the dimer (b) (unfavorable with a K_d = 4.6 M) and the trimer (still less favorable than polymerization, with a K_d = 0.6 mM). Beyond the trimer, however, the association and dissociation rate constants are essentially equal to the polymerization values, indicating that, in this case, the trimer is the critical nucleus. The polymerization plots that result from these kinetic rates agree very well with the experimental curves for a choice of G₀ = 20.5 kcal/mol (Fig. 4). Because we have one more nucleation step in this pathway, the overall nucleation rate has a stronger dependence on the monomer concentration and the resulting plots have a wider separation for the same concentrations.

The rate constants and equilibrium constants that we arrived at are in fairly good agreement with previous estimates from basic kinetic modeling (Wegner and Engel, 1975; Tobacman and Korn, 1982; Frieden and Goddette 1983; Frieden, 1983; Buzan and Frieden, 1996). All of these studies had different nucleation schemes, and, in some cases, it had to be assumed that the rate constants for each of the nucleation steps were identical, but still the general conclusion was that the critical nucleus size was a trimer (summarized in Cooper et al., 1983). Frieden (1983) used different rate constants for the nucleation steps and found equilibrium constants of 0.8 M and 5 µM for the formation of the dimer and trimer, respectively. Based on experimental differences in pH, ionic strength, and the type of actin used (yeast versus muscle), the deviation between these values and ours is understandable.

Interpretation of G₀ and  values

Our assumption in the equation for our binding free energy in Eq. 2, was that all of the components in the binding energy could be grouped into the three terms G_elec, A, and G₀. This is a great simplification in terms of the detail of the interactions that we are able to capture, but the results appear to support this model. We know that G_elec captures the electrostatic interactions and the effect of desolvation, but it is also possible to at least partially assign the other terms to specific contributions. The term A is intended to account for many different factors, but the main contributions are most likely the result of a combination of hydrophobic interactions and the removal of bound waters from the binding region. The value of 10.9 cal/mol/Ų for  is consistent with previous estimates that have found a wide range of values for these two interactions (e.g., Sharp et al., 1991; Horton and Lewis, 1992; Giesen et al., 1994; Simonson and Brünger, 1994; Fukunishi and Suzuki, 1996; Hermann, 1996; Hummer et al., 1998). The primary contributions to G₀ will be the loss of translational and rotational entropy of the monomer that is binding. Estimates of the translational and rotational entropy of a protein are quite variant, but the value of 20.5 kcal/mol for G₀ is certainly consistent with theoretical estimates (Erickson, 1989; Brady and Sharp, 1997; Tamura and Privalov, 1997) and experimental measurements made for actin polymerization (Kinosian et al., 1991).

Effect of the nucleotide and divalent cation

There are differences in both the nucleation and polymerization properties of actin depending on the nucleotide (ATP, ADP, or no nucleotide) or metal ion (Ca²⁺ or Mg²⁺) that is bound (e.g., Estes et al., 1992). These differences almost certainly arise from changes in the conformation or dynamics that affect the interaction between the monomers (Moraczewska et al., 1999), but we have very little structural information that we can use to support this theory. By changing the bound nucleotide or cation, or even altering environmental conditions of such pH or ionic strength, we would arrive at different rate constants for the nucleation and polymerization steps. However, based on the structural arguments presented earlier, it seems unlikely that the nucleus size could ever be larger than a trimer. It is feasible, however, that cation and nucleotide changes could lead to the cross-filament dimer (a) being more favorably formed. This could introduce another nucleation pathway and possibly decrease the effective size of the critical nucleus. Recent experiments showing a decrease in the lag phase of nucleotide-free actin polymerization are but one possible demonstration of this effect (De La Cruz et al., 2000), but more structural information is needed.

Implications for nucleation within the cell

The nucleation of actin filaments in vivo is of utmost importance because this is the only method the cell has of controlling when and where actin filaments are formed. Spontaneous nucleation may not play a large role in the cell, but actin polymerization is often triggered by some other nucleating factor. Recently, significant interest has been directed toward the study of the Arp 2/3 complex and its ability to initiate filament assembly. It is most tempting to think that the complex of Arp 2 and Arp 3 would mimic an actin dimer, thereby removing the most unfavorable nucleation step, and, by simply binding one actin monomer, a stable nucleus could be formed. Studies using purified Arp 2/3 complex (Mullins et al., 1998) do not support this notion, but, when combined with other proteins from the WASp/Scar family, Arp 2/3 complex significantly increases the amount of nucleation (Higgs and Pollard, 2000; Higgs et al., 1999; Machesky et al., 1999; Rohatgi et al., 1999; Winter et al., 1999; Yarar et al., 1999). Another study involving Arp 2/3 complex and ActA does appear to remove the lag phase of polymerization (Welch et al., 1998), but the structural details of this mechanism again are not known. If Arp 2/3 complex does not mimic an actin dimer but instead simply stabilizes the nucleus as it is formed (by reducing the k for one or both of the nucleation steps), this would explain its ability to promote polymerization without the complete removal of the lag phase. This is an area that obviously requires much further investigation.

Limitations of the model

Although the results of our model appear to agree very well with experimental results, there are several details about the methods that need to be pointed out. The Brownian dynamics simulations assume that the formation of each protein-protein complex is controlled by diffusion and electrostatic interactions. We know this to be the case for barbed-end actin polymerization, but here this assumption also applies to the nucleation phase. For the binding free energy calculations, we assumed that we could represent the energies using Eq. 2. This is admittedly a simplified representation, but it captures the essential components: electrostatic and hydrophobic interactions, desolvation and configurational entropy. It also has the advantage that it introduces only one free parameter into the model because it is constrained by the binding energy for the polymerization step. Including more terms in the energy expansion (e.g., van der Waals, polar and apolar contributions), could increase the accuracy of our free energy calculations, but it would also introduce additional free parameters in our model, which do not appear to be required. We are also limited by the fact that we must deal with rigid protein structures. There is no doubt that conformational changes occur during the nucleation and polymerization phases, but, currently, we have no information about the differences between G-actin and F-actin structures, or how these structures compare with the actin-DNaseI structure used by Holmes et al. (1990). A recently reported structure for G-actin (R. Dominguez, Boston Biomedical Research Institute, personal communication) should offer new insights into how important these difference are in actin nucleation and polymerization.

	CONCLUSIONS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION CONCLUSIONS REFERENCES

Our study of the spontaneous nucleation of actin filaments leads us to the conclusion that the trimer is the critical nucleus size. Through the combination of BD simulations and free energy calculations, we were able to estimate the kinetic rate constants for each of the nucleation steps by scaling with known values for actin polymerization. The predicted time course of polymerization arising from these rate constants agrees very well with experimental results over a range of actin monomer concentrations. Future work combining such calculations with additional factors, such as the Arp 2/3 complex and other associated proteins, could give more insight into nucleation and polymerization within the cell.