[
Abstract
In this work we examine thermodynamics of fluid with ‘‘molecules’’ represented by two fused hard spheres, decorated by the attractive squarewell sites. Interactions between these sites are of shortrange and cause association between the fusedsphere particles. The model can be used to study the nonspherical (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 isobaricisothermal conditions. In particular, we are interested in the liquidliquid 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 nonspherical shape is to reduce the critical density and temperature. This finding is consistent with experimental observations for the antibodies of nonspherical shape. \keywordsnonspherical proteins, liquidliquid 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
Abstract
В ц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я збурень
201619223801
\doinumber10.5488/CMP.19.23801
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
1 Introduction
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, biopharmaceutical 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 [1].
The class of proteins we are interested in here are the socalled globular proteins. A typical representative of this class is lysozyme, which was extensively studied both experimentally and theoretically (see for example [1], 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 nonspherical 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 shortrange.
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 [9]. (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, Yshaped antibodies exhibit the shift toward much lower values, way down to 100 g/L [12, 13]. (iii) The hydrodynamic radii of the nonspherical objects are different, therefore, their hydrodynamic and transport coefficients [14], as well as, kinetic parameters [15] 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 [8] to analyze experimental results for the cloudpoint temperatures in aqueous protein solutions with added salts [10, 11]. We modelled the solution as a onecomponent system; the protein molecules were represented as perfect spheres, embedded in the continuum solvent composed of water, buffer, and various simple salts. The attractive shortrange interactions between the proteins were due to the squarewell sites located on the surface of protein molecules. The model was numerically evaluated using Wertheim’s perturbation theory [17, 18, 19]. The shortrange and directional nature of the interactions among proteins led to good agreement with the experimental data for the liquidliquid phase diagram in case of lysozyme and crystallin solutions [10, 11]. With knowledge of the experimental cloudpoint temperature as a function of composition of electrolyte present in the system, the model gave predictions for the liquidliquid 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 [20] 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 liquidliquid 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 interpenetration (‘‘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 coworkers [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 multidensity [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 fusedsphere molecules with the attractive shortrange binding sites, which may cause the association of the newly formed molecules. Such an extension of the protein model follows from our previous work [8]. Here, we wish to explore the effects of the nonspherical 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 squarewell bonding site, placed on the surface of one of the hardsphere 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 hardsphere monomers fused at a distance , which is less than the contact distance and (ii) the molecules upon association can form a threedimensional network. Due to the former feature of the model, a straightforward application of Wertheim’s multidensity 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
2.1 Model
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:
(1) 
where is the pair potential for hard spheres, and and are interparticle sitesite potentials. The vector () connects the centers of hard spheres, and denotes the interparticle vector connecting two sites of the type . As mentioned above, is the orientation dependent squarewell potential between the sites , defined as follows:
(2) 
The site A causes interpenetration 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 .
2.2 Theory
An appropriate theoretical approach to be used is the firstorder 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:
(3) 
where is the free energy of the reference system represented by the hardsphere system [36] and is the contribution due to A–A and B–B interactions. Following Chapman et al. [19], we have:
(4)  
(5)  
(6) 
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 massaction law [19]
(7)  
(8) 
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:
(9) 
Expression for the solidangle averaged Mayer function
(10) 
was initially derived by Wertheim [37] and further generalized here to be
(11) 
To suppress the cross interactions A–B, we set , which finally yields two independent equations, written in a quadratic form
(12)  
(13) 
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 rewrite the association parameter and introduce the cavity correlation function to obtain
(14) 
where . Note that, as already mentioned before, . By applying the sticky limit approximation [37], that is by assuming the constant value of within the integration domain, we obtain
(15)  
(16) 
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 (unfused) hard spheres . These fractions are related as follows:
(17)  
(18) 
where is the packing fraction of hard spheres and is reduced A–A bonding distance. Using the sticky limit approximation [37] for [equation (9)], we have:
(19)  
(20) 
The integral in can be evaluated analytically. We have used the CarnahanStarling approximation for the contact value of at the effective packing fraction of fused spheres
(21) 
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 TildesleyStreett expression for pressure of the hard dumbbell fluid [38]. 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 [37]
(22) 
Following Wertheim [37], we get the expression for the excess pressure in the form:
(23) 
We are now in position to obtain the cavity correlation function of hard sphere system. We use the CarnahanStarling equation of state [36] for the reference system () and the TildesleyStreett equation of state [38] for the perturbed hard dumbbell system ().

CarnahanStarling EOS:
(24) 
TildesleyStreett EOS:
(25)
where is the number density of hard dumbbells. The set of numerical parameters , , , , , is given in table 1.
0.37836  1.07860  1.30376  1.80010  2.39803  0.35700 
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
(26) 
The set of equations which determine and [] are as follows:
(27)  
(28)  
(29) 
with the arrays
It is obvious , therefore equation (26) can be easily integrated to yield:
(30) 
The integral was checked to be nonsingular 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
(31) 
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
(32)  
(33)  
(34) 
The expression is obtained from equation (26), while the second derivative is obtained analytically at from equation (21)
(35) 
The excess chemical potential is obtained through the relation
(36) 
The logarithmic term in equation (5) is divergent for the complete association limit (), therefore we rewrite this term by using equation (15) as follows:
(37) 
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:
(38)  
(39) 
where and are the two coexisting densities. At this step, the phase diagram can be constructed by applying equations (38)–(39) as it is in more detail explained in the previous work [8].
Another thermodynamic quantity is the excess internal energy , obtained as
(40)  
(41)  
(42) 
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 [36].
2.3 Monte Carlo simulation
To validate the accuracy of the modified TPT1 approach, we performed Monte Carlo computer simulations in the ensemble [39]. 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 [40], 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, unfused) 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.
:  0.1 
:  5.0 
:  0.5 
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 subfigures). 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 liquidliquid phase diagram
To illustrate the influence of protein shape on the liquidliquid 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 [8]. Other parameters and relations between examples (I) and (II) are listed in table 3.
two fused spheres at (I)  equivalent sphere (II)  

diameter:  
:  
:  
:  0.1  0.1 
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 [41], 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 liquidliquid phase diagram. The shift toward lower critical densities (or packing fractions) is consistent with experimental studies of the Yshaped antibodies [12, 13].
4 Conclusions
Proteins come in many shapes, from ellipsoidal to Ylike 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 interparticle interactions are directional and of shortrange. In the present study, we modify the firstorder thermodynamic perturbation theory for associating fluids to be applicable to the models allowing hardsphere particles to interpenetrate. 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 liquidliquid phase diagram for protein solutions. We are interested in the effects of protein shape on the phase coexistence curve. We show that the fused hardsphere model reduces the critical density of the system in comparison with the same quantity calculated for the hardsphere model. This finding is consistent with experimental observations for Yshaped antibodies. Using the cloudpoint 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 hardsphere 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 nonspherical geometry.
Acknowledgements
This study was supported by the Slovenian Research Agency fund (P10201), 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в, Україна
Footnotes
 thanks: We dedicate this contribution to our friend and coworker Professor A.D.J. Haymet on occasion of his 60 birthday.
References
 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;
\bibdoi10.1073/pnas.96.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., HribarLee 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/9789400712232.
 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;
\bibdoi10.1016/j.bpj.2010.10.040.  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/acp1211892012.
 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., LlanoRestrepo 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;
\bibdoi10.1016/00092614(95)00106E.  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/S00092614(01)003049.

Duda Yu., Lee L.L., Kalyuzhnyi Yu.V., Chapman W.G., Ting P.D., J. Chem.
Phys., 2001, 114, No. 19, 8484;
\bibdoi10.1063/1.1363667.  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;
\bibdoi10.1103/PhysRevLett.97.168301.