# Metastability and nucleation in the dilute fluid phase of a simple model of globular proteins

## 1 Introduction

Metastability is the persistence for a long time of a phase which is not the equilibrium phase. It can be both a blessing and a curse. In protein solutions it is a curse. Protein crystals are required for X-ray crystallography to determine their full structure  but protein solutions at concentrations well in excess of the solubility of the crystalline phase are often stable essentially indefinitely; the rate of nucleation of the crystalline phase is essentially zero. Here, we consider a crude model of a globular protein and we find that depending on the parameters of the model, the dilute fluid phase may be stable indefinitely with respect to crystallisation. If the solution is cooled at some low density, it is stable with respect to crystallisation down to temperatures at which the solution undergoes a fluid-fluid transition. This transition has been observed in protein solutions . The agreement in the phase and nucleation behaviour between the simple model and experiment is encouraging. It is clear what underlies the behaviour of the model and we may hope that similar physics underlies the behaviour of globular proteins. The fluid is metastable for a wide range of parameters and temperatures because as the attractions are directional and short-ranged the crystal is only stable when these attractions are strong relative to the thermal energy $`kT`$. These strong attractions mean that the interfacial tension between the dilute fluid and crystalline phases is high and it is this that inhibits crystallisation. This suggests that to increase the nucleation rate the attractions should be modified to become more like that in argon or other simple atoms and molecules, i.e., to become less anisotropic and longer ranged. The phase diagram will correspondingly become more like that of simple atomic fluids.

The interactions between globular proteins are rather poorly understood but it seems clear that many of the attractive interactions are directional and quite short-ranged ; two protein molecules must not only be close to each other to attract each other but they must also be correctly oriented. An example is the attraction between hydrophobic patches on the surfaces of globular proteins; only if the proteins are oriented so that these parts of their surfaces face each other is there an attraction. So, our model, specified in section 2, contains directional attractions; in fact for simplicity it contains only directional attractions. The model was introduced by us and its bulk phase behaviour calculated in Ref. . Here, we carefully define metastability and derive an approximate theory to tell us when the dilute fluid is metastable and when nucleation occurs. We then present and discuss results, and finish with a conclusion.

## 2 Model

Our model is exactly the same as in Ref. . The potential is a pair potential $`\varphi `$ which is a sum of two parts: a hard-sphere repulsion, $`\varphi _{hs}`$, and a set of sites which mediate short-range, directional attractions. There are $`n_s`$ sites, where $`n_s`$ is an even integer. In order to keep the model as simple as possible there are no isotropic attractions and all the directional attractions are of the same strength. The sites come in pairs: a site on one particle binds only to the other site of the pair on another particle. The two sites of a pair are numbered consecutively so that an odd-numbered site, $`i`$, binds only to the even-numbered site, $`i+1`$. This is the only interaction between the sites, an odd-numbered site, $`i`$, does not interact at all with sites other than the $`(i+1)`$th site. The orientation of site number $`i`$ is specified by means of a unit vector $`𝐮_i`$. We can write the interaction potential between a pair of particles as

$$\varphi (r_{12},\mathrm{\Omega }_1,\mathrm{\Omega }_2)=\varphi _{hs}(r_{12})+\underset{i}{\overset{^{}}{}}\left[\varphi _{ii+1}(r_{12},\mathrm{\Omega }_1,\mathrm{\Omega }_2)+\varphi _{ii+1}(r_{12},\mathrm{\Omega }_2,\mathrm{\Omega }_1)\right],$$
(1)   
where the dash on the first sum denotes that it is restricted to odd values of $`i`$. The interactions between the sites on the two particles are $`\varphi _{ii+1}(r_{12},\mathrm{\Omega }_1,\mathrm{\Omega }_2)`$, which is the interaction between site $`i`$ on particle 1 and site $`(i+1)`$ on particle 2, and $`\varphi _{ii+1}(r_{12},\mathrm{\Omega }_2,\mathrm{\Omega }_1)`$, which is the interaction between site $`i`$ on particle 2 and site $`(i+1)`$ on particle 1. These are functions of $`r_{12}`$, $`\mathrm{\Omega }_1`$ and $`\mathrm{\Omega }_2`$, which are the scalar distance between the centres of particles 1 and 2, the orientation of particle 1 and the orientation of particle 2, respectively. The particle is rigid, but not axially symmetric, so its position is completely specified by the position of its centre and its orientation $`\mathrm{\Omega }`$, which may be expressed in terms of the three Euler angles.

The hard-sphere potential, $`\varphi _{hs}`$, is given by

$$\varphi _{hs}(r)=\{\begin{array}{cc}\mathrm{}\hfill & r\sigma \hfill \\ 0\hfill & r>\sigma \hfill \end{array},$$
(2)   
where $`\sigma `$ is the hard-sphere diameter. The conical-site interaction potential $`\varphi _{ii+1}`$ is given by 

$$\varphi _{ii+1}(r_{12},\mathrm{\Omega }_1,\mathrm{\Omega }_2)=\{\begin{array}{cc}ϵ\hfill & r_{12}r_c\text{and}\theta _{1i}\theta _c\text{and}\theta _{2i+1}\theta _c\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array},$$
(3)   
where $`\theta _{1i}`$ is the angle between a line joining the centres of the two particles and the unit vector $`𝐮_i`$ of particle 1, and $`\theta _{2i+1}`$ is the angle between a line joining the centres of the two particles and the unit vector $`𝐮_{i+1}`$ of particle 2. The conical-site potential depends on two parameters: the range, $`r_c`$, and the maximum angle at which a bond is formed, $`\theta _c`$. Of course, as the attractions are directional, $`\theta _c`$ will be small, no more than about $`30^{}`$. The attractions are also short ranged, $`r_c`$ no more than 10% larger than $`\sigma `$.

The angles between the site orientations, the vectors $`𝐮_i`$ will determine which crystal lattice is formed. For simplicity, we will take the sites to be arranged such that they are compatible with a simple cubic lattice. Then if we express the unit vectors $`𝐮_i`$ in Cartesian coordinates, $`(x,y,z)`$, then when we have four sites, $`n_s=4`$, the set of vectors $`𝐮_1=(1,0,0)`$, $`𝐮_2=(1,0,0)`$, $`𝐮_3=(0,1,0)`$ and $`𝐮_4=(0,1,0)`$ would describe our model. For six sites then we add two additional sites at orientations $`𝐮_5=(0,0,1)`$ and $`𝐮_6=(0,0,1)`$.

Later on when we discuss the metastable fluid we will discuss rates. In order to do this we let $`\tau `$ be the characteristic time for the dynamics in the dilute fluid. We will not need to specify $`\tau `$ exactly but it is of order the time a molecule takes to diffuse the average separation between the molecules.

## 3 Theory for the bulk phases

The free energies of the fluid and bulk phases were derived in our previous paper, Ref. , so we only outline their derivations here. See Ref.  for details.

### 3.1 Theory for the fluid phase

The theory for the fluid phase of particles interacting via a hard-core and directional attractions mediated by sites is well established . Our theory is based on the generalisation of Chapman, Jackson and Gubbins  of Wertheim’s perturbation theory . The perturbation theory gives for the Helmholtz free energy per particle of the fluid phase, $`a_f`$, 

$$\beta a_f(\eta ,T)=\beta a_{hs}(\eta )+n_s\beta \mathrm{\Delta }a(\eta ,T),$$
(4)   
where $`a_{hs}`$ is the Helmholtz free energy per particle of a fluid of hard spheres, and $`\mathrm{\Delta }a`$ is the change in free energy per bonding site due to bonding,

$$\beta \mathrm{\Delta }a=\mathrm{ln}X+\frac{1}{2}(1X).$$
(5)   
We use an accurate expression derived from the equation for the pressure of Carnahan and Starling  for $`a_{hs}`$. The volume fraction $`\eta =(N/V)(\pi /6)\sigma ^3`$ is a reduced density, it is the fraction of the solution’s volume occupied by the molecules. $`N`$ and $`V`$ are the number of molecules and the volume, respectively. $`\beta =1/kT`$, where $`k`$ is Boltzmann’s constant and $`T`$ is the temperature. $`X`$ is the fraction of sites which are not bonded to another site. As all site-site interactions are equivalent the fraction of each type of site which is not bonded is the same. The fraction of sites which are bonded and the fraction which are not bonded must, of course, add up to one. Thus we can simply write down a mass-action equation for $`X`$, 

$$1=X+\rho X^2Kg_{hs}^c(\eta )\mathrm{exp}(\beta ϵ),$$
(6)   
where $`g_{hs}^c`$ is the contact value of pair distribution function of a fluid of hard spheres, and $`\rho =(N/V)\sigma ^3`$. The volume of phase space (both translational and orientational coordinates) over which a bond exists is $`K`$ ,

$$K=\pi \sigma ^2(r_c\sigma )(1\mathrm{cos}\theta _c)^2.$$
(7)   
The mass-action equation, Eq. (6), is a quadratic equation for $`X`$ and it can be solved for $`X`$. Inserting this solution in Eq. (5) and then the result into Eq. (4) yields the Helmholtz free energy as a function of density and temperature. The state of our single component fluid is specified by the ratio of the site energy to the thermal energy, $`\beta ϵ`$, and the volume fraction, $`\eta `$. Note that Eq. (6) is not quite the same as the equivalent equations in Refs. . In those references $`\mathrm{exp}(\beta ϵ)`$ is replaced by $`\mathrm{exp}(\beta ϵ)1`$. As $`\beta ϵ`$ is quite large, five or more, the difference between the two is very small. Also, our $`K`$ is $`4\pi `$ times the $`K_{AB}`$ of Ref. .

The second virial coefficient $`B_2`$ was obtained in Ref. . It is

$$B_2=B_2^{hs}\frac{n_s}{2}K\mathrm{exp}(\beta ϵ),$$
(8)   
where $`B_2^{hs}=(2\pi /3)\sigma ^3`$ is the second virial coefficient of hard spheres.

### 3.2 Theory for the crystalline phase

At low temperature, crystallisation is driven by the attractive interactions, not packing effects as it is with hard spheres. In Ref.  we used a cell theory to describe the free energy of the crystalline phase of our model . The theory is a low-temperature theory and we will use it only at low temperatures. Vega and Monson  used a cell theory to describe the solid phase of a very similar model, a simple model of water. They avoid a couple of the approximations used here at the cost of not having an analytical free energy. Within a cell theory for a solid phase, the Helmholtz free energy per particle, $`a_s`$, is given by

$$\beta a_s(\eta ,T)=\mathrm{ln}q_P$$
(9)   
where $`q_P`$ is the partition function of a single particle trapped in a cage formed by the requirements that all its $`n_s`$ sites bond to neighbouring particles, and that its hard core not overlap with any of these neighbours. If the lattice constant is $`b`$, then the particle can move a distance distance $`b\sigma `$ in the direction of any of its neighbours, without overlapping with the neighbour. In order for the bonds to not be broken the particle must always be within $`r_c`$ of the surrounding particles. This fixes the lattice constant, $`b`$, at a little less than $`r_c`$. It is a little less as when the particle moves off the lattice site it will be moving towards some of its neighbours and away from others. Thus it can explore regions where it is further than $`b`$ from some of its neighbours. The exact value of the maximum lattice constant for which the particle can move about, constrained by the hard-sphere interactions, without breaking any bond, is difficult to estimate; as is the volume available to the centre of mass of the particle . Therefore, we approximate the lattice constant $`b`$ by $`r_c`$ and the volume to which the particle is restricted by $`(r_c\sigma )^3`$. The requirement that no bonds be broken also severely restricts the orientations of the particle. When a non-axially symmetric particle is free to rotate it explores an angular phase space of $`8\pi ^2`$. However, in the crystal its rotations will be restricted to those which are small enough not to violate the requirement that the orientations of its site vectors are within $`\theta _c`$ of the lines joining the centre of the particle with the those of the neighbouring particles. Again the exact value of angular space available to the particle is complex, and it also depends on the position of the particle. We approximate this angular space by assuming that each of the three angular degrees of freedom can vary independently over a range of $`2\theta _c`$. The normalised angular space available to a particle in the solid phase is then $`(2\theta _c)^3/8\pi ^2=\theta _c^3/\pi ^2`$. The energy per particle is, of course, $`(n_s/2)ϵ`$, and so the partition function, $`q_P`$, is then just the volume available to the centre of mass of the particle times the angular space available times $`\mathrm{\Lambda }^1\mathrm{exp}[(n_s/2)\beta ϵ]`$, where $`\mathrm{\Lambda }^1`$ is the integral over the momentum degrees of freedom. Thus, we have for $`q_P`$,

$$q_P=v_P\mathrm{\Lambda }^1\mathrm{exp}\left(\frac{n_s}{2}\beta ϵ\right),$$
(10)   
where

$$v_P=(r_c\sigma )^3\left(\frac{\theta _c^3}{\pi ^2}\right).$$
(11)   
Inserting Eq. (10) for $`q_P`$ into Eq. (9),

$$\beta a_s=\mathrm{ln}\left(v_P/\mathrm{\Lambda }\right)\frac{n_s}{2}\beta ϵ=\beta \mu _s.$$
(12)   
This is the free energy at a lattice constant of $`r_c`$. The maximum possible density of a simple-cubic lattice is when the lattice constant $`b=\sigma `$, then the density is $`\sigma ^3`$. This density corresponds to a volume fraction $`\eta =\pi /6`$. When the lattice constant is $`r_c`$, the density is $`r_c^3`$ and the volume fraction is $`(\pi /6)(\sigma /r_c)^3`$.

We are interested in finding coexistence between the crystal phase and the fluid phase at low temperature, when our assumption that no bonds are broken in the solid phase will be accurate. Then the pressure at coexistence will be low and the solid will be near its minimum possible density, $`r_c^3`$. The chemical potential $`\mu _s=a_s+p_s/\rho `$ where $`p_s`$ is the pressure and $`\rho `$ is the density. At low pressure $`p_s/\rho `$ contributes a negligible amount to the chemical potential, which enables us to equate $`a_s`$ and $`\mu _s`$ as we have done in Eq. (12). The coexisting fluid density at the fluid-solid transition is then found by equating the chemical potentials in the two phases. The density of the coexisting solid phase, when the temperature is low enough that solidification is driven by the attractive interactions not packing effects, is assumed constant at $`r_c^3`$. See Ref.  for details.

## 4 Crystalline clusters

We derive a simple but rather crude approximation for the equilibrium density of crystalline clusters in a dilute fluid. The approximations used are similar in spirit to our calculation of the interfacial tension between the crystal and dilute fluid phases of the spheres with a short-range isotropic attraction . We will assume that the interface between the cluster and the surrounding dilute fluid is sharp and that the interaction between the crystalline cluster and the surrounding fluid is weak. Both these assumptions are reasonable if the fluid is dilute but not if it is dense or near a fluid-fluid critical point . Thus we will only be able to predict the densities of crystalline clusters and therefore the nucleation rate of the crystalline phase in the dilute fluid.

We require the density of crystalline clusters of $`n`$ particles, $`\rho _c(n)`$, in a dilute gas. To find this we start from the $`n`$-particle distribution function, $`\rho ^{(n)}(1\mathrm{}n)`$ in the grand-canonical ensemble 

$$\rho ^{(n)}(1\mathrm{}n)=\frac{_N\frac{z^N}{(Nn)!}\mathrm{d}(n+1)\mathrm{}\mathrm{d}(N)\mathrm{exp}\left(\beta U\right)}{_N\frac{z^N}{N!}\mathrm{d}(1)\mathrm{}\mathrm{d}(N)\mathrm{exp}\left(\beta U\right)},$$
(13)   
where $`(i)`$ is a compact form for the positional, $`𝐫_i`$, and orientational, $`\mathrm{\Omega }_i`$, coordinates of molecule $`i`$, $`(1\mathrm{}n)`$ indicates that $`\rho ^{(n)}`$ is a function of the set of $`n`$ coordinates of the $`n`$ molecules. $`U`$ is the total energy of the fluid and depends on all $`N`$ coordinates.

$$z=\mathrm{\Lambda }^1\mathrm{exp}(\beta \mu )$$
(14)   
is the activity.

Equation (13) gives the density of an $`n`$-tuple of particles with coordinates $`(1\mathrm{}n)`$ in the fluid. We want the density of an $`n`$-tuple of molecules which are in a configuration which is compatible with the $`n`$ molecules being part of a single compact crystalline cluster. Therefore we integrate over all the positions of the $`n`$ particles which are consistent with the $`n`$ particles forming a crystalline cluster, and over no other positions. Integration over all $`n`$ coordinates will give us the total number of crystalline clusters, to obtain the number density $`\rho _c(n)`$ (here $`(n)`$ indicates the dependence of $`\rho _c`$ on $`n`$ the number of molecules in the cluster, not that $`\rho _c`$ depends on the coordinates of the $`n`$th molecule) we divide by the volume,

$$\rho _c(n)=\frac{1}{n!V}^{}\mathrm{d}(1)\mathrm{}\mathrm{d}(n)\rho ^{(n)}(1\mathrm{}n),$$
(15)   
where the dash on the integration sign indicates that the integration is restricted to those configurations of the $`n`$ particles which are consistent with them forming a cluster. The factor of $`1/n!`$ is present because the particles are indistinguishable and so the integral integrates over configurations which differ only by the exchange of indistinguishable particles.

As we are assuming that the cluster is in an ideal gas Eq. (13) simplifies as we set the energy of interaction to be zero except for the energy of interaction between the $`n`$ particles in the cluster. Then the integral in the denominator of Eq. (13) is simply $`V^N`$ and that in the numerator is $`V^{Nn}\mathrm{exp}(\beta u(1\mathrm{}n))`$, where $`u(1\mathrm{}n)`$ is the energy of interaction of $`n`$ molecules. So Eq. (13) becomes

$$\rho ^{(n)}(1\mathrm{}n)=\frac{_N\frac{z^NV^{Nn}}{(Nn)!}\mathrm{exp}\left(\beta u(1\mathrm{}n)\right)}{_N\frac{z^NV^N}{N!}}.$$
(16)   
Substituting this in Eq. (15),

$$\rho _c(n)=\frac{_N\frac{z^NV^{Nn}}{(Nn)!}^{}\mathrm{d}(1)\mathrm{}\mathrm{d}(n)\mathrm{exp}\left(\beta u(1\mathrm{}n)\right)}{n!V_N\frac{z^NV^N}{N!}}.$$
(17)   
We can take $`z^n`$ times the integral out of the sum in the numerator leaving the sum in the numerator identical to that in the denominator. They cancel leaving

$$\rho _c(n)=\frac{z^n}{Vn!}^{}\mathrm{d}(1)\mathrm{}\mathrm{d}(n)\mathrm{exp}\left(\beta u(1\mathrm{}n)\right).$$
(18)   
The density of crystalline clusters of $`n`$ molecules in an ideal gas is simply $`z^n/Vn!`$ times the configurational integral of $`n`$ molecules in a cluster.

As in the cell theory for a bulk crystal we factorise the integration of Eq. (18) into a product of $`n`$ integrals and delete the factor of $`1/n!`$ as once the molecules are restricted to lie in cells they are distinguishable. Now, one of the $`n`$ integrations is over the whole volume $`V`$ of the fluid, the other $`(n1)`$ are just over the rattling motion as in the bulk and they each give a factor of $`v_P`$. The energy is taken to be the ground state energy as in the bulk and so is $`nn_sϵ/2+u_s(n)`$ where $`u_s`$ is the increase in energy due to broken bonds at the surface of the cluster. So, we have that Eq. (18) becomes

$$\rho _c(n)=z^nv_P^{n1}\mathrm{exp}\left[\frac{nn_s}{2}\beta ϵ\beta u_s\right].$$
(19)   
The spheres at the faces of the cluster do not interact with the full $`n_s`$ other spheres and this increases the energy of a cluster. If we assume that the cluster of $`n`$ molecules is cubic then it has 6 faces, each of area $`n^{2/3}\sigma ^2`$, i.e., with $`n^{2/3}`$ molecules in each face. For $`n_s=6`$ there are sites pointing in all 6 directions and a sphere at any of the 6 faces but not at an edge or corner has one bond broken. So assuming that the cluster is cubic, neglecting the fact some spheres are at edges and some at corners and therefore have 2 or 3 bonds not 1 bond broken and treating $`n`$ as a continuous variable, results in the approximation that there are $`6n^{2/3}`$ bonds broken on the surface of the cluster. Each broken bond costs an energy $`ϵ/2`$ — the energy of a bond is $`ϵ`$ with $`ϵ/2`$ assigned to each of the two particles forming the bond. Thus, for $`n_s=6`$, the increase in energy $`u_s=3n^{2/3}ϵ`$. With this expression for $`u_s`$ Eq. (19) becomes

$$\rho _c(n)=z^nv_P^{n1}\mathrm{exp}\left[\frac{n_sn}{2}\beta ϵ3n^{2/3}\beta ϵ\right]n_s=6.$$
(20)   
The approximation $`u_s=3n^{2/3}ϵ`$ becomes worse as $`n`$ decreases but it is never seriously wrong. Indeed for the smallest cluster we consider, that of 8 spheres, there is cancellation of errors and there are exactly $`6\times 8^{2/3}=24`$ bonds broken. For $`n=9`$ we predict 26.0 bonds broken when in fact there 28 broken bonds but this is not a large error. For $`n_s=4`$ only 4 of the 6 faces involve broken bonds, because there are no attraction sites on 2 faces. So, instead the energy cost is only two thirds that for 6 sites and the increase in energy is $`2n^{2/3}ϵ`$. For $`n_s=4`$ or 6 the increase in energy is given by $`(n_s/2)n^{2/3}ϵ`$.

So far we have assumed that the cluster does not interact with any of the surrounding spheres. This is reasonable for a very dilute fluid. However for a fluid which is not very dilute and is at low temperature, spheres in the surrounding fluid will tend to bond to the spheres in the faces of the cluster. We can take this into account approximately by treating the sites on the faces of a crystalline cluster as if they were sites in the fluid, then for each site there is a free energy change given by Eq. (5) — which reflects the fact that it can bond to one of the surrounding spheres. The change to the configurational integral is then of course $`\mathrm{exp}(\beta \mathrm{\Delta }a)`$ per surface site. Then the configurational integral is

$$\rho _c(n)=z^nv_P^{n1}\mathrm{exp}\left[\frac{n_sn}{2}\beta ϵn_sn^{2/3}\left(\frac{\beta ϵ}{2}+\beta \mathrm{\Delta }a\right)\right],$$
(21)   
where $`X`$ in Eq. (5) for $`\mathrm{\Delta }a`$ is the same as in the surrounding fluid.

## 5 Metastability and nucleation

Consider the density of clusters $`\rho _c(n)`$ of Eq. (21). For large $`n`$, $`\rho _c`$ is dominated by the part $`(zv_P\mathrm{exp}[(n_s/2)\beta ϵ])^n`$ as the other parts vary only as the $`n^{2/3}`$ power or are constants. Using Eqs. (12) and (14), we obtain

$$zv_P\mathrm{exp}[(n_s/2)\beta ϵ]=\mathrm{exp}(\beta \mu )\mathrm{exp}(\beta \mu _s).$$
(22)   
But as we are within the fluid-crystal coexistence region the chemical potential of the crystal $`\mu _s`$, is less than that of the fluid phase, $`\mu `$. So, Eq. (22) is greater than 1 and hence $`\rho _c(n)`$ diverges as $`n\mathrm{}`$. This is actually an automatic consequence of the fact that the crystal is more stable than the fluid.

So, our Eq. (21) predicts that in the fluid there are high densities of large crystalline clusters. This is of course not what is observed in a metastable fluid. This is because our calculation of Eq. (21) assumed that the densities of all clusters were at equilibrium, whereas in a metastable fluid the system is by definition not at equilibrium. In order to describe a metastable fluid, a fluid which is out of true equilibrium, we must apply a constraint; see Refs  for definitions and discussions of the application of constraints to study metastable fluids. This constraint must eliminate the large crystalline clusters to leave us with a fluid. We choose the constraint which eliminates all clusters above a size $`n_{min}`$:

$$\rho _c(n)=0n>n_{min},$$
(23)   
where $`n_{min}`$ is defined by

$$\rho _c(n_{min})=\underset{n}{\mathrm{min}}\left\{\rho _n\right\},$$
(24)   
i.e., $`n_{min}`$ is the number of molecules in the cluster with the lowest density, as predicted by Eq. (21). So, our constrained distribution of cluster densities is

$$\rho _c(n)=\{\begin{array}{cc}z^nv_P^{n1}\mathrm{exp}\left[\frac{n_sn}{2}\beta ϵn_sn^{2/3}\left(\frac{\beta ϵ}{2}+\beta \mathrm{\Delta }a\right)\right],\hfill & nn_{min}\hfill \\ 0\hfill & n>n_{min}\hfill \end{array}.$$
(25)   
We set the constraint so as to eliminate all clusters above the size $`n_{min}`$ because this constraint is in a specific sense the least restrictive. It is the least restrictive because if we start with the constrained equilibrium distribution of clusters, which is given by Eq. (25) and then remove the constraint, i.e., allow clusters with $`n>n_{min}`$ to form, then the initial rate at which these clusters with $`n>n_{min}`$ form is minimised. This assumes that clusters only grow one molecule at a time; that a cluster with $`(n+1)`$ molecules is formed by a cluster of $`n`$ molecules adsorbing an additional molecule. This is a reasonable assumption in a dilute fluid in which the density of single molecules is much larger than the density of clusters of 2 or more molecules. With this assumption of growth one molecule at a time the initial rate at which clusters with $`n>n_{min}`$ appear is just equal to the rate at which clusters of $`n_{min}`$ molecules acquire an additional molecule to become clusters of $`(n_{min}+1)`$ molecules, which is approximately

$$\text{rate}\rho _c(n_{min})\tau ^1.$$
(26)   
Therefore, with our choice of constraint the initial rate at which the distribution of clusters changes when the constraint is removed is minimised. This is what we meant by the constraint being least restrictive. When the constraint is removed the distribution will tend towards the equilibrium one with its crystalline-cluster densities which diverge in the $`n\mathrm{}`$ limit, i.e., the fluid will crystallise. If we neglect the fact that not all the clusters with $`n_{min}`$ molecules which gain an extra molecule will grow all the way into a crystallite, then the rate of nucleation of the crystalline phase is given by Eq. (26). In view of the highly approximate nature of our theory this neglect is reasonable so Eq. (26) is our approximation for the nucleation rate. If $`\rho _c(n_{min})`$ is very small then if the constraint is removed the distribution of cluster densities will change only very slowly. Therefore the unconstrained fluid will persist for a long time, much longer than $`\tau `$, and so the unconstrained fluid phase is observable: it is metastable. However, if $`\rho _c(n_{min})`$ is not very small then as soon as the constraint is removed the unconstrained fluid starts to crystallise. The unconstrained fluid does not last long enough to be observable: it is unstable. What constitutes a very small density is of course rather arbitrary but we will try to quantify it when we discuss our results in the next section.

## 6 Results

Experiments on globular proteins have found metastable  fluid–fluid transitions , i.e., a fluid-fluid transition which lies within the fluid-solid coexistence region. The crystallisation of proteins is often slow, taking several days, which allows the protein solution to be cooled into a region of the phase diagram where the fluid phase separates into two fluid phases of differing densities. Therefore, we show phase diagrams, in Figs. 2 and 3, in which the fluid-fluid transition lies within the fluid-solid coexistence region. For other values of the parameters of the models, $`n_s`$, $`\theta _c`$ and $`r_c`$, there is a stable fluid-fluid transition . Fig. 2 is the phase diagram of a model protein with 4 sites and Fig. 3 is the phase diagram for a model with 6 sites and a larger value of $`\theta _c`$. These two models were chosen as their phase diagrams were calculated and discussed in Ref.  and they differ markedly in how deep the fluid-fluid transition is into the fluid-solid coexistence region. In Fig. 2, the fluid-fluid critical point is at a volume fraction $`\eta =0.090`$ and at reciprocal temperature $`ϵ/kT=10.24`$. We can use as a measure of how deep the fluid-fluid transition is into the fluid-solid coexistence region the ratio of the temperature at the critical point to that of a fluid of the same density which coexists with the solid. For the model of Fig. 2 fluid at a volume fraction $`\eta =0.090`$ coexists with the solid at $`ϵ/kT=8.37`$. The ratio of the temperatures is then 0.82. For the model of Fig. 3 the critical point is at $`\eta =0.154`$ and $`ϵ/kT=7.18`$. A fluid with this density coexists with the solid phase at $`ϵ/kT=4.54`$. The ratio of the temperatures is now 0.63. Note that our temperature is a reduced temperature, a dimensionless ratio $`kT/ϵ`$. We have plotted our phase diagrams as a function of $`kT/ϵ`$ but this scale is not directly related to the real temperature of a protein solution as the protein-protein interactions (which determine $`ϵ`$) vary with the temperature of the experiment.

In Figs. 2 and 3 we have shown as a dot-dashed curve an estimate of where percolation occurs in the fluid. At percolation the association of the molecules is sufficiently strong that an infinite cluster appears , that is to say that there are an infinite number of the molecules which are joined to each other via pathways of bonds. The percolation curve gives us an indication of when the density is too high or the interactions too strong for our approximation that the crystalline clusters interact weakly with the surrounding fluid to be valid. We will not use our approximation for the cluster densities, Eq. (25), beyond (i.e., to the right of) the percolation curve. See Ref.  for an introduction to percolation. If we neglect loops of bonds we obtain what is called the classical theory of percolation which predicts that percolation occurs at a fraction of bonds $`(1X_p)`$ given by 

$$1X_p=\frac{1}{n_s1}\text{or}X_p=1\frac{1}{n_s1},$$
(27)   
$`X_p`$ is the fraction of sites not bonded when percolation occurs.

Now we will use Eq. (25) to calculate the cluster densities within the dilute fluid part of the fluid-solid coexistence region of the phase diagrams in Figs. 2 and 3. For the phase diagram of Fig. 2, the 4-site model, we have calculated cluster densities in the region of the phase diagram bounded at the right by curve where percolation occurs, from below by the curve describing the density of the fluid phase which coexists with the crystal and from above by the density of the dilute fluid phase which coexists with the dense fluid phase. The approximations we used to calculate the cluster densities, $`\rho _c(n)`$, are only reasonable at low densities and away from a critical point. The region is bounded from above by the fluid-fluid coexistence curve as we expect the fluid to become unstable with respect to condensation a little inside the coexistence curve and so our calculated cluster densities are meaningless there. We expect condensation to occur only a little into the fluid-fluid coexistence region as we expect the interfacial tension between the two fluid phases will be small and therefore that nucleation of the dense fluid phase will be rapid except very near the coexistence curve.

Throughout this region the densities of crystalline clusters of all sizes $`n=8`$ and up are tiny. For example, at $`\eta =0.1`$ and $`\beta ϵ=9`$ the density of crystalline clusters of 8 spheres is $`\rho _c(8)=𝒪(10^{21}\sigma ^3)`$ and as $`n`$ increases the density rapidly decreases. So, the density of even small crystalline clusters is negligible. The nucleation rate, Eq. (26), is effectively zero and the dilute fluid phase will be stable with respect to crystallisation effectively indefinitely: it is metastable. This finding that the crystal cannot nucleate from a dilute fluid is interesting as experiments on solutions of many proteins find it difficult or impossible to find crystallisation.

The nucleation rate is so low because the nuclei, the crystalline clusters have extremely low densities. This can be traced to the interfacial term in our expression for $`\rho _c`$, Eq. (21). This is the second term in the exponential which varies as the number of molecules at the surface, as $`n^{2/3}`$. It is large because under conditions that the crystal coexists with a dilute fluid the ratio between the attraction energy and the thermal energy $`ϵ/kT`$ is large. At the surface of the cluster bonds are broken and each broken bond decreases the density of a nucleus by $`\mathrm{exp}(\beta ϵ/2)`$, which is rather large. In the language of classical nucleation theory  the barrier to nucleation is high because the surface tension is high. The surface tension $`\gamma `$ here comes from the energy of the broken bonds, $`\gamma (1/2)ϵ\sigma ^2+\mathrm{\Delta }a\sigma ^2`$, where $`\mathrm{\Delta }a`$ is small, $`𝒪(0.1kT)`$.

In view of the extremely small numbers we have not plotted cluster densities for the model parameters of Fig. 2. However, the fluid-fluid transition is deeper in the fluid-solid coexistence region in Fig. 3 so larger cluster densities are achievable. Plots of $`\rho _c(n)`$ against $`n`$ for three points in the dilute phase of Fig. 3 are shown in Fig. 4. The three points are chosen to be at roughly the highest densities at which the theory is reliable and the fluid is outside the fluid-fluid coexistence region the solid, dashed and dotted curves the supersaturations $`\beta (\mu \mu _s)`$ are 3.71, 4.77 and 5.52, respectively. An approximation to the nucleation rate is given by Eq. (26) which is proportional to the densities at the minima of the curves in Fig. 4. We can get an estimate of what the numbers mean for a protein solution. Protein molecules are a few nms in diameter so in a sample 1mm across there are of order $`10^{16}`$ protein molecules. At $`\beta ϵ=7.5`$, $`\eta =0.05`$, $`\rho _c(n_{min})=𝒪(10^{16}\sigma ^3)`$ so in a sample 1mm across we have $`𝒪(1)`$ crystallites nucleating in the sample per time $`\tau `$. Muschol and Rosenberger  estimate diffusivities for lysozyme (a well studied globular protein) of order $`10^{10}`$m<sup>2</sup>s<sup>-1</sup>. The characteristic time of the dynamics $`\tau `$ should be of order the time a protein takes to diffuse its own diameter, this time is the square of the diameter, $`10^{17}`$m<sup>2</sup> divided by the diffusion constant $`10^{10}`$m<sup>2</sup>s<sup>-1</sup>, so we have $`\tau =𝒪(10^7\text{s})`$. So we end up with the rough estimate of $`10^7`$ crystallites nucleating in the sample per second. Nucleation is therefore rapid. In common with classical nucleation theory our approximation for the nucleation rate is the calculation of a very small number and so the errors are typically large, easily several orders of magnitude . Bearing this in mind our theory can only tell us that the model parameters of Fig. 3 lie close to the dividing line between parameter values for which nucleation of crystalline phase from a dilute fluid phase is not achievable on experimental time scales and parameter values for which it is.

George and Wilson  determined the second virial coefficients of a number of globular proteins under the conditions for which they crystallised. They found that the values of the second virial coefficients lay within a small range, which they called the ‘crystallisation slot’. Using Eq. (8) for the second virial coefficient, $`B_2`$, we can determine the values of $`B_2`$ at the 3 temperatures for which we plotted the cluster densities in Fig. 4. They are $`B_2=0.21\sigma ^3`$, $`3.03\sigma ^3`$ and $`6.35\sigma ^3`$ for $`\beta ϵ=6`$, 7 and 7.5, respectively. So, although at all 3 temperatures we have (at different volume fractions) similar densities of the minimum-density cluster the second virial coefficient varies over a large range, it even changes sign. Thus, our results for the nucleation rate do not offer an explanation of George and Wilson’s finding.

## 7 Conclusion

We have studied a simple model of a globular protein molecule in solution. The phase diagram and the densities of crystalline clusters in the dilute fluid phase have been calculated. The phase diagram predicted by our bulk free energy includes fluid-fluid coexistence within the fluid-crystal coexistence region. When this fluid-fluid coexistence region is not too deep into the fluid-crystal coexistence region, as in Fig. 2, we find that the dilute fluid phase outside of the fluid-fluid coexistence region is metastable, i.e., the rate of nucleation of the crystalline phase is negligible. It is not possible to produce a crystal directly from the dilute fluid for this model. When the fluid-fluid coexistence region is deeper into the fluid-crystal coexistence region, as in Fig. 3, the nucleation rate becomes large enough to be observable within the dilute fluid.

Essentially, we defined the dilute fluid as being the fluid at densities below the percolation threshold. This means that the fluid-fluid critical point is not included in our definition of the dilute fluid. Ten Wolde and Frenkel  have shown that near a fluid-fluid critical point the interface between the crystalline nucleus and the surrounding fluid is diffuse and that this enhances the nucleation rate dramatically. The diffuse interface is very different from the sharp interface we had to assume to obtain approximations for the cluster densities and hence the nucleation rate. If we consider the (highly inaccurate) predictions of our theory near the critical points of Figs. 2 and 3, we find that $`\rho _c(n_{min})=𝒪(10^{135}\sigma ^3)`$ and $`𝒪(10^{14}\sigma ^3)`$, respectively. So, nucleation is certainly rapid near the critical point of Fig. 3. However, the density $`\rho _c(n_{min})`$ is predicted to be so low near the critical point of Fig 2, that even taking into account the very large errors in our theory we would not expect nucleation. Although the nucleation rate is enhanced by the nature of the fluid near its critical point, as the model parameters are varied to move the critical point toward the fluid-crystal coexistence curve, the rate will tend to zero. In the limit that the critical point touches the fluid-crystal coexistence curve, i.e., at the point where the fluid-fluid transition goes from being metastable to being stable, the supersaturation at the critical point tends to zero, reducing the nucleation rate to zero.

As this work has been motivated by the difficult and important problem of crystallising globular proteins it is interesting to speculate on how the model of Fig. 2 could be crystallised. The nucleation rate is far too low in the dilute fluid so in order to increase the rate the fluid must either be made more dense or the attractions strengthened. Both of these may result in equilibrium being difficult to reach with the result that fluid could become gel-like. Also, if the fluid undergoes a fluid-fluid transition its density and hence its nucleation rate jumps . There is an optimum nucleation rate to obtain good, i.e., large with few defects, crystals. Now, if there were no fluid-fluid transition then the crystalline cluster densities and hence the nucleation rate vary continuously but at condensation the densities will jump so there is a risk that the nucleation rate will jump over the optimum one making good crystals hard to obtain. Crystallisation would be facilitated if the free energy cost of the surface of the cluster, the second term in the exponential of Eq. (21), was less. If the interactions were less directional then the crystal would be stable at higher temperatures, i.e., smaller values of $`\beta ϵ`$, where the surface would have a lower free energy.