Soft Matter

Please note that technical editing may introduce minor changes to the text and/or graphics, which may alter content. The journal’s standard Terms & Conditions and the Ethical guidelines still apply. In no event shall the Royal Society of Chemistry be held responsible for any errors or omissions in this Accepted Manuscript or any consequences arising from the use of any information it contains. Accepted Manuscript


Introduction
Much of biology depends on proteins interacting with each other -pairwise or in form of oligomers -mediated by solvent and different ligands, such as salts and excipients. 1,2We need to understand how protein molecules interact with water, ions and inbetween them.The self-assembly of proteins into various structures plays a crucial role in biology.A better understanding of these processes is expected to be reached by modern computer simulations, which have become indispensable for insights into physics and chemistry.Yet, the systems of several protein species in water with salts and other co-solvents, which are relevant to biology and pharmacy, are often too big to be handled well enough by current explicit-water molecular simulations.A somewhat different, less detailed-, approach, based on the ideas of condensed matter physics has been forwarded in last decades.Interactions of mixtures of oppositely charged proteins have been studied experimentally 3 and theoretically 4 .Important previous studies rele-vant to our work are listed as refs.5 and 6.Two excellent reviews of developments in this area of research have been published recently. 7,80][11] First, the self-assembly is a step in protein crystallization; good quality crystals are needed to determine protein structure through diffraction experiments.The latter pathway is the following: soluble protein aggregates initialize the coexistence of low-and high density liquid phases, differing in protein concentration.The so-called liquid-liquid phase separation is a meta-stable stage and may further lead to crystallization.For this reason, many experimental and theoretical efforts, were invested to improve understanding of the liquid-liquid phase separation.Pioneering work of George and Wilson 12 , and later study of Vliegenthart and Lekkerkerker 13 revealed strong correlation of liquid-liquid phase boundary with the second virial coefficient and thus inter-protein interactions.Secondly, protein self-assembling in living cells plays a key role in diseases, such as Alzheimer's, Parkinson's, and others. 1,14,15hirdly, a key step in developing biotech drugs 16 , is to formulate proteins that do not form aggregates.The challenge is to deliver formulations with high concentrations of solutes, low vis-cosity and maintaining the protein therapeutics properties for at least two years.Self-assembly decreases therapeutic efficacy and increases its viscosity.
As atomistic level molecular simulations are not practical for studying phase separations in such solutions, the traditional approach is to use the colloidal model of Deryaguin-Landau-Verwey-Overbeek (DLVO). 17In this approach, proteins are treated as spheres that interact via van der Waals and screened Coulomb potentials, acting between centers of molecules.The approach proved to be unsatisfactory in several respects. 19][20][21] This strategy was successfully used to model a wide class of globular proteins and to analyze the liquid-liquid phase separation 22,23 , crystallization 24,25 , osmotic pressures 26 , distribution of aggregates 27,28 , percolation threshold 27 , and other physico-chemical properties. 21eras and co-workers 29,30 extended these ideas to binary mixtures of equally-sized spherical particles.They examined an influence of pair potential characteristics on the phase behavior under isobaric conditions.A special interest for us is the work of Roldnán-Vargas 31 and coworkers.These authors treated binary mixtures of spherical particles of different sizes and with different surface patterns.The fact that the shape of protein molecules affects the phase separation, has recently been confirmed by us.A distinct shift of liquid-liquid critical density toward lower values has been obtained in case of the associating dumbbells. 32Such a behavior was before observed experimentally for the system of Y-shaped antibodies. 33For proteins in dense mixtures near the phase separation a realistic modeling of protein shape seems to be important.
In the previous studies 34,35 we applied simple models with anisotropic interactions to analyze aqueous protein solutions in presence of low-molecular-mass salts.We treated proteins as hard spheres, decorated with the square-well-energy binding sites, using Wertheim's thermodynamic perturbation theory. 36he model provided good fits to the cloud-point temperature curves of lysozyme in buffer-salt mixtures as a function of the type and concentration of salt.The necessary condition required for such modeling to be realistic is that proteins in solution during the experiment remain in their compact form. 37,38Unfortunately, many experimentalists often do not provide this information.In the next contribution 32 the approach above was extended to examine the thermodynamics of fluid with "molecules" represented by two fused hard spheres, again decorated by the attractive square-well sites.Interaction between these sites was of shortrange and caused association between the fused-sphere particles.The model served to estimate the effects of non-spherical geometry of proteins on solution properties.
Here we continue with applications of simple models to protein systems.Going one step further, we treat a binary mixture of proteins with molecules of each of the respective species represented as dumbbells.In the first step we need to form the molecule (dumbbell) by fusing two hard spheres and further decorate it with attractive sites on the surface.These sites allow association between the same type of molecules (self-association) as also with molecules of the other protein species present in the system.Thermodynamic properties of the model mixture are obtained via Wertheim's thermodynamic perturbation theory. 36e utilized the McMillan-Mayer approach, where the principal variables are concentration, chemical potential of composed solvent, and temperature.We calculated the liquid-liquid phase diagram and present the results as the volume packing fractioncomposition -temperature diagram.Under such conditions the osmotic pressure follows from the equation of state.When it comes to comparison with experimental data, we chose to model γD-β B1 crystallin mixture, thoroughly studied by Wang and coworkers. 39Both proteins show significant deviations from spherical shape and resemble dumbbells 40 (see Fig. 1 therein, and Fig. 6 of this paper).These proteins are of special interest: an improper functioning of γD crystallin may lead to cataract in human eyes. 14For this reason crystallins were intensively studied.The investigations include stability measurements 41 , influence of mutation studies 42,43 and other surface modifications 44 , and the analysis of surface residues forming inter-protein salt bridges. 457][48][49] Investigation of this subtle phenomenon is beyond the scope of this work.In this initial study we focus to the liquid-liquid phase equilibrium, and how the latter is governed by the interactions between two protein species.We use a simple "dumbbell" protein model, coupled with the thermodynamic perturbation theory of Wertheim 36 , to analyze the experimental data of Wang and co-workers. 39

Protein molecules as dumbbells decorated by attractive sites
In the first step of the calculation we compose the dumbbell-like molecules of type i with number density ρ hard-spheres through the sticky sites B (see Fig. 1).
The sticky BB interaction, u BB (z BB ), is written as:

Soft Matter Accepted Manuscript
where r ij denotes the distance between centers of spheres of the type i and j and Ω i , Ω j their orientations.Furthermore, z BB is the distance between sites B, . . .Ω 1 Ω 2 denote orientation average, δ (. ..) is the Dirac delta function, and δ ij is the Kroneker delta symbol.The dumbbells are formed when K (11)   BB and K (22)   BB → ∞.For simplicity we assume that the hard spheres forming the protein species are of equal size σ .In the experimental case analyzed here this is not far from being true.The number of the attractive sites of type A α , located on a sphere belonging to the molecule i is m i .Further, α runs over 1, . . ., m i , while index i is assuming values 1 and 2. The number density of the system is ρ Protein molecules i and j are interacting via the pair potential, which consists of the sum of four terms u ij (r ij , Ω i , Ω j ).Note again that hard-spheres labeled 1 form molecules 1, while the hardspheres denoted 2 form molecules of type 2. These terms describe the interaction between decorated hard spheres i and j, belonging to molecules i and j: Here u hs (r ij ) is the hard-sphere potential, z αβ is the distance between sites A α and A β , which belong to the spheres of the type i and j, ε ij and ω ij are the square-well depth and width, respectively.In summary, we have four types of, i.e. 1-1, 1-2, 2-1, and 2-2 interaction between the decorated hard spheres denoted 1 and 2.

Thermodynamic perturbation theory
Numerical results for the model described above are obtained by Wertheim's thermodynamic perturbation theory (TPT). 36The basic quantity, Helmholtz free energy A, is composed of the ideal (id), hard-sphere (hs), and association term (ass): These terms are written as: where β = 1/k B T and T is the absolute temperature, Λ i is the de-Broglie thermal wavelength 50 , hs = (1 + η/2)/(1 − η) 2 is the Percus-Yevick expression for the contact value of the hard-sphere radial distribution function 51 , and X i is the fraction of the particles not bonded at any site A α of the hard sphere of type i.In addition we assume two types of steric incompatibilities: (i) the site A α on one molecule cannot bond simultaneously to two A α sites on another molecule; and (ii) double bonding between molecules is not allowed; see also Fig. 1 and the ref.52 (Fig. 2 therein).X i follows from the statistical-mechanical analogue of the mass action law 53 : where Here fij (r) is the orientational average of the Mayer function for the square-well site-site interaction: Wertheim's TPT uses the orientation average of f ij (r) and accordingly placement of the sites does not enter the theory.There is computer evidence in literature that this is not a severe approximation. 54ce the free energy A is known, the chemical potential of the protein species i, µ (p) i , and the equation of state (pressure P) can be obtained from the standard thermodynamic relations 50

Conditions of phase equilibrium
Binary system with protein densities {ρ (p) 2 (0)} may undergo phase decomposition into Q phases, denoted here as (1), (2) . . .(γ) . . .(Q).According to the Gibbs phase rule Q ≤ 4. The notation {ρ (p) 2 (0)} denotes the (unstable) "mother" phase, while the notation {ρ (p) 2 (γ)} applies to the (stable) "daughter" phases (γ).Thermodynamic conditions to be satisfied in equilibrium read: In words: chemical potential of species i = 1, 2 within phase (γ), µ (p) i (γ), must be at given temperature and pressure equal in all phases.In addition to these conditions, we need also to consider the conservation of the volume fractions of phases X(γ) = V (γ)/V (0) and the number densities of proteins {ρ (p) Eqs. (13-16) form the set of expressions for unknowns {ρ 2 (0)} of the mother phase.In order to determine the coexistence region we need to numerically solve the set of equations defined above.This task may be accompanied by several problems.First, the solution of the system of nonlinear equations is not unique and the exact number of solutions is not known in advance.Second, numerical methods for solving nonlinear equations are usually based on iterative algorithms, requiring good initial approximation to reach the convergence.Third, in a process of searching for the solution one can get stranded at a saddle point, away from the global free energy minimum.In general, additional restrictions 55 must be considered, what makes such an approach time consuming.
To avoid these complications, we used a somewhat different approach, searching for the global minimum directly.The appropriate functional to be minimized is the free energy of the system per volume, A(0)/V (0), To account for the conservation conditions given by Eqs. ( 15) and ( 16), we insert these expressions into Eq.( 17), and thus reduce the dimensionality of the problem.Next we use the Boltzmann simplex simulated annealing method (BSSA), described in Chapters 10.4 and 10.9 of ref. 56, to find the global minimum.When this minimum is reached, the conditions given by Eqs.(13-14), are automatically fulfilled.Numerically, the equations are satisfied within 10 −2 %.More detailed description of the search for the global minimum is given as a part of the Supporting Information.

Model parameters shape the phase diagram
The initial case to be examined is a mixture where both proteins have four sites per sphere, m 1 = m 2 = 4, equal square well depths ε 11 = ε 22 = 13.7ε, and corresponding ranges ω 11 = ω 22 = 0.061σ .In other words, we are mixing the same type of molecules.Note that in this chapter (and here only), we use the reduced units: energy is measured in ε, temperature as T * = k B T /ε, and length in diameter of sphere σ .Results for such a "symmetric" mixture are shown on panel A of Fig. 2. The states for temperatures below the equilibrium surface are unstable and decompose into lowand high-(volume) packing fraction (η) liquid phases.The "mix-ture" studied above represents the simplest possible example; it translates the symmetry of the pair potentials into topology of the phase diagram.The fact that two protein species have equal interaction parameters and sizes causes for the curves belonging different x 2 values to be parallel to each other.When the phase separation occurs, the decomposition into low-and high-packing fractions is purely density driven.This is shown on the right site of panel A (two-and one-phase regions are labeled by numbers "2" and "1", respectively).In other words, the lines connecting coexisting phases, the so-called tie lines are in this case parallel to the η axis.
To study the deviations from the basic case examined above, we place unequal numbers of sites on each protein (m 1 = m 2 ), while keeping other parameters unchanged.We examine the case where m 1 is increased by one (m 1 = 5), and m 2 decreased by one (m 2 = 3).The corresponding phase diagram is shown on panel B of Fig. 2. In contrast to the previous example, the x 2 curves are not parallel to each other and a distinct asymmetry along x 2 axis is observed.Envelopes for pure components, denoted by black solid lines at x 2 = 0 and x 2 = 1, are different: protein 1 has higher critical density and temperature, since it has five attractive sites (per sphere) than the protein species labeled 2. This observation is consistent with our previous study 34 , as well as, with observations of other authors studying one-component systems. 57,58n addition to this, the unstable states exhibit higher complexity: as inferred from the right site of panel B, T * = 1.20, the high-packing fraction phase tends to retain the composition of the unstable state (density driven), while the low-packing fraction phase is enriched on an account of the second component (composition driven), to satisfy the conservation of particles.At higher temperature, T * = 1.40, this effect is less pronounced.
In the next example the differences between protein species are introduced on account of modifying their interaction energies ε ij .It is important to stress out that overall protein-protein interaction is controlled by ε ij and ω ij , as well as, by the number of binding sites m i .Here the square well depth ε 11 is increased for 10 %, while keeping the number of sites m 1 = 5 and m 2 = 3, and ranges ω ij unchanged.These results are shown on panel C of Fig. 2. On the first sight there seems to be little difference to panel B, except for higher critical point of the envelope at x 2 = 0.A closer look, however, reveals the coexistence region of three liquid phases (black shaded area), embedded into the two-phase region.Interestingly, the three-phase region disappears for the temperature equal to the critical value of pure protein 2 (x 2 = 1) that is at T * = 1.225.To demonstrate a simultaneous presence of two-and three-phase regions, we on the same panel show also the isotherms for temperatures below and above T * = 1.225.The isotherm for T * = 1.20 reflects the presence of three-phase region, labeled by number "3", and surrounded by closed dashedblue loop.Unstable states within this region (denoted by empty squares) decompose into the phases with packing fraction and composition as indicated on the figure.We may call these three phases as low-, intermediate-, and high-packing fraction liquid phases.According to the Gibbs phase rule, in a two-component mixture the region where three phases are in equilibrium has only one degree of freedom, in our case this is temperature.The three- Finally, in contrast with the previous example, we keep ε 11 = 13.7ε(ε = k B T /T * ) unchanged and increase the ε 22 for 10 %.As before m 1 = 5 and m 2 = 3.In this way we suppress -in comparison with the previous case -the interaction potential difference between the two protein species.The resulting phase diagram is shown on panel D of Fig. 2. In this case the phase behavior exhibits higher symmetry along x 2 axis than in panels B and C: the shape of the envelope at x 2 = 1.0 resemble the one at x 2 = 0.0.Interestingly, the three-phase region shown in panel C disappears under such conditions.Because of similarity in interactions among the two protein species, the system does not need a third phase for stabilization.Furthermore, the high packing fraction (daughter) phases, derived from unstable states (mother phases), tend to retain the composition of the mother phases, similarly as before shown on panel B. The best-fit parameters, collected in Table 1, were determined by hand; initially we pick m 1 , ε 11 , and ω 11 to fit the envelope of pure γD protein (x 2 = 0.0).Experimental results for pure β B1 (x 2 = 1.0) are not available and to adjust m 2 , ε 22 , and ω 22 values, we fit the isotherm x 2 vs η at T = 268.5K (see Fig. 3).In this range of x 2 values the T cloud data are readily available. 39This is all the input information we use.We also assume that conventional Lorentz-Berthelot mixing rules provide a reasonable measure of the cross interactions; it is known that an excessive deviations in parameters may cause an un-physical instability of protein mixtures. 59arameters listed in Table 1 are then used to calculate (i) the isotherm at T = 270.5K (including the tie-lines), for which the T cloud measurements are available (see Fig. 3), (ii) T cloud vs η curves for x 2 equal to 0.08, 0.16, and 0.29, shown in Fig. 4 as colored curves, (iii) the isotherms at 255 K, 265 K, and 275 K presented in Fig. 5, as also (iv) the equilibrium surface for the model γD-β B1 mixture (see Fig. 4) for a range of temperatures, x 2 , and η values.Circles show experimental data, while the solid lines represent calculations.In addition, we compare the experimentally determined tie-lines 39 (dashed green/orange lines), with those calculated (black solid lines).Experimental data for the isotherm at 268.5 K (without the tie-lines) were together with the measurements for x 2 = 0 (Fig. 4) used to fix the model parameters, collected in Table 1.Other results, including the isotherm at 270.5 K and the tie-lines for both temperatures, were than calculated using these values.

Analysis of experimental data for γD-β B1 mixtures
Results for the calculated tie-lines and isotherm at 268.5 K are, as shown in Fig. 3, in fair agreement with measurements.fact that we obtained correct slopes for the tie-lines supports the proposed model; it is known that the latter are very sensitive on the protein-protein interactions.As concluded by Dorsaz 6 et al; the difference in cross interactions for less then k B T may already invert the slope of the tie-line.Notice that uncertainties in measurements of the packing fraction and composition of the coexisting phases can be up to ±0.02. 39In Fig. 4 we show the agreement between the T cloud measurements and theoretical calculations (lines) for x 2 values equal to 0.00, 0.08, 0.16, and 0.29.In addition to these results we also show the equilibrium surface, which signals the onset of cloudiness and the coexistence of low-and high-packing fraction (η) liquid phases.Large difference between m 1 (equal to eight), and m 2 (equal to two) for γD and β B1, we picture this situation in Fig. 1, is needed for the model to correctly describe experimental results.This set of the m i values leads to the presence of three-phase region (black shaded area) at temperatures below 264 K, where no experimental results are available.Notice that the two-phase coexistence region is present even for packing fractions as high as η > 0.4, what was not the case when the difference between m 1 and m 2 was smaller (compare Fig. 2).β B1 protein exhibits a slightly longer range ω 22 as also a somewhat stronger association energy ε 22 .The measurements of the association enthalpies 60 , based on the thermodynamic analysis of protein oligomerization, reveal much larger association enthalpy in

Soft Matter Accepted Manuscript
case of β B1 crystallin.This finding is consistent with our calculations.It is worth mentioning that the potential parameters (ε 11 , ε 22 and ω 11 , ω 22 ) applied in this calculation fall in a broad range of hydrogen bond magnitudes 61 .This result is consistent with our previous findings. 34 4. We show the η-x 2 projections of the equilibrium surface for three different temperatures: 255 K, 265 K, and 275 K.The graphs show an evolution of the phase coexistence curves with temperature.For the highest one, T = 275 K, the one-phase region is dominant; the two-phase regions only exist around x 2 =0 (η between 0.1 and 0.25) as also for high η values with the composition x 2 between roughly 0.4 and 0.8.For somewhat lower temperature, T = 265 K, the twophase region is much broader.The lines indicate how the unstable points (squares) phase separate into the phases having different x 2 , η values.For the lowest isotherm, T = 255 K, the graph is quite complex.Here we have a domain, denoted by "3", where three liquid phases coexist.In addition both one-and two-phase regions are present.The complexity is very likely due to fact that interprotein interactions gain on their importance if temperature is lowered.See also the discussion of results presented on panel C of Fig. 4.
Although the choice of the set of parameters is not unique, there seems to be a limited number of such sets providing reasonable results.In several examples we tested the sensitivity of the parameters listed in Table 1.In one of such cases we assumed for the number of sites to be equal m 1 = m 2 = 8, and ranges ω 11 = ω 22 = 0.18 nm, while at the same time, we increased the difference between ε 11 and ε 22 for 12 %.With these modifications, we have not been able to reproduce the experimental isotherm at T = 268.5 K, shown in Fig. 3: the result for x 2 ≈ 0.3 at η ≈ 0.14 was in our best fitting attempt too large by a factor of two.In yet another example we choose m 1 = 8 ("fixed" by fitting the ex-periment at x 2 = 0) and m 2 = 4, keeping the same temperature and ω 11 and ω 22 values as above.Again, even upon considerable adjustment of ε 22 (we took for 15 % more negative value) the agreement with experimental x 2 vs η curve was poor, much poorer than in Fig. 3.
Authors of the experimental work provided data in the domain η × x 2 ∈ [0, 0.3] × [0, 0.3]; the state points out of this range were not accessible to measurements due to the tendency of the system to solidify.We calculated the system properties for η and x 2 values out of this range -no convergence problems were encountered.An example is the liquid-liquid envelope at x 2 = 1, i.e. for pure β B1 protein.Unfortunately, this curve could not be obtained by direct measurements, because of the low critical temperature below −10 o C. 60 Addition of polyethylene glycol to the protein buffer mixture increases the critical point and allows to measure the cloud point temperatures in a wide range of protein concentrations.Such an indirect measurements 60 , in conjunction with extrapolation to zero concentration of polyethylene glycol, suggest for the critical point to be around 250 K. Our theoretical model yields for x 2 = 1 (see Fig. 4) the value equal to 255 K. Notice again that this prediction is based on fitting the results for x 2 = 0 (absence of β B1 protein) and the isotherm for the mixture at T = 268.5 K.
An alternative (and to a certain degree complementary) presentation of the information given by phase diagram in Fig. 4 and isotherms in Fig. 5 is provided by the cloud and shadow curves. 62These curves are fundamental for studying mixtures, because they show the partitioning of protein species into lowand high-packing fraction phases.These results are shown in SI as Fig. (S1).

Relating model parameters to structural information
Crystallographic data for both proteins are available from the Protein Data Bank (PDB) server.While for γD (1hk0) complete structure is provided, for β B1 only a truncated structure is available.The truncation applies to the hydrophobic tail (41 amino acids) at the N-terminus, named β B1∆N41 (1oki).Even though the absence of hydrophobic tail may modify phase behavior 60 , the quasi-elastic light scattering and circular dichroism measurements of the wild-type β B1 and mutant β B1∆N41 revealed no significant differences. 60e surface distributions of charges and insights into tertiary structures of γD and β B1∆N41 proteins are shown in Fig. 6.Since all the measurements were performed slightly below isoelectric points of proteins (7.1), at pH equal to 7.0, the net charge of proteins is around zero (blue and red surface areas are roughly of the same size).From the Fig. 6 it is difficult to draw any firm conclusions about why the numbers of sites per sphere (m 1 , m 2 ), needed to adequately describe crystalline mixture, differ from one protein to another so much.In search for some clues we turned to the Protein Dipole Moments Server 66 (http://dipole.weizmann.ac.il/) for information about higher multipole terms for 1hk0 and 1oki.The results are summarized in Table 2.The two proteins are having dipole moments of similar magnitude, while β B1∆N41 has more than nine Soft Matter times larger quadrupole moment than γD.This enhances the  The comparison of the model calculations with experimental data for binary mixture of γD and β B1 crystallins yields the following conclusions: (i) for good agreement between theory and experiment a large difference in number of sites per sphere of protein (m 1 = 8 and m 2 = 2) is needed.The results are consistent with experimental findings indicating that β crystallins form small clusters with only few proteins involved.We trace this result to the large difference in the quadrupole moments between the two proteins.(ii) Our model predicts correct slopes for the tie-lines in x 2 vs η graphs.(iii) The calculation yields good agreement with experimentally determined temperature-packing fraction curves for mixtures.(iv) Model provides the results for pure protein β B1, for which no direct measurements are available, and accordingly no these data were included into fitting procedure.Calculation for the critical temperature is in good agreement with the value suggested from the extrapolation of experimental data for β B1 protein.(v) The "dumbbell" model predicts a region where three liquid phase are in equilibrium -the region is just out of the domain where measurements are currently available.The calculation demonstrates, in agreement with articles listed in the recent review papers 7,8 , that methods of the condensed matter physics may be useful in analysis of experimental data of protein systems.

Fig. 1
Fig. 1 Hard spheres of equal size form dumbbells (model for protein molecule) through the sticky site B. Further interaction between protein molecules is possible through A α sites.In this scheme m 1 = 8 and m 2 = 2.
r n a l N a me , [ y e a r ] , [ v o l .] ,
Encouraged these results we turned to analyze experimental data for γD-β B1 protein mixtures, published in ref. 39.The au-thors measured the cloud point temperatures, T cloud , as a function of composition of the mixture.The crystalline proteins γD and β B1 are approximately of the same size, what makes our analysis somewhat simpler.Experimental data are converted from the volume packing fraction values φ 39 to corresponding η units via the relation: η= πφ N A σ 3 /[3Ω(M 1 + x 2 (M 2 − M 1 ))],where N A is Avogadro number, M 1 = 20, 607 g mol −1 and M 2 = 27, 892 g mol −1 molar masses of proteins, and Ω = 7.1 • 10 −4 L g −1 the specific volume of proteins.39

Fig. 3
Fig. 3 Isotherms at 268.5 K (green-bottom), and 270.5 K (orange-top).Circles show experimental data, while the solid lines represent calculations.In addition, we compare the experimentally determined tie-lines39 (dashed green/orange lines), with those calculated (black solid lines).Experimental data for the isotherm at 268.5 K (without the tie-lines) were together with the measurements for x 2 = 0 (Fig.4) used to fix the model parameters, collected in Table1.Other results, including the isotherm at 270.5 K and the tie-lines for both temperatures, were than calculated using these values.

Fig. 4
Fig. 4 Panel A: phase separation of the γD-β B1 mixture in 0.1 M phosphate buffer at pH = 7.0.We show the equilibrium surface and its projections on the η-x 2 plane.Experimental results are denoted by points and calculations by lines.The comparison is shown for x 2 values equal to 0, 0.08, 0.16, and 0.29.The three-phase region (black shaded area) is embedded into two-phase region for temperatures below 264 K. Panel B: for easier visualization we re-plotted the three-phase region.Coexisting phases are denoted by black lines, see also the isotherm at T = 255 K in Fig. 5.
r n a l N a me , [ y e a r ] , [ v o l .] ,

Fig. 5
Fig. 5 Isotherms for the temperatures 255 K, 265 K, and 275 K, extracted from Fig. 4. The color notation corresponds to the temperature code in Fig 4. The one-, two-, and three-phase regions are denoted by the appropriate numbers.Squares denote mixtures, which phase separate while the lines are connecting the mother phases with the coexisting (daughter) phases.

Fig. 5
Fig.5presents the results for the three isotherms, which are part of the full phase-diagram shown in Fig.4.We show the η-x 2 projections of the equilibrium surface for three different temperatures: 255 K, 265 K, and 275 K.The graphs show an evolution of the phase coexistence curves with temperature.For the highest one, T = 275 K, the one-phase region is dominant; the two-phase regions only exist around x 2 =0 (η between 0.1 and 0.25) as also for high η values with the composition x 2 between roughly 0.4 and 0.8.For somewhat lower temperature, T = 265 K, the twophase region is much broader.The lines indicate how the unstable points (squares) phase separate into the phases having different x 2 , η values.For the lowest isotherm, T = 255 K, the graph is quite complex.Here we have a domain, denoted by "3", where three liquid phases coexist.In addition both one-and two-phase regions are present.The complexity is very likely due to fact that interprotein interactions gain on their importance if temperature is lowered.See also the discussion of results presented on panel C of Fig.4.Although the choice of the set of parameters is not unique, there seems to be a limited number of such sets providing reasonable results.In several examples we tested the sensitivity of the parameters listed in Table1.In one of such cases we assumed for the number of sites to be equal m 1 = m 2 = 8, and ranges

Fig. 6
Fig.6 Comparison between γD (left) and β B1∆N41 (right) crystallins.Panel a: tertiary structures embedded into solvent accessible area, generated by Schrödinger software63 .Panel b: charged amino acid residues on surfaces of proteins; positive -blue, negative -red, neutral -white.The location of residues is shown in the Hammer projection64 , having the polar and the azimuthal angles as the principle coordinate axes.An implementation of this projection was proposed by Koromyslova65 et al and modified by us, as described in SI.
We use a simple model to calculate the phase diagram for the mixture of two proteins.Both protein species are modeled as dumbbells (two fussed hard spheres) decorated by different numbers of attractive square-well sites allowing inter-protein attraction.Major parameters of the model are: number of attractive sites per sphere forming the protein (m i ; i=1,2), their attractive range and square-well depth.In the fist part of the study we investigated how the principal model parameters shape the phase diagram.This initial calculations showed that number of the attractive sites per sphere of protein plays an important role.In the second part of the study we turned to analyze the behavior of the lens proteins.The mixture of γD and β B1 crystallins was thoroughly studied experimentally by Wang39 et al, and the authors provided the cloud point temperatures for various experimental conditions.To fix the principal parameters of our model -in particular m 1 and m 2 and the attractive square-well depths -we fitted experimen-tal results for pure γD protein and the isotherm for the mixture at T = 268.5 K.The cross interaction between proteins was described by the usual Lorentz-Berthelot rule.Numerical evaluation of this model was performed by Wertheim's thermodynamic perturbation theory, which allows calculation of the measurable properties, including the complete phase diagram.
11 , ω 22 .For the cross interactions we use standard Lorentz-Berthelot mixing rules: (ε 12 ) 2 = ε 11 ε 22 and 2ω 12 = ω 11 +ω 22 .Initially we perform a systematic analysis of how main model parameters affect the phase diagrams, presented as volume packing fraction (η) -composition (x 2 = ρ Our goal here is to study stability of mixtures of two protein species.The main parameters of the model are: (i) number of binding sites m 1 and m 2 per sphere (for protein species 1 and 2); (ii) the attractive square well depths ε 11 , ε 22 ; and (iii) the corresponding interaction ranges ω

Table 1
Model parameters used to analyze experimental data of Wang 39 et al.In agreement with crystallographic data we use σ 2.95 nm.