## Abstract

Numerical simulations of homogeneous crystal nucleation with a model for globular proteins with short-range attractive interactions showed that the presence of a metastable fluid-fluid critical point drastically changes the pathway for the formation of a crystal nucleus. Close to this critical point, the free-energy barrier for crystal nucleation is strongly reduced and hence, the crystal nucleation rate increases by many orders of magnitude. Because the location of the metastable critical point can be controlled by changing the composition of the solvent, the present work suggests a systematic approach to promote protein crystallization.

As a result of rapid developments in biotechnology, there has been an explosive growth in the number of proteins that can be isolated. However, the determination of the three-dimensional structures of proteins by x-ray crystallography remains a time-consuming process. One bottleneck is the difficulty of growing protein crystals good enough for analysis. In his book on this subject, McPherson wrote, “The problem of crystallization is less approachable from a classical analytical standpoint, contains a substantial component of trial and error, and draws more from the collective experience of the past century. . . . It is much like prospecting for gold” (*1*, p. 6). The experiments clearly indicate that the success of protein crystallization depends sensitively on the physical conditions of the initial solution (1, 2). It is therefore crucial to understand the physical factors that determine whether a given solution is likely to produce good crystals.

Studies have shown that not just the strength but also the range of the interactions between protein molecules is crucial for crystal nucleation. In 1994, George and Wilson (3) demonstrated that the success of crystallization appears to correlate with the value of*B*
_{2}, the second osmotic virial coefficient of the protein solution.

The second virial coefficient describes the lowest order correction to the van't Hoff law for the osmotic pressure Π:(1)where ρ is the number density of the dissolved molecules, *k*
_{B} is Boltzmann's constant, and*T* is the absolute temperature. For macromolecules,*B*
_{2} can be determined from static light-scattering experiments (4). Its value depends on the effective interaction between a pair of macromolecules in solution (5).

George and Wilson measured *B*
_{2} for a number of proteins in various solvents. They found that for those solvent conditions that are known to promote crystallization,*B*
_{2} was restricted to a narrow range of small negative values. For large positive values of*B*
_{2}, crystallization did not occur at all, whereas for large negative values of *B*
_{2}, protein aggregation rather than crystallization took place.

Rosenbaum, Zamora, and Zukoski (6, 7) established a link between the work of George and Wilson and earlier studies of the phase behavior of spherical, uncharged colloids (8-11). Since the theoretical work of Gast, Russel, and Hall (8), it has been known that the range of attraction between spherical colloids has a drastic effect on the appearance of the phase diagram. If the range of attraction is long in comparison with the diameter of the colloids, the phase diagram of the colloidal suspension resembles that of an atomic substance, such as argon. Depending on the temperature and density of the suspension, the colloids can occur in three phases (Fig.1A): a dilute colloidal fluid (analogous to the vapor phase), a dense colloidal fluid (analogous to the liquid phase), and a colloidal crystal phase. However, when the range of the attraction is reduced, the fluid-fluid critical point moves toward the triple point, where the solid coexists with the dilute and dense fluid phases. If the range of attraction is made even shorter (less than 25% of the colloid diameter), two stable phases remain, one fluid and one solid (Fig. 1B). However, the fluid-fluid coexistence curve survives in the metastable regime below the fluid-solid coexistence curve (Fig.1B). This is indeed found in experiments (11, 12) and simulations (10). This observation is relevant for solutions of globular proteins, because they often have short-range attractive interactions. A series of studies (13-15) showed that the phase diagram of a variety of proteins is of the kind shown in Fig. 1B. Moreover, the range of the effective interactions between proteins can be changed by the addition of nonadsorbing polymer (such as polyethylene glycol) (11, 16) or by changing the pH or salt concentration of the solvent (1, 2).

Rosenbaum, Zamora, and Zukoski (6, 7) observed that the conditions under which a large number of globular proteins can be made to crystallize map onto a narrow temperature range, or more precisely, a narrow range in the value of the osmotic second virial coefficient of the computed fluid-solid coexistence curve of colloids with short-range attraction (10). Several authors had already noted that a similar crystallization window exists for colloidal suspensions (17). Here our aim was to use simulation to gain insight into the physical mechanism responsible for the enhanced crystal nucleation. We found that the presence of a metastable fluid-fluid critical point is essential.

The rate-limiting step in crystal nucleation is the crossing of a free-energy barrier. In our simulations we computed the barrier for homogeneous crystal nucleation for a model globular protein. Our model system has some of the essential features needed to get a proteinlike phase diagram of the type shown in Fig. 1B—the particles repel strongly at short distances and attract at slightly larger distances.

For the interaction between the particles, we chose a suitable generalization of the Lennard-Jones potential:
(2)where σ denotes the hard-core diameter of the particles, *r* is the interparticle distance, and ɛ is the depth of the potential well. The potential in Eq. 2 provides a simplified description of the effective interaction between real proteins in solution: It accounts both for direct and for solvent-induced interactions between the globular proteins. The width of the attractive well can be adjusted by varying the parameter α. We found that for α ≈ 50, the system reproduced the phase behavior of protein solutions studied in (13-15), that is, the fluid-fluid coexistence curve was located in the metastable region ∼20 % below the equilibrium crystallization curve.

Conventional molecular dynamics simulations cannot be used to study crystal nucleation under realistic conditions. The reason is that the formation of a critical nucleus is a rare event; crystallization in real protein solutions may take days or weeks. In a simulation, the situation would be worse because the volume that can be studied by simulation is orders of magnitude smaller, and the probability of crystal formation is decreased by the same amount. Moreover, the computational cost of molecular dynamics simulations that cover more than 10^{−8} s becomes prohibitive. Hence, the problem has to be approached in a different way. In the present work, we used an approach (18) based on the Bennett-Chandler numerical scheme to compute the rate of activated processes (19). The rate Γ can be written as the product of two factors:(3)where Δ*G** is the height of the (Gibbs) free-energy barrier separating the metastable fluid from the crystal phase. The factor exp(−Δ*G**/*k*
_{B}
*T*) denotes the probability that a spontaneous fluctuation will result in the formation of a critical nucleus; it depends strongly on the degree of supercooling. The kinetic prefactor ν is a measure of the rate at which critical nuclei, once formed, transform into larger crystallites. It depends only weakly on supercooling, unless the system is close to gelation (11). Below, we discuss simulations under conditions away from the gelation curve. Hence, the variation of the nucleation rate is dominated by the variation in the barrier height Δ*G**. A rough estimate of Δ*G** is given by classical nucleation theory (20),(4)where γ is the free-energy density per unit area of the solid-fluid interface, ρ is the number density of the solid phase, and Δμ is the difference in chemical potential between the fluid and the solid. It is the thermodynamic driving force for crystallization.

To compute the free-energy barrier that the system must overcome to form a critical crystal nucleus, we used the umbrella-sampling Monte Carlo scheme of Valleau (18, 21). This is a biased Monte Carlo scheme that makes it possible to concentrate the sampling in the region of the free-energy barrier. Without the bias, this region of configuration space is effectively inaccessible. To compute the free energy as a function of the size of the incipient crystal nucleus, we must be able to identify which particles belong to the crystal nucleus. As discussed in (18), each particle can be classified as either solidlike or liquidlike by analyzing the local symmetry of its surroundings. If the distance between two solidlike particles is less than 1.5σ, they are considered to be connected and belong to the same cluster. We denote the number of solidlike particles belonging to a given crystal nucleus by *N*
_{crys}. Crystallization near the metastable fluid-fluid critical point is strongly influenced by the large density fluctuations that occur in the vicinity of such a critical point; therefore, the free-energy barrier associated with formation of dense, liquidlike droplets should also be considered. Hence, for every particle, we determined not just the structure of its surroundings but also the local density. In the same way, the size of a high-density cluster (be it solid or liquidlike) can be defined as the number of connected particles, *N*
_{ρ}, that have a significantly denser local environment than particles in the remainder of the system. In our simulations, the free-energy “landscape” of the system was determined as a function of the two coordinates*N*
_{crys} and *N*
_{ρ}. In a crystal nucleation event, the beginning is from the homogeneous liquid (*N*
_{crys} ≈ *N*
_{ρ} ≈ 0). The free energy then increases until it reaches a saddle-point (the critical nucleus). From there on, the crystal will grow spontaneously.

We first computed the phase diagram shown in Fig. 1B using the techniques described in (10) and then studied the free-energy barrier for crystal nucleation at four different points in the phase diagram: one well above the metastable critical point (*T* = 2.23 *T*
_{c}), one at*T*
_{c}, and the remaining two at 0.89*T*
_{c} and 1.09 *T*
_{c}. The degree of supercooling was chosen such that classical nucleation theory would predict the same value of Δ*G**/(*k*
_{B}
*T*) for all systems. To estimate Δ*G** from Eq. 4, we used(5)where Δ*H* is the enthalpy of melting at the coexistence temperature *T*
_{m}. Turnbull's empirical rule was used to estimate the surface free-energy γ, which states that γ is proportional to Δ*H* (20). For all points, we studied the free-energy landscape and the lowest free-energy path to the critical nucleus.

We found that away from *T*
_{c} (above and below), the path of lowest free energy is one where the increase in*N*
_{ρ} is proportional to the increase in*N*
_{crys} (Fig. 2A). Such behavior is expected if the incipient nucleus is a small crystallite. However, around *T*
_{c}, critical density fluctuations lead to a striking change in the free-energy landscape (Fig. 2B). First, the route to the critical nucleus leads through a region where *N*
_{ρ} increases while*N*
_{crys} is still essentially zero. In other words, the first step toward the critical nucleus is the formation of a liquidlike droplet. Then, beyond a certain critical size, the increase in *N*
_{ρ} is proportional to*N*
_{crys}, that is, a crystalline nucleus forms inside the liquidlike droplet.

Clearly, the presence of large density fluctuations close to a fluid-fluid critical point effects the route to crystal nucleation. But, more importantly, the nucleation barrier close to*T*
_{c} is much lower than at either higher or lower temperatures (Fig. 3). The observed reduction in Δ*G** near *T*
_{c} by ∼30*k*
_{B}
*T* corresponds to an increase the nucleation rate by a factor 10^{13}. One interpretation of this observation is that near the fluid-fluid critical point, the wetting of the crystal nucleus by a liquidlike layer results in a value of the interfacial free energy γ, and therefore of the barrier height Δ*G**, that is much lower than would be estimated on the basis of Turnbull's rule. In fact, Haas and Drenth (22) noted that the experimentally determined interfacial free energy of small protein crystals (23) is much smaller than the value predicted on the basis of their version of Turnbull's rule.

The conventional way to lower the crystal nucleation barrier is to prepare a more supersaturated solution. However, highly supersaturated solutions tend to form aggregates rather than crystals (3, 6, 7,11, 12, 15). Moreover, in such a solution, the thermodynamic driving force for crystallization (μ_{liq} − μ_{cryst}) is also enhanced. As a consequence, the crystallites that nucleate will grow rapidly and far from perfectly (2). One of the implications of our finding that the crystal nucleation barrier is reduced near *I*
_{c} is that one can selectively speed up the rate of crystal nucleation, without increasing the rate of crystal growth, or the rate at which amorphous aggregates form. This can be achieved by adjusting the solvent conditions (for instance, by the addition of nonionic polymer) and thereby changing the range of interaction, such that a metastable fluid-fluid critical point is located just below the sublimation curve.

Our description of the early stages of protein crystallization is highly simplified, yet we suggest that the mechanism for enhanced crystal nucleation described here is quite general. The phase diagram shown in Fig. 1B is likely to be the rule rather than the exception for compact macromolecules. Moreover, it occurs both in the bulk and in (quasi) two-dimensional systems (such as membranes). It is therefore tempting to speculate that nature already makes extensive use of critical density fluctuations to facilitate the formation of ordered structures.