In this work we examine thermodynamics of fluid with ‘‘molecules’’ represented by two fused hard spheres, decorated by the attractive square-well sites. Interactions between these sites are of short-range and cause association between the fused-sphere particles. The model can be used to study the non-spherical (or dimerized) proteins in solution. Thermodynamic quantities of the system are calculated using a modification of Wertheim’s thermodynamic perturbation theory and the results compared with new Monte Carlo simulations under isobaric-isothermal conditions. In particular, we are interested in the liquid-liquid phase separation in such systems. The model fluid serves to evaluate the effect of the shape of the molecules, changing from spherical to more elongated (two fused spheres) ones. The results indicate that the effect of the non-spherical shape is to reduce the critical density and temperature. This finding is consistent with experimental observations for the antibodies of non-spherical shape. \keywordsnon-spherical proteins, liquid-liquid transition, directional forces, aggregation, thermodynamic perturbation theory \pacs87.15.A-, 87.15.ad, 87.15.km, 87.15.nr, 82.60.Lf, 64.70.Ja
В цiй роботi ми дослiджуємо термодинамiку плину з ‘‘молекулами’’, представленими двома спаяними твердими сферами, якi декорованi вузлами з притягувальними потенцiалами типу квадратної ями. Взаємодiя мiж цими вузлами є короткодiюча i спричиняє асоцiацiю мiж частинками спаяних сфер. Модель може бути використана для дослiдження несферичних (чи димеризованих) протеїнiв у розчинi. Термодинамiчнi величини системи розраховуються за допомогою модифiкацiї термодинамiчної теорiї збурень Вертгайма, i результати порiвнюються з новими симуляцiями методом Монте Карло при iзобарично-iзотермiчних умовах. Зокрема, нас цiкавить фазове розшарування рiдина-рiдина в таких системах. Модельний плин використовується для оцiнки ефекту форми молекул, що змiнюється вiд сферичної до бiльш видовженої (двi спаянi сфери). Результати вказують, що ефект несферичної форми має зменшувати критичну густину i температуру. Це узгоджується з експериментальними спостереженнями для антитiл iз несферичною формою.\keywords
несферичнi протеїни, перехiд рiдина-рiдина, напрямляюча сила, агрегацiя, термодинамiчна теорiя збурень
Fluid of fused spheres]
Fluid of fused spheres as a model for protein solution
1 Svientsitskii St., 79011 Lviv, Ukraine \authorcopyrightM. Kastelic, Yu.V. Kalyuzhnyi, V. Vlachy, 2016
Aggregation of proteins in solution is both desired and undesired. It represents the first step in the downstream processing, i.e., salting out of the proteins for the purpose of cleaning and application. It is also one of the intermediate steps in the process of protein crystallization. The unwanted, pathological, protein aggregation is known to cause several diseases. Very importantly, bio-pharmaceutical drugs should be free of aggregates, otherwise they may be harmful. To increase the stability of protein in aqueous solutions is, therefore, an important technical problem. For an excellent review of the theoretical and experimental studies of protein solutions see reference .
The class of proteins we are interested in here are the so-called globular proteins. A typical representative of this class is lysozyme, which was extensively studied both experimentally and theoretically (see for example , Chapter 9). Despite their complicated 3D structure, many properties of protein solutions can be successfully described using relatively simple models [2, 3, 4, 5, 6, 7, 8]. Globular proteins are in their native form (we assume that during the experimental treatment they do not denature) most often pictured as perfectly spherical objects. This naïve representation is in reality never satisfied, it is used merely to simplify the calculations. There is a large list of non-spherical proteins, for example the shape of lysozyme mentioned above is ellipsoidal, including antibodies, lactoferrin, and others, where more complex geometry of particles would need to be taken into account to generate realistic results. This is important because the interactions leading to protein aggregation are directional and of short-range.
The shape of protein molecules influence their mutual interaction directly and indirectly. For example, (i) depletion interaction is largely dependent on the shape of particles . (ii) Experimentally, it is observed that many of proteins with roughly spherical shape exhibit upper critical solution temperature at protein concentration equal to 240 g/L [1, 10, 11]. In contrast to that, Y-shaped antibodies exhibit the shift toward much lower values, way down to 100 g/L [12, 13]. (iii) The hydrodynamic radii of the non-spherical objects are different, therefore, their hydrodynamic and transport coefficients , as well as, kinetic parameters  are modified. It is also known that classical nucleation theory has difficulties in describing the crystallization of other than spherical (for example ellipsoidal) particles [1, 16].
Recently, we used a simple spherical model  to analyze experimental results for the cloud-point temperatures in aqueous protein solutions with added salts [10, 11]. We modelled the solution as a one-component system; the protein molecules were represented as perfect spheres, embedded in the continuum solvent composed of water, buffer, and various simple salts. The attractive short-range interactions between the proteins were due to the square-well sites located on the surface of protein molecules. The model was numerically evaluated using Wertheim’s perturbation theory [17, 18, 19]. The short-range and directional nature of the interactions among proteins led to good agreement with the experimental data for the liquid-liquid phase diagram in case of lysozyme and -crystallin solutions [10, 11]. With knowledge of the experimental cloud-point temperature as a function of composition of electrolyte present in the system, the model gave predictions for the liquid-liquid coexistence curves, the second virial coefficients, and other similar properties under such experimental conditions.
One weakness of the model presented above was its simplified geometry. Neither lysozyme nor other proteins are spherical, and some of them for example, lactoferrin  look more like two fused spheres. The other weakness was that we assume for protein molecules to exist in form of monomers, which is not true. Even in very dilute solutions, proteins can be at least partially dimerized. The purpose of the present work is to investigate how the relaxing of these two basic assumptions influence the liquid-liquid coexistence curve.
The models for the association of spherically symmetric particles into dimer molecules are of considerable interest to science and technology and have been actively studied earlier. Of particular interest for us are the models where there is an inter-penetration (‘‘fusing’’ of cores) of the spherical particles upon association so that the bonding length is less than the core diameter . The ‘‘shielded sticky shell’’ and the ‘‘shielded sticky point’’ models of Stell and co-workers [21, 22, 23] and their extensions [24, 25, 26], belong to this group of models and are the starting point for our work. These types of the models were studied using regular [21, 22, 24] and multi-density [27, 28, 29, 23, 25, 26] integral equation theories.
In the present study, we use spherical particles as building blocks, which are fused together to form a new species. In this way, we compose the molecule with arbitrary spacing between the centers of spheres. Next, we decorate the surfaces of fused-sphere molecules with the attractive short-range binding sites, which may cause the association of the newly formed molecules. Such an extension of the protein model follows from our previous work . Here, we wish to explore the effects of the non-spherical shape on various thermodynamic properties.
Different versions of the model of dimerizing particles, represented by the tangentially bonded chain molecules, have been studied earlier [30, 31]. In this type of the model, dimerization occurs due to square-well bonding site, placed on the surface of one of the hard-sphere terminal monomer of each chain. Theoretical description of the model was carried out using first order thermodynamic perturbation theory (TPT1) of Wertheim [17, 18]. There are two major features of our model that set it apart from the models studied earlier, i.e., (i) in our model the molecules are represented by the two hard-sphere monomers fused at a distance , which is less than the contact distance and (ii) the molecules upon association can form a three-dimensional network. Due to the former feature of the model, a straightforward application of Wertheim’s multi-density approach fails to produce accurate results [29, 32, 33]. In the present work, we use a modified version of the TPT1, which takes into account the change of the overall packing fraction of the system due to the association forces [29, 34, 35]. The accuracy of our modified TPT1 approach is checked by the newly generated Monte Carlo simulation data.
2 Model, theory, and simulations
We introduce a one component model of spherical particles, decorated with additional binding sites of two different types, A and B. The binding site A is placed within the sphere, with the displacement , while an arbitrary number of binding sites of type B is located on the surface of the sphere (the displacement ). The model is visualized in figure 1. We consider a special case, where we exclude the cross interactions among sites A and B. The total pair potential is written as follows:
where is the pair potential for hard spheres, and and are inter-particle site-site potentials. The vector () connects the centers of hard spheres, and denotes the inter-particle vector connecting two sites of the type . As mentioned above, is the orientation dependent square-well potential between the sites , defined as follows:
The site A causes inter-penetration of particles (see figure 3). Note that we need the term to compensate for the hard sphere repulsion. For the periphery sites B, we do not need such a term, therefore . To fuse hard cores at separation , we choose and take the limit .
An appropriate theoretical approach to be used is the first-order Wertheim’s thermodynamic perturbation theory (TPT1) [17, 18]. According to this theory, the Helmholtz free energy of the system can be written as a sum of several terms:
Here, and is Boltzmann’s constant as usual, is the absolute temperature, and is the number of spheres. Further, defines the average number fraction of particles, which are not bonded through the binding site . Parameters and are determined by the mass-action law 
where is the number density of spheres and connects the pair distribution function of hard spheres (reference system) and the binding potential for sites and . The corresponding parameters are:
Expression for the solid-angle averaged Mayer function
was initially derived by Wertheim  and further generalized here to be
To suppress the cross interactions A–B, we set , which finally yields two independent equations, written in a quadratic form
Association parameters and
The association parameter is related to via equation (12) and to the free energy contribution due to A–A binding, by equation (5). For the complete association limit, i.e., fusing of hard cores at separation , no monomer spheres are present, so . We re-write the association parameter and introduce the cavity correlation function to obtain
where . Note that, as already mentioned before, . By applying the sticky limit approximation , that is by assuming the constant value of within the integration domain, we obtain
The integral given by equation (16) is not used in further calculations and, accordingly, it will not be considered in more detail here.
Parameter determines the degree of association of fused spheres and the free energy contribution due to the B–B binding, see equations (6) and (13). Notice that due to the association between -type of the sites, the packing fraction of fused spheres is different from the packing fraction originally present (un-fused) hard spheres . These fractions are related as follows:
The integral in can be evaluated analytically. We have used the Carnahan-Starling approximation for the contact value of at the effective packing fraction of fused spheres
Cavity correlation function
The last unknown quantity in equation (15) is the cavity correlation function of hard sphere fluid, . It is calculated by using the Tildesley-Streett expression for pressure of the hard dumbbell fluid . By choosing , and applying sticky limit conditions, i.e., , while keeping the integral in equation (14) finite, our model reduces to the hard dumbbell fluid. We modify the mass action law [equation (12)], by inserting equation (15) with , where stands for the number density of spheres, not bonded through binding site A (monomers). The result represents a different form of equation (110) of Wertheim’s paper 
Following Wertheim , we get the expression for the excess pressure in the form:
We are now in position to obtain the cavity correlation function of hard sphere system. We use the Carnahan-Starling equation of state  for the reference system () and the Tildesley-Streett equation of state  for the perturbed hard dumbbell system ().
where is the number density of hard dumbbells. The set of numerical parameters , , , , , is given in table 1.
Within the framework of Wertheim’s theory, we must set in equation (23) to recover the fluid of hard dumbbell particles (no monomers present). Next, we use equations (17), (18), (23) and equations of state [(24) and (25)] to obtain the derivative
The set of equations which determine and  are as follows:
with the arrays
It is obvious , therefore equation (26) can be easily integrated to yield:
The integral was checked to be non-singular for all investigated and values. Numerical results for are for a few fluid densities shown in figure 2.
Other thermodynamic properties
Next, we calculate the excess pressure and the excess chemical potential needed in phase diagram calculations as also excess internal energy , due to association. Starting with the pressure, we have
By inserting the appropriate derivatives from equations (4), (12), and (13), we get the final expression for the excess pressure. The second term is evaluated at , see equation (21), therefore upon differentiation we get an additional factor
The excess chemical potential is obtained through the relation
The second term in equation (37) is independent of density and, accordingly, does not contribute to the pressure. The expression for is the same as derived before [equations (32)–(34)]. The equilibrium conditions require the equality of chemical potential at a constant temperature (see equations below), so the second term in equation (37) cannot affect the coexistence curve. The equilibrium conditions read:
Another thermodynamic quantity is the excess internal energy , obtained as
Since is divergent, the only relevant part is . Derivative is obtained analytically from equations (19)–(20), since is independent. Thermodynamic functions for the reference system of hard spheres, , , and , can be found elsewhere .
2.3 Monte Carlo simulation
To validate the accuracy of the modified TPT1 approach, we performed Monte Carlo computer simulations in the ensemble . We assumed fused spheres with one and two binding sites on each sphere, where the prescribed arrangement of sites was preserved during the simulation. Simulated molecules are schematically shown in figure 3. We adopted the sampling method suggested by Tildesley and Streett , where a single displacement parameter was needed to describe the translation and rotation of fused spheres. The simulation box contained 250 fused spheres (molecules), which is equivalent to 500 penetrating (original, un-fused) spheres. We defined the cycle with 250 attempts to move the object and by 1 attempt to change the volume box. Next, we defined the block to be equal to cycles. Initially, we performed 1 block, to equilibrate the system, while 4 independent blocks were needed to calculate thermodynamic properties via the block averaging. Simulations were performed for three values: 0.2, 0.6, 1.0, and four different pressures : 0.5, 1.0, 2.0, and 4.0, for each model object visualized in figure 3. The acceptance rate of trial configurations was between 0.2 and 0.6.
3 Results and discussion
3.1 Thermodynamic properties: Theory against Monte Carlo simulations
To test the accuracy of TPT1 we performed Monte Carlo simulations for values of equal to 1 and 2 (three values for each ). We chose to compare the pressure and the internal energy due to B–B binding, . Since we used the complete A–A association limit within TPT1, the latter quantity, , was the one that could be directly compared to computer simulations. We fixed the temperature , while the pair potential characteristics are given in table 2.
The comparison between the theory and simulations is presented in figure 4. We found very good agreement for the pressure, while the theoretical predictions for were less accurate. In case of we obtained very good agreement for and fair agreement for (black lines and corresponding symbols in sub-figures). If we reduced the values (blue and red lines, symbols), the deviations became larger, though the qualitative picture remained correct. Deviations at low could be caused by the facts that: (i) fusing of two spheres at small is not a small perturbation regarding the reference system of hard spheres, and (ii) the arrangement of binding sites B is fixed during the simulation, which is not the case in TPT1, where the orientation average over all geometries was assumed.
3.2 Effect of protein’s shape on the liquid-liquid phase diagram
To illustrate the influence of protein shape on the liquid-liquid phase behavior, we compared the phase diagrams for two versions of the model: model (I) of two fused hard spheres and the model (II) of equivalent hard sphere, which was defined as a limiting case of (I), when and . The latter was chosen in such a way, that the volume of two fused spheres in case (I) was equal to that of the equivalent sphere (II): . Example (II) might be interpreted as the usual hard sphere model of diameter , with of sites B, i.e., the same number as on the two fused spheres. In such an interpretation, the A–A contributions to the physical properties can be neglected. Describing the aggregation of fused spheres of diameter within the limiting conditions and , where , led us to the model examined in reference . Other parameters and relations between examples (I) and (II) are listed in table 3.
|two fused spheres at (I)||equivalent sphere (II)|
|number of sites B:||(per building block)|
In figure 5, we show phase diagrams for the variants (I) and (II) described above, at three values and for different numbers of sites B. As observed before , an increase of the number of sites B shifts the critical density toward higher values. What is more interesting here is the effect of the separation distance parameter on the phase behavior. For a sufficiently small , i.e., — see figure 5 (a), the difference between the phase diagrams for two versions of the model becomes negligible, regardless of the value. If centers of fused spheres are located at larger distance, — see figure 5 (b), the difference becomes more pronounced: both critical temperature and density are lowered. Deviations become the strongest for the limiting example of two spheres fused in contact, that is for , cf. figure 5 (c). In this case, the larger number of sites (larger area available for interaction) additionally affects the liquid-liquid phase diagram. The shift toward lower critical densities (or packing fractions) is consistent with experimental studies of the Y-shaped antibodies [12, 13].
Proteins come in many shapes, from ellipsoidal to Y-like and are never perfectly spherical as treated by most theoretical models. Further, even in dilute solutions they have a tendency to form dimers and can be represented by two fused spheres. For dense systems close to precipitation, the actual geometry of the protein molecules is important; the inter-particle interactions are directional and of short-range. In the present study, we modify the first-order thermodynamic perturbation theory for associating fluids to be applicable to the models allowing hard-sphere particles to inter-penetrate. These particles can further aggregate. We confront theoretical predictions for thermodynamic properties of the proposed model with predictions of the corresponding Monte Carlo simulations. We obtain an excellent agreement for the pressure and fair agreement for the excess internal energy. Next, we use this model to predict the liquid-liquid phase diagram for protein solutions. We are interested in the effects of protein shape on the phase coexistence curve. We show that the fused hard-sphere model reduces the critical density of the system in comparison with the same quantity calculated for the hard-sphere model. This finding is consistent with experimental observations for Y-shaped antibodies. Using the cloud-point temperature measurements, we currently investigate the influence of various salts on the stability of lactoferrin solutions in water. The latter protein has a shape of two fused spheres, and the hard-sphere model is not a good representation of it. Theoretical approach developed in this paper will be used to analyze experimental data for lactoferrin and some other proteins of non-spherical geometry.
This study was supported by the Slovenian Research Agency fund (P1-0201), NIH research grant (GM063592), and by the Young Researchers Program (M.K.) of the Republic of Slovenia.
Плин зi спаяних сфер як модель розчину протеїнiв М. Кастелiч\refaddrlabel1, Ю.В. Калюжний\refaddrlabel2, В. Влахi\refaddrlabel1 \addresses \addrlabel1 Факультет хiмiї i хiмiчної технологiї, Унiверситет Любляни, вул. Вечна, 113, 1000 Любляна, Словенiя \addrlabel2 Iнститут фiзики конденсованих систем НАН України, вул. I. Свєнцiцького, 1, 79011 Львiв, Україна
- thanks: We dedicate this contribution to our friend and coworker Professor A.D.J. Haymet on occasion of his 60 birthday.
- Gunton J.D., Shiryayev A., Pagan D.L., Protein Condensation: Kinetic Pathways to Crystallization and Disease, Cambridge University Press, Cambridge, 2007.
- Abramo M.C., Caccamo C., Costa D., Pellicane G., Ruberto R., Wanderlingh U., J. Chem. Phys., 2012, 136, No. 3, 035103; \bibdoi10.1063/1.3677186.
Lomakin A., Asherie N., Benedek G.B., Proc. Natl. Acad. Sci. U.S.A.,
1999, 96, No. 17, 9465;
- Ruppert S., Sandler S.I., Lenhoff A.M., Biotechnol. Progr., 2001, 17, No. 1, 182; \bibdoi10.1021/bp0001314.
- Rosch T.W., Errington J.R., J. Phys. Chem. B, 2007, 111, No. 43, 12591; \bibdoi10.1021/jp075455q.
- Gogelein C., Nagele G., Tuinier R., Gibaud T., Stradner A., Schurtenberger P., J. Chem. Phys., 2008, 129, No. 8, 085102; \bibdoi10.1063/1.2951987.
- Carlsson F., Malmsten M., Linse P., J. Phys. Chem. B, 2001, 105, No. 48, 12189; \bibdoi10.1021/jp012235i.
- Kastelic M., Kalyuzhnyi Yu.V., Hribar-Lee B., Dill K.A., Vlachy V., Proc. Natl. Acad. Sci. U.S.A., 2015, 112, No. 21, 6766; \bibdoi10.1073/pnas.1507303112.
- Lekkerkerker H.N.W., Tuiner R., Colloids and the Depletion Interaction, Lecture Notes in Physics Vol. 833, Springer, Dordrecht, 2011; \doi10.1007/978-94-007-1223-2.
- Taratuta V.G., Holschbach A., Thurston G.M., Blankschtein D., Benedek G.B., J. Phys. Chem., 1990, 94, No. 5, 2140; \bibdoi10.1021/j100368a074.
- Broide M.L., Berland C.R., Pande J., Ogun O.O., Benedek G.B., Proc. Natl. Acad. Sci. U.S.A., 1991, 88, No. 13, 5660; \bibdoi10.1073/pnas.88.13.5660.
Mason B.D., Enk J.Z., Zhang L., Remmele R.L.J., Zhang J., Biophys. J.,
2010, 99, No. 11, 3792;
- Wang Y., Lomakin A., Latypov R.F., Laubach J.P., Hideshima T., Richardson P.G., Munshi N.C., Anderson K.C., Benedek G.B., J. Chem. Phys., 2013, 139, No. 12, 121904; \bibdoi10.1063/1.4811345.
- Yatsenko G., Schweitzer K.S., J. Chem. Phys., 2007, 126, No. 1, 014505; \bibdoi10.1063/1.2405354.
- Liu B.T., Hsu J.P., J. Chem. Phys., 1995, 103, No. 24, 10632; \bibdoi10.1063/1.469849.
- Wheeler M.J., Bertram A.K., Atmos. Chem. Phys., 2012, 12, No. 2, 1189; \bibdoi10.5194/acp-12-1189-2012.
- Wertheim M.S., J. Stat. Phys., 1986, 42, No. 3–4, 459; \bibdoi10.1007/BF01127721.
- Wertheim M.S., J. Stat. Phys., 1986, 42, No. 3–4, 477; \bibdoi10.1007/BF01127722.
- Chapman W.G., Jackson G., Gubbins K.E., Mol. Phys., 1988, 65, No. 5, 1057; \bibdoi10.1080/00268978800101601.
- Li W., Persson B.A., Morin M., Behrens M.A., Lund M., Oskolkova M.Z., J. Phys. Chem. B, 2014, 119, No. 2, 503; \bibdoi10.1021/jp512027j.
- Cummings P.T., Stell G., Mol. Phys., 1984, 51, No. 2, 253; \bibdoi10.1080/00268978700100861.
- Stell G., Zhou Y., J. Chem. Phys., 1989, 91, No. 6, 3618; \bibdoi10.1063/1.456894.
- Kalyuzhnyi Yu.V., Stell G., Llano-Restrepo M.L., Chapman W.G., Holovko M.F., J. Chem. Phys., 1994, 101, No. 9, 7939; \bibdoi10.1063/1.468221.
- Pizio O.A., J. Chem. Phys., 1994, 100, No. 1, 548; \bibdoi10.1063/1.466971.
Kalyuzhnyi Yu.V., Stell G., Holovko M.F., Chem. Phys. Lett., 1995,
235, No. 3–4, 355;
- Duda Yu.J., Kalyuzhnyi Yu.V., Holovko M.F., J. Chem. Phys., 1996, 104, No. 3, 1081; \bibdoi10.1063/1.470763.
- Wertheim M.S., J. Stat. Phys., 1984, 35, No. 1–2, 19; \bibdoi10.1007/BF01017362.
- Wertheim M.S., J. Stat. Phys., 1984, 35, No. 1–2, 35; \bibdoi10.1007/BF01017363.
- Kalyuzhnyi Yu.V., Stell G., Mol. Phys., 1993, 78, No. 5, 1247; \bibdoi10.1080/00268979300100821.
- Ghonazgi D., Chapman W.G., Mol. Phys., 1993, 80, No. 1, 161; \bibdoi10.1080/00268979400101881.
- Ghonazgi D., Chapman W.G., Mol. Phys., 1994, 83, No. 1, 145; \bibdoi10.1080/00268979400101141.
- Duda Yu., Lee L.L., Kalyuzhnyi Yu.V., Chapman W.G., Ting P.D., Chem. Phys. Lett., 2001, 339, No. 1–2, 89; \bibdoi10.1016/S0009-2614(01)00304-9.
Duda Yu., Lee L.L., Kalyuzhnyi Yu.V., Chapman W.G., Ting P.D., J. Chem.
Phys., 2001, 114, No. 19, 8484;
- Urbic T., Vlachy V., Kalyuzhnyi Yu.V., Dill K.A., J. Chem. Phys., 2007, 127, No. 17, 174511; \bibdoi10.1063/1.2784124.
- Kalyuzhnyi Yu.V., Hlushak S.P., Cummings P.T., J. Chem. Phys., 2012, 137, No. 24, 244910; \bibdoi10.1063/1.4773012.
- Hansen J.P., McDonald I.R., Theory of Simple Liquids, 3rd Edn., Elsevier, Boston, 2006.
- Wertheim M.S., J. Chem. Phys., 1986, 85, No. 5, 2929; \bibdoi10.1063/1.451002.
- Tildesley D.J., Streett W.B., Mol. Phys., 1980, 41, No. 1, 85; \bibdoi10.1080/00268978000102591.
- Frenkel D., Smit B., Understanding Molecular Simulations: From Algorithms to Applications, Academic Press, London, 1996.
- Streett W.B., Tildesley D.J., Proc. R. Soc. London, Ser. A, 1976, 348, No. 1655, 485; \bibdoi10.1098/rspa.1976.0051.
Bianchi E., Largo J., Tartaglia P., Zaccarelli E., Sciortino F.,
Phys. Rev. Lett., 2006, 97, No. 16, 168301;