Open Access Article
This Open Access Article is licensed under a
Creative Commons Attribution 3.0 Unported Licence

Effective patchiness from critical points of a coarse-grained protein model with explicit shape and charge anisotropy

Jens Weimar , Frank Hirschmann and Martin Oettel *
Institute for Applied Physics, University of Tübingen, Auf der Morgenstelle 10, 72076 Tübingen, Germany. E-mail: martin.oettel@uni-tuebingen.de

Received 15th July 2024 , Accepted 7th October 2024

First published on 8th October 2024


Abstract

Colloidal model systems are successful in rationalizing emergent phenomena like aggregation, rheology and phase behaviour of protein solutions. Colloidal theory in conjunction with isotropic interaction models is often employed to estimate the stability of such solutions. In particular, a universal criterion for the reduced second virial coefficient at the critical point image file: d4sm00867g-t1.tif is frequently invoked which is based on the behavior of short-range attractive fluids (Noro–Frenkel rule, image file: d4sm00867g-t2.tif). However, if anisotropic models for the protein–protein interaction are considered, e.g. the Kern–Frenkel (KF) patchy particle model, the value of the image file: d4sm00867g-t3.tif criterion is shifted to lower values and explicitly depends on the number of patches. If an explicit shape anisotropy is considered, as e.g. in a coarse-grained protein model, the normalization of image file: d4sm00867g-t4.tif becomes ambiguous to some extent, as no unique exclusion volume can be defined anymore. Here, we investigate a low-resolution, coarse-grained model for the globular protein bovine serum albumin (BSA) and study effects of charge-anisotropy on the phase diagram (determined by simulations) at the isoelectric point. We present methods of assigning an “effective patchiness” to our protein model by comparing its critical properties to the KF model. We find that doubling the native charges increases the critical temperature Tc by ≈14% and that our BSA model can be compared to a 3 to 5 patch KF model. Finally, we argue that applying existing image file: d4sm00867g-t5.tif criteria from colloidal theory should be done with care, due to multiple, physically plausible ways of how to assign effective diameters to shape-anisotropic models.


1 Introduction

Liquid–liquid phase transitions (LLPS) in protein solutions are a common phenomenon, they are related to biological function,1 but they also play an important role in the pathogenesis of some human diseases2,3 and also affect pathways to crystallization.4,5 Concepts from soft matter science (liquid state and colloidal theory) have been useful to see general and universal trends in protein phase behavior.6 This has been possible despite the fact that proteins are in general very anisotropic in shape and mutual interactions. Furthermore, interactions can be specific and non-specific and may sensitively depend on solvent conditions, and additionally proteins are flexible to a highly varying degree. In the simplest approach, isotropic colloidal models (hard spheres with short-ranged attractions) are used to explain features of protein phase diagrams. Notably these are the metastability of liquid–liquid demixing with respect to crystallization and an overall criterion for the necessary strength of attraction for the onset of this demixing. The latter has been proposed in ref. 7 and 8 and is based on an approximately universal phase behavior of short-range attractive colloidal models. It is formulated in terms of a reduced second virial coefficient image file: d4sm00867g-t6.tif at the critical temperature which takes the approximate value image file: d4sm00867g-t7.tif for the investigated colloidal models and has been suggested also to hold for protein solutions. In the following we refer to this finding as “Noro–Frenkel rule” or “law of corresponding states”.

A comprehensive study from the group of the late Stefan Egelhaaf9 finds a rather good agreement with the Noro–Frenkel rule if one scales the low-concentration part of the protein binodal with the corresponding part from a short-range square well fluid. On the other hand, the anisotropy of the attractive interactions between proteins frequently gives rise to clustering behavior which is reminiscent of the one seen in patchy particle models. For the latter (provided the interactions are only patchy and hard core), however, the location of the gas–liquid critical point is not consistent with the Noro–Frenkel rule and mainly the number of patches determine the reduced second virial coefficient, e.g.image file: d4sm00867g-t8.tif for three patches, image file: d4sm00867g-t9.tif for four patches etc.10 Thus one would expect also a deviation of image file: d4sm00867g-t10.tif from the Noro–Frenkel value for protein solutions if there is enough anisotropy in the interactions. Real protein–protein interactions are neither purely isotropic nor purely patchy, but always a combination of both. A study about patchy square well fluids, which represent such a combination, did not observe a law of corresponding states, i.e. in the form of image file: d4sm00867g-t11.tifvs. patch number.11 This implies that using universal image file: d4sm00867g-t12.tif criteria should only be applicable if the interactions of the model at hand are considered either mostly isotropic or mostly patchy.

Anisotropy is of course fully present if an atomistic protein model in explicit solvent is used. From a simulation point of view, it is still difficult to investigate phase separation in all-atom models for biomolecules and currently is restricted to small proteins and peptides.12 Here it is better to resort to coarse-grained (CG) models in which the protein is built from pseudo-atoms that unite multiple atoms into a single, bigger one. There exist numerous coarse-graining approaches (for reviews see e.g. ref. 12–14), with each of them operating at a certain resolution level. High-resolution CG models like the well-known MARTINI15 force field are often used with explicit water (e.g. for intrinsically disordered proteins16 or nanoparticles17), whereas low-resolution models in which only a few effective beads represent the molecule interactions are modeled by implicit-solvent potentials, as e.g. of DLVO type. These low-resolution models have been widely employed in the literature, e.g. in studies on antibodies18–21 but also for more globular proteins.22,23

Here we investigate the phase behavior of a coarse-grained model for the globular protein bovine serum albumine (BSA) consisting of 6 beads, using grand canonical Monte Carlo (GCMC) simulations. Interactions between beads are of DLVO type, the total BSA charge is zero but there is a fixed distribution of relative bead charges with one multiplicative parameter scaling the individual bead charges. In this model, the shape anisotropy (coming from the arrangement of the beads) determines also an asymmetry of the van der Waals attractions. Additionally, the bead charges add a further charge asymmetry to these interactions. We study the influence of charge on the critical point and the bead–bead correlations and attempt to quantify an effective patchiness from the image file: d4sm00867g-t13.tif-value for the investigated system.

The motivation for studying a BSA model is linked to the existence of LLPS of BSA in solutions with trivalent cations.24 On the one hand, these LLPS have been rationalized using patchy particle models,25,26 but on the other hand, scattering studies near the critical point in these systems show only small deviations from the Noro–Frenkel rule, thereby suggesting a small influence of patchiness on the phase diagram. Furthermore, BSA is well characterized by scattering methods in terms of its shape, dynamics and diffusion properties,27–31 it appears to be a rather rigid molecule suitable for modelling.

The paper is structured as follows: we first introduce our coarse-grained BSA model and its intermolecular interactions in Section 2, where in total six different parameter combinations are defined. Subsequently, in Section 3 the determination of critical and binodal points via GCMC simulation is described. In Section 4 the resulting phase diagrams are shown, along with the effective potentials at criticality, the second virial coefficients and effective molecule diameters, which are calculated via four different approaches. Additionally, we present bead–bead correlation functions and discuss their charge dependence. In Section 5 we discuss how an “effective patchiness” can be deduced from the critical point of our models with paying special attention to the ambiguity in defining effective hard diameters and thus in the normalization of the second virial coefficient (needed for defining image file: d4sm00867g-t14.tif). Furthermore the relationship of our results to experiments is discussed, before we conclude in Section 6.

2 Model description

The BSA model investigated in this work has previously been derived and analysed via Brownian dynamics (BD) simulations in ref. 32. In brief, the authors conducted single-molecule all-atom MD simulations of a crystal-structure of BSA and subsequently used direct Boltzmann inversion to obtain coarse-grained force field parameters, yielding a flexible, six-beaded, low-resolution CG model, in which each bead represents a biological domain, see Fig. 1(a) for illustration. The study assessed the influence of flexibility in self-crowded conditions on static and dynamic solution properties like center-of-mass pair-correlations and diffusion. However, the inter-molecular interactions were taken to be purely repulsive and soft, accounting only for excluded volume of the molecules, corresponding to a high-salt regime, where electrostatic interactions are sufficiently screened.
image file: d4sm00867g-f1.tif
Fig. 1 (a) Atomistic structure of BSA29 superimposed onto the domain-based, six-bead coarse-grained model of ref. 32 with bead diameter σ = 29.52 Å. Missing hydrogens have been added with the pdb2gmx33 tool. (b) Resulting coarse-grained charge distribution in units of the electron charge e of BSA at the isoelectric point (pH ≈ 5.8). Figures have been rendered using the VMD software.34

In this work, we consider more realistic interactions between molecules by incorporating dispersion forces and screened electrostatics which give rise to gas–liquid phase transitions and coexistence. More specifically, we set the interaction potential (bead–bead potential) between two coarse-grained sites on two distinct molecules to:

 
Utot(rij) = Urep(rij) + UvdW(rij) + Uel(rij) + Ucont(rij)(1)
where indices i and j refer to the ith and jth bead on the first and the second molecule (respectively) and rij represents the center-to-center distance of these beads. The individual contributions are as follows.

Volume exclusions are modelled by the repulsive part of the Lennard-Jones potential Urep, which is cut-off at its minimum:

 
image file: d4sm00867g-t15.tif(2)
where ε is the energy scale. Dimensionless energies U* and temperatures T* can be defined by U* = U/ε and T* = kBT/ε, where kB is the Boltzmann constant and T is the temperature. The asterisks are dropped immediately and for the remainder of this work, energies and temperatures are to be understood in their respective reduced units. The cutoff distance is rcutrep = 21/6σ and the bead diameter σ = 29.52 Å has been derived in ref. 32.

Dispersion interactions are modelled by the van der Waals potential UvdW for two equally sized spheres with:35

 
image file: d4sm00867g-t16.tif(3)
where AH is the Hamaker constant (in units of ε) and the cutoff distance has been set to rcutvdW = 2.5σ.

Electrostatic interactions are approximated via a screened Yukawa potential Uel:

 
image file: d4sm00867g-t17.tif(4)
where λB is the Bjerrum length, Zi and Zj are the electric charges on bead i and j (respectively, in units of the elementary charge e) and κ−1 is the Debye screening length. In order to have a noticeable influence of charges, we choose a moderate screening length of κ−1 = σ, which corresponds to a monovalent salt concentration of c ≈ 11 mmol L−1. Moreover, for simplicity, we ignore the temperature dependence of the Bjerrum length and set it to the textbook value for water at room temperature, λB ≈ 0.7 nm. Furthermore, rcutel = rcutvdW = 2.5σ. The charge distribution (i.e. Zi in eqn (4)) of all coarse-grained molecules within simulations is identical and determined through processing the atomistic structure of BSA (PDB ID: 4F5S, reported in ref. 29) with the APBS webserver,36 employing the PROPKA model37 to calculate the titration states of amino acids at the isoelectric point at pH ≈ 5.8 (experiments determine a slightly deviating value of pH ≈ 4.7).38 Subsequently, the partial charges of all amino acid residues belonging to the same coarse-grained site are summed up to calculate the total charge of the bead, see Fig. 1(b) for the resulting charge distribution. We deliberately chose our molecules to have an overall neutral net charge, since in experiments for various globular proteins in solution with multivalent ions, aggregation of molecules and liquid–liquid phase separation (LLPS) is observed around the isoelectric point.39,40 Additionally, in order to artificially alter the charge anisotropy of the CG molecules, we introduce a charge multiplier (CM), which uniformly scales all bead charges up or down according to:
 
Z → CM[thin space (1/6-em)]Z,(5)
such that CM = 0 would result in molecules with no electric charges on all beads and CM = 2 would result in a model where all charges depicted in Fig. 1(b) are doubled.

The last term Ucont in eqn (1) compensates the usage of cutoffs by making the total potential continuous. It reads:

 
image file: d4sm00867g-t18.tif(6)

For illustration, the complete interaction potential specified by eqn (1) is depicted in Fig. 2 for some typical parameters.


image file: d4sm00867g-f2.tif
Fig. 2 Total bead–bead potential according to eqn (1) for typical parameters.

In our simulations the CG molecules are rigid bodies without internal degrees of freedom. Relative bead positions are fixed, corresponding to the average structure present in coarse-grained BD simulations at infinite dilution, as described in ref. 32. In this study, we investigate different parameter compositions, i.e. the charge anisotropy is varied by choosing CM ∈ [0, 1, 2] and it is either AH = 4 or T = 1 fixed (while keeping either T or AH variable). Overall, this yields six different systems for which phase behavior can be studied, see Table 1 for an overview. In the simulations of the AH = const.-systems a physical critical temperature is determined. In contrast, in the T = const.-systems one determines the needed bead–bead attraction strength for criticality at ambient temperatures.

Table 1 Summary of investigated system parameters and results for the critical Hamaker constants AH,c, temperatures Tc, densities ρc (in units of σ−3) and B2 values (in units of σ3). Also shown are critical temperatures Tisoc, densities ρisoc and B2 values of the effective, isotropic systems (see Section 4.2). We estimate the relative statistical error for Tc and AH,c to be ±0.1% and for ρc ±0.4%. The relative total error of B2 is ±0.7%
System CM A H T A H,c T c ρ c A H,c/kBTc B 2(Tc) T isoc ρ isoc B 2 (Tisoc)
I 0 4 Variable 4 0.405 0.059 9.88 −28.9 1.12 0.18 −6.41
II 0 Variable 1 8.59 1 0.061 8.59 −27.5 2.95 0.21 −5.63
III 1 4 Variable 4 0.411 0.059 9.74 −28.6 1.02 0.17 −7.29
IV 1 Variable 1 8.46 1 0.061 8.46 −27.4 2.61 0.19 −6.55
V 2 4 Variable 4 0.470 0.054 8.51 −35.6 1.00 0.12 −10.7
VI 2 Variable 1 7.44 1 0.054 7.44 −35.7 2.25 0.13 −9.9


The value of AH = 4 appears to be a reasonable estimate for the bead–bead attraction strength for BSA solutions with monovalent salt, since an experimental structure factor at pH ≈ 7.0, cBSA = 174 mg mL−1 and cNaCl = 18 mmol L−1 at room temperature41 is roughly similar to the center-of-mass structure factor from our simulations at equivalent conditions (the methodology for obtaining experimental structure factors is described in ref. 28). On the other hand in the experiments with BSA solutions with trivalent salts the attraction strength between molecules is a tunable parameter at constant temperature,25 in our model this would be equivalent to a variable bead–bead attraction strength AH and a variable CM.

3 Simulation methods

The CG molecules were simulated in a grand canonical Monte Carlo scheme with textbook42 insertion, deletion and roto-translation moves (with 5%, 5% and 90% attempt probability, respectively). Orientations of molecules were represented via quaternions.43 The simulations were performed inside a cubic simulation box of side length L = 10σ. The main goal of our calculations is to obtain the critical points of the BSA systems. We additionally compute some binodal points to estimate the shape of the phase diagram. In this respect a finite size study was not conducted, since the expected finite size effects on the critical points are estimated not to be significant for the purpose of the present work.42 However, we expect that our binodal points are affected by finite size effects to some extend, especially those close to the critical points.

For obtaining the binodal points in the finite system close to Tc, histograms of the instantaneous molecule number N in the simulation box at different temperatures and values for the chemical potential μ were recorded and scanned for the appearance of double peaks. At the coexistence value for the chemical potential the area below these peaks must be equal. In order to obtain equal areas a fitting and optimization process using histogram reweighting44–46 was employed. Reweighting towards equal area gave new values for the chemical potential at which new histograms were simulated, whose peaks are possibly modified compared to the reweighted histogram. The procedure was repeated until the area difference below the peaks was below a chosen limit. The final peak positions give the average number of molecules in the coexisting gas and liquid states from which the corresponding number densities ρ = N/L3 are determined (binodal points).

Critical points of the different BSA systems (listed in Table 1) were determined by matching their double-peak histograms to the universal Ising distribution, following the methodology of ref. 42 and 47. In short, near the critical temperature and the critical chemical potential, where the histograms already looked similar to the universal Ising distribution, more simulations with slightly adjusted μ and T were performed. Second, by performing histogram reweighting on this data for new values μ′, T′ and the mixing parameter s (which connects the system energy E and N to the Ising-like order parameter X = NsE) as fit parameters, we determine the critical temperature Tc as the value of T′ which gives the best fit result on the universal Ising distribution. In systems where the critical Hamaker constant was determined, the simulation and reweighting process was performed with different values of AH (estimated by searching for equal area in the histograms) until the fitted value for the critical temperature Tc was very close to unity and then AH,c was computed as AH,c = AH/Tc (the correction due to Tc ≠ 1 was always below 0.2%). Fig. 3 shows fits on the universal Ising curve for the six different parameter combinations. The resulting critical densities ρc, temperatures Tc and Hamaker constants AH,c are listed in Table 1. Bead–bead correlations (see Section 4.4) were recorded in a larger box with the same GCMC code but with insertion and deletion moves disabled, which corresponds to a canonical NVT ensemble.


image file: d4sm00867g-f3.tif
Fig. 3 Fits of the order parameter X on the universal Ising curve for all six different subsystems. X connects particle number and energy through a mixing parameter s as X = NsE and was scaled on the Ising curve (zero mean, unitary variance, cf. ref. 42 and 47). (a) shows fits for systems with constant AH = 4 (with the critical temperature Tc as fit parameter), and (b) shows fits for the systems with a constant temperature T = 1 (with the critical value of the Hamaker constant AH,c as fit parameter).

For obtaining the critical points 50 to 60 independent systems have been simulated, each of them with 109 GCMC moves in total. For the binodal points and the bead–bead correlations 10 resp. 80 independent systems have been employed, with 1 × 109 to 2 × 109 and 109 GCMC moves, respectively. We also note that obtaining binodal points proved difficult, especially for low temperatures, due to systems remaining in one phase even for long simulation runs. To this end, more elaborate simulation techniques could be applied, e.g. successive umbrella sampling,48 which was out of the scope of this study.

4 Results

4.1 Phase diagrams and critical points

Our results for the binodals of the three differently charged systems with CM = 0, 1, 2 and AH = 4, as well as the three systems with T = 1 (see Table 1) are depicted in Fig. 4.
image file: d4sm00867g-f4.tif
Fig. 4 Binodals of differently charged systems (see Table 1) with fixed Hamaker constant AH = 4 (a) and fixed temperature T = 1 (b). Critical points are indicated by colored stars and state points where bead–bead correlation functions gij(rij) have been recorded are indicated by colored crosses (see Fig. 8 in Section 4.4).

In general, the resulting phase diagrams obtained for systems with either variable T or variable AH are in shape very similar to each other, but not exactly identical, since even for the completely uncharged systems the Hamaker constant with our definition is not just an inverse temperature. Starting from the uncharged molecules (CM = 0) and adding the native charges (CM = 1) to our molecules introduces additional attractive interactions to the systems and shifts Tc upwards by (only) ≈1.5% (AH,c downwards by ≈1.5%), and also inflicts minute changes to the slope of the coexistence curves. However, more pronounced effects are seen if one changes CM from 1 to 2, as this shifts Tc upwards by another ≈14% and AH,c downwards by ≈12%. We attribute these additional attractions to the growing, leading order, electrical dipole moment of our molecules, and note that this physical quantity has been observed to play an important parameter for loci of coexistence regions in protein-like fluids.49 Overall, our findings indicate that the stability in our BSA systems for CM = 1 is mainly governed by the van der Waals interactions and not so much by electrostatics.

From the results in Table 1 we infer that the effective temperature-like variable kBT/AH is not even approximately the same at the critical point for the different systems (AH = const. and T = const.) at the same CM. The two types of systems cannot be mapped onto each other due to the soft repulsions. However the dependence of kBT/AH on CM is similar, and indeed it turns out that only CM is a relevant variable (see Section 4.3).

4.2 Effective potentials at the critical points

An isotropic effective potential Ueff(r) between two molecules can be defined in the usual way by the angular average:50
 
image file: d4sm00867g-t19.tif(7)
where r is the center-of-mass distance between two coarse-grained molecules, Ωi specifies the orientation of molecule i (e.g. by three Euler angles), U(r, Ω1, Ω2) is the total orientation-dependent pair potential between two molecules (in our case image file: d4sm00867g-t20.tif) and the integral over orientations is normalized image file: d4sm00867g-t21.tif. If U does not depend on orientations, then U = Ueff of course. We calculate the six-dimensional integral via two methods: first, we employ Monte Carlo sampling with 108 uniformly distributed orientations per distance r with a spatial resolution of Δr = 0.01σ. Additionally, we compute eqn (7)via Gauss–Legendre quadrature and 20 grid points per angle. The obtained Ueff show negligible differences, the resulting B2 values differ on average by ≈0.22% for the two methods. For the Gauss–Legendre method using 15 grid points instead of 20 yields a relative error of obtained B2 values of ≈0.5%.

Through examination of molecule–molecule configurations in eqn (7) the lowest-energy state of two molecules can be determined: for the CM = 0, 1 systems this is a very compact state at r = 1.2, which maximizes bead contacts. For CM = 2 it is a staggered state at r = 1.4, in which one molecule sits on top of the other but slightly rotated, such that bead 4 is near bead 3. These states, however, do not contribute to Ueff in a noticeable manner. Consistent with that is the absence of a noticeable amount of compact dimers in simulation snapshots. In these, dimers are on average further apart with varying bonding configurations.

We utilize the obtained effective potentials in additional GCMC simulations, in which isotropic particles interact via the six critical Ueff(r) obtained by the previously described procedure. Subsequently, via the same protocol as described in Section 3 we determine critical temperatures Tisoc for these systems and the obtained values are listed in Table 1. In the following sections of this paper these systems are referred to as “isotropic systems”. The obtained Tisoc are approximately two to threefold higher as the corresponding Tc of the anisotropic models, pointing to a significant over-representation of attractions and/or under-representation of repulsions in Ueff(r). An example are higher-order multi-body interactions not being present in Ueff(r): it could be favorable for two CG molecules to form a low-energy dimer, which enters the effective potential. In the isotropic system three spheres can easily arrange to a trimer, with a threefold binding energy of the dimer, whereas three CG molecules can not easily arrange to a trimer with an equally low binding energy due to steric incompatibilities.

4.3 Second virial coefficients and effective molecule diameters at the critical points

Using the effective potential Ueff(r) of eqn (7), the second virial coefficient can be written as:
 
image file: d4sm00867g-t22.tif(8)
The obtained values at the respective critical points of all models are depicted in Fig. 6(a) and listed in Table 1. Increasing CM from 1 to 2 introduces additional attractions between molecules, as B2 noticeably decreases, whereas almost no influence is visible if comparing the uncharged molecule to the native CM = 1 state. Within anisotropic systems sharing the same CM no significant difference is visible. The values for the isotropic systems on the other hand show a stronger charge dependence and values of systems with CM = 0, 1 are more separated to systems with CM = 2. Note that the reason for the B2 values of the isotropic systems being higher than for the anisotropic ones, is that eqn (8) is evaluated at T = Tisoc and that Tisoc/Tc ≈ 2…3, as pointed out before.

For obtaining reduced second virial coefficients image file: d4sm00867g-t23.tif it is necessary to determine an exclusion volume which represents steric incompatibilities, i.e. an effective spherical diameter σeff of the molecule in terms of which B2,ex = (2π/3)σeff3. Here, we present four approaches with a certain theoretical preference for the first one and three additional methods for comparison, which however also appear natural from a physical point of view.

(i) σNFeff and σNF,isoeff. In the original Noro–Frenkel (NF) analysis of various short-range attractive, isotropic systems8σNFeff is computed from the repulsive part (vrep(r)) of the isotropic pair potential (determined by the WCA rule, with the pair potential v(r) cut-off at the minimum and shifted such that the minimum is zero) via the Barker–Henderson (BH) rule image file: d4sm00867g-t24.tif.51 In our case, we use the effective potential for the anisotropic systems, i.e. v(r) → Ueff(r). For the isotropic systems σNF,isoeff is calculated exactly the same way, only with T = Tisoc (see Fig. 6(b) and Table 2 for resulting values).

(ii) σrepeff. Excluded volume results from the repulsive part of the bead–bead potentials. First we define a bead–bead potential only from repulsive interactions, Utot(rij) = Urep(rij), as defined in eqn (2). Then Ueff is computed from eqn (7). The resulting repulsive Ueff is charge-independent (if evaluated at the same T) and it is shown in comparison to the full Ueff in Fig. 5 (black, dashed line). Like before, from the BH rule the diameter σrepeff of an equivalent hard sphere system can be determined. The values σrepeff ≈ 2.1σ (see Fig. 6(b) and Table 2) are much larger than σNFeff from the effective potential, this is also visible in the comparison of the effective potentials in Fig. 5.


image file: d4sm00867g-f5.tif
Fig. 5 Effective potential βUeff(r) (eqn (7)) obtained via MC integration for the critical systems with Tc = 1 (colored) and for a purely repulsive model for which AH = CM = 0 and T = 1 (black, dashed line). Effective potentials of the critical system with AH,c = 4 (not shown) are very similar to the ones depicted with the same CM.

image file: d4sm00867g-f6.tif
Fig. 6 (a) B2 values according to eqn (8) of all critical BSA models presented in Table 1, together with B2 values of the corresponding critical isotropic systems. (b) Effective molecule diameters σeff at the critical points, calculated via the four different approaches (i)–(iv) discussed in the main text.

(iii) σdryeff. In ref. 32 the bead diameter σ = 29.52 Å was derived by distributing the molecule volume of BSA on six equally sized spheres. The molecule volume was obtained by computing the protein specific volume with the molecular mass of BSA. Here, we conduct the calculation in reverse: first, an effective bead diameter σbeadeff is determined by applying the BH rule to the repulsive part of the bead–bead potential Urep(rij) (eqn (2)), its value is always slightly above 1.0σ (see Table 2). Then, the 6 beads of the model constitute the occupied, “dry” volume 6·(π/6)(σbeadeff)3 of the protein without hydration shell, which is approximately the BSA volume in protein powder form (and is nearly identical the molecule volume calculated in ref. 32). The dry volume can be converted to an effective sphere diameter σdryeff of an equivalent sphere with the same volume, i.e. σdryeff = 61/3σbeadeff. Hence, in contrast to the previous methods no Ueff(r) is computed here.

(iv) σcpeff. Having calculated the critical molecule number densities ρc of our systems (see Table 1), we can define an effective volume packing fraction ϕceff at the critical point via ϕceff = (π/6)(σcpeff)3ρc. Subsequently, we equate ϕceff = ϕcSW, where ϕcSW = 0.249 is the well defined volume packing fraction at the critical point of a short-range, hard sphere square well fluid (taken from ref. 9, Table 2 in that reference and the system with reduced well width δ = 0.05), which is approximately representative of the short-range potentials to which the Noro–Frenkel rule applies. This analysis yields σcpeff ≈ 1.98…2.06σ, see Fig. 6(b) and Table 2. Furthermore, we calculate the reduced second virial coefficients for the binodal points presented in Fig. 4, using ϕeff = (π/6)(σcpeff)3ρ and B2,ex = (2π/3)(σcpeff)3. The resulting phase diagrams in the image file: d4sm00867g-t25.tif plane are shown in Fig. 7, along with the phase diagram of the aforementioned, short-range, reference square well system (black crosses). The packing fractions of our binodals approximately match the binodals of the reference, which is in effect similar to the construction9 used by the group of Stefan Egelhaaf mentioned in the introduction and further discussed in Section 5.4. Note that if a longer ranged hard sphere square well fluid with δ ≈ 0.25 would be used as reference, then ϕcSW ≈ 0.2,9 and therefore σcpeff would decrease by approx. 7%.


image file: d4sm00867g-f7.tif
Fig. 7 Reduced second virial coefficient image file: d4sm00867g-t26.tif with B2,ex = (2π/3)(σcpeff)3 as a function of the volume packing fraction ϕ = ϕeff = (π/6)(σcpeff)3ρ for the six different parameter combinations. The cyan shaded area represents the range of critical image file: d4sm00867g-t27.tif values for short-range (δ ≤ 0.25) and long-range (δ ≥ 0.25) hard sphere square well fluids calculated by Gazzillo et al.52 Black crosses represent the reference hard sphere square well binodal taken from Platten et al.9 (Fig. 2 and Table 2 in that reference), which was used to obtain σcpeff as described in Section 4.3 by matching the critical points. Dashed lines are guides to the eye.

Fig. 7 also shows B2/B2,ex-ranges in which short-range (SR, δ ≤ 0.25) and long-range (LR, δ ≥ 0.25) hard sphere square well systems are critical (cyan shaded area).52 This total range is nearly identical to the range of the isotropic Noro–Frenkel fluids depicted in Fig. 9(a) (also cyan shaded area), but allows for a distinction between different ranges of attraction.

All previously discussed effective diameters are listed in Table 2 and shown in Fig. 6(b). Evidently, method (i) yields the smallest σeff compared to the other methods, whereas methods (ii)–(iv) yield higher values, which are generally relatively close to each other for all investigated systems. Systems with constant T roughly yield identical results as systems with constant AH for the same CM. Adding charges to the molecules significantly increases their effective diameter only for method (i), as the increased electrostatic repulsion at short distances is directly present in Ueff, whereas methods (ii) and (iii) are (for the systems with variable temperature) only indirectly sensitive to the parameter CM, i.e. solely via the shifted inverse temperature βc which enters the BH rule. The latter is also responsible for all σNF,isoeff being smaller than their respective σNFeff, as Tisoc/Tc ≈ 2…3.

Note that methods (i) and (ii) are in spirit similar to each other, as they derive the effective diameter only from the two-body interaction and the pair potential we attribute to the model. Method (iii) is based on an experimentally measured quantity, while method (iv) is based on the critical point, a collective property.

4.4 Bead–bead correlations

In order to evaluate how charge-anisotropy influences the preferred bead contacts in our simulations we compute bead–bead correlation functions gij in Fig. 8 (note that gij = gji) for our coarse-grained molecules. We choose a relatively low density ρ = 5.75 × 10−3σ−3 outside the coexistence region in the gas region of the phase diagram (colored crosses in Fig. 4(b)) to mitigate packing effects, which would complicate a separation with effects stemming from possible patchy interactions.
image file: d4sm00867g-f8.tif
Fig. 8 Bead–bead correlation functions gij(rij) with rij representing the center-to-center distance of bead i to bead j, for systems with constant temperature T = 1 and variable charge ((a): CM = 0, (b): CM = 1, (c): CM = 2) at a dilute molecule number density ρ = 5.75 × 10−3σ−3 (for the locations in the phase diagram see the colored crosses in Fig. 4(b)). Only selected (i, j) combinations that exhibit the highest and lowest contact peak values are shown. Attractive, repulsive and neutral charge configurations are indicated in the legend. Error bars show the standard deviation of the mean.

Starting from the uncharged molecule state CM = 0 one notices first-peak contacts are most pronounced for combinations (i, j) of exposed beads at the edges of the molecule, i.e. beads 1, 4 and 6 and least pronounced for beads that are closer to the center of the molecule, i.e. beads 2, 3 and 5 (compare to Fig. 1). Adding the native charges to the BSA models (CM = 1) does not introduce new (i, j) combinations to the top first-peak bead–bead correlations, indicating that for this system the van der Waals and not the electrostatic interaction determines association of the molecules. This is additionally supported by the observation that only bead combinations with like charge are within the top correlation functions. However, increasing CM to 2 finally switches the dominant interaction mechanism from van der Waals to electrostatics, as now the strongest charged beads exhibit either the highest first-peak correlations, i.e. (3, 6) and (3, 4) if oppositely charged, or the lowest first-peak correlations, i.e. (3, 3) and (3, 5) if like charged. Intriguingly, despite the electrostatic dominance the bead–bead combinations (4, 4) and (1, 4) favoured by the molecule's geometry combined with the dispersion force still exhibit the third and fourth highest correlation functions. We thus conclude that for our model there is a crossing-over between CM = 1 and CM = 2, where the mechanism of preferred contact selection between beads switches from van der Waals to electrostatic interactions, where the former's influence still persist to some extent. Additionally, we infer that this analysis provides arguments that our BSA model effectively exhibits three patches (within a KF picture), due to the three sites at the edges of the CG molecule (beads 1, 4 and 6), that show preferred contact propensity in gij.

5 Discussion: effective patchiness from critical points

5.1 Reduced second virial coefficient image file: d4sm00867g-t28.tifvia method (i)

In their seminal paper,8 Noro and Frenkel have shown that systems with isotropic, short-range attractive pair potentials can be mapped onto hard sphere square well systems (with the same B2 and the effective hard sphere diameter σNFeff determined with the Barker–Henderson rule), with the effect that the reduced virial coefficient image file: d4sm00867g-t29.tif at the critical temperature is approximately constant over quite a range of pair potentials and given by image file: d4sm00867g-t30.tif. The Noro–Frenkel rule does not hold for patchy particles, though. Foffi and Sciortino10 have shown that for Kern–Frenkel (KF) particles53 (hard spheres with attractive surface patches defined via narrow cones with opening angle θ0 and with reduced attraction range δKF) with M = 3, 4, 5 patches image file: d4sm00867g-t31.tif is approximately a function of M only, but the actual values are very different from image file: d4sm00867g-t32.tif (≈−30 for M = 3, ≈−4…−5 for M = 4 and ≈−2…−3 for M = 5).

In the KF model, let εKF be the attractive interaction strength of two interacting KF patches. The corresponding effective potential is a square well potential with well depth ε′, which is related to εKF by:

 
exp(βε′) = χ2(exp(βεKF) − 1) + 1,(9)
where χ = (M/2)(1 − cos[thin space (1/6-em)]θ0) < 1. It is interesting to ask for the relation of critical temperatures in the orientation-dependent KF system vs. in the corresponding effective square well system. For the inverse critical temperature βc,KF of KF we use the data taken from ref. 10. The inverse critical temperature image file: d4sm00867g-t33.tif of the effective square well system can be safely estimated using the Noro–Frenkel rule and the reduced B2 in the square well system, i.e.:
 
image file: d4sm00867g-t34.tif(10)
where δ is the range of the square well attraction in units of the hard sphere diameter. Here, ε′ is determined from εKF from eqn (9), with β = βc,KF (the inverse critical temperature of the KF system) and we set δ = δKF. From the data in ref. 10 for various KF systems we find a ratio image file: d4sm00867g-t35.tif (M = 3) and image file: d4sm00867g-t36.tif (M = 4) and image file: d4sm00867g-t37.tif (M = 5). It is clear that the critical temperature of the effective system (particles interacting with Ueff) is always noticeably higher than the one of the KF system, but the ratio is not only a function of the number M of patches. There is a clear separation between the values for all M = 3 systems from the rest, but there is a considerable overlap in the image file: d4sm00867g-t38.tif-ranges for M = 4 and 5, see Fig. 9(b). All values are, however, significantly smaller than 1 which is the limit of isotropic potentials.


image file: d4sm00867g-f9.tif
Fig. 9 (a) Reduced second virial coefficients image file: d4sm00867g-t39.tif = B2/B2,ex, with B2,ex = 2πσeff3/3 and where σeff is taken as σNF(,iso)eff, σrepeff, σdryeff or σcpeff (see Section 4.3 and Fig. 6(b)) for the anisotropic, critical BSA systems together with image file: d4sm00867g-t40.tif values of the corresponding critical isotropic systems, as well as regions in which isotropic fluids (Table 1 in ref. 8, cyan shaded area) and Kern–Frenkel particles with M = 4, 5 patches are critical (Table 1 in ref. 10, violet and green shaded areas). (b) Ratios image file: d4sm00867g-t41.tif for all critical BSA systems together with regions in which image file: d4sm00867g-t42.tif falls for different numbers of patches M = 3, 4, 5 in the Kern–Frenkel model (yellow, violet and green shaded areas). Latter regions have been calculated as explained in the main text with values taken from ref. 8 and 10.

Note that in the KF analysis above the anisotropy only arose from the patchy attractions, while repulsions are isotropic hard sphere in both the KF and the effective system. This is different in our system where excluded volume interactions are very anisotropic and the effective isotropic repulsions are very soft.

These considerations can be applied to our model since Ueff and the corresponding critical temperatures in the effective systems have been determined in Section 4.2, we refer to the results for σNF,isoeff and image file: d4sm00867g-t43.tif in Fig. 6(b) and 9(a) (respectively) and also in Table 2. For all investigated isotropic systems, image file: d4sm00867g-t44.tif which is different from image file: d4sm00867g-t45.tif. This is due to the rather long range in the attractions of Ueff and matches very well results in Noro and Frenkel's work8 for wide square wells: image file: d4sm00867g-t46.tif changes from −1.65 to −2.26 if the reduced radial range δ of the attractive square well changes from 0.5 to 1. The effective diameters σNF,isoeff ≈ 1.2…1.4σ for all systems (see Fig. 6(b)), and it is interesting to note that these diameters are only ≈20…40% larger than the bead diameter σ.

For the assessment of “patchiness”, the values of image file: d4sm00867g-t47.tif, as well as image file: d4sm00867g-t48.tif and σNFeff at the critical temperature of the anisotropic system are relevant. As can be seen in Fig. 9(b) and Table 2, the values of image file: d4sm00867g-t49.tif are around 0.34…0.47 with a moderate spread and all values of image file: d4sm00867g-t50.tif are between −4 and −5 with a likewise moderate spread. Regarding charge dependence, σNFeff exhibits a systematic increase for a rising CM value (see Fig. 6(b), corresponding to the observation of a larger repulsive core in Ueff). Hence, although the magnitude of B2 is growing with increasing CM (see Fig. 6(a)), the charge dependence of σNFeff (and consequently also of B2,ex) leads to the monotonic and moderate charge dependence in image file: d4sm00867g-t51.tif (similar to image file: d4sm00867g-t52.tif).

Comparing this to the KF results, these values correspond to 4 patches (from image file: d4sm00867g-t53.tif) or 3 patches for the CM = 2 systems (from image file: d4sm00867g-t54.tif), which points to a considerable effective anisotropy in the BSA model.

5.2 Alternative determinations of B2,exvia methods (ii)–(iv) and their influence on image file: d4sm00867g-t55.tif

The normalization of B2 with B2,ex = 2πσeff3/3, the second virial coefficient for the excluded volume interaction of the equivalent sphere system, leaves a certain ambiguity in the definition of B2,ex which becomes manifest in a system where the repulsive interactions are explicitly anisotropic. It is interesting to explore resulting image file: d4sm00867g-t56.tif if the additional σeff of Section 4.3 are used.

The values for σrepeff, σdryeff and σcpeff are around 1.9…2.1σ (approximately charge-independent) and relatively close to each other, see also Fig. 6(b) and Table 2. All three effective diameters are considerably larger than σNFeff with the immediate consequence that image file: d4sm00867g-t57.tif is reduced in magnitude and they all exhibit a stronger charge dependence than before: for method (ii) (with σeff = σrepeff) we find image file: d4sm00867g-t58.tif for the systems with CM = 0 and 1, and image file: d4sm00867g-t59.tif for the systems with CM = 2. For method (iii) (with σeff = σdryeff) the difference between CM = 0, 1 on the one side and CM = 2 on the other side can be seen as well: image file: d4sm00867g-t60.tif (CM = 0, 1) and ≈−2.5…−2.7 (CM = 2). Method (iv) yields approximately values that fall in the middle of (ii) and (iii): image file: d4sm00867g-t61.tif (CM = 0, 1) and image file: d4sm00867g-t62.tif (CM = 2). Compared with the values from the KF model, we would argue that according to methods (ii)–(iv) our BSA model is fairly isotropic for CM = 0, 1. However, the deviating values of image file: d4sm00867g-t63.tif for the CM = 2 systems are intriguing. Here method (iii) yields results compatible with a 5 patch KF model.

The ambiguity in defining B2,ex drops out if we form the ratio of critical B2-values

 
image file: d4sm00867g-t64.tif(11)
which yields ≈1.24…1.29. Using the results for different systems in Table 1 of ref. 10, we estimate for the KF model r ≈ 4.12…7.35 for the ratio between a 3-patch and a 4-patch model, and r ≈ 1.09…2.56 for the ratio between a 4-patch and a 5-patch model. The latter range includes our results, such that one could conclude that our BSA model is moderately patchy, around 4 patches for CM = 2 and 5 patches for CM = 0.

As a remark, we note that the term “patchiness” used in this study refers solely the kind present in the KF model, i.e. a purely attractive one. There are other kinds of interesting model systems without uniform van der Waals attraction exhibiting patches that repel each other, for example inverse patchy colloids54 (IPCs), lattice charge models55 or continuous charge patchy models.56 Such repulsive patches are obviously also present in our CG model, due to same-charge interactions, and they are overall less dominant than the opposite-charge interactions (otherwise Tc would not increase by increasing CM). By the choice of our KF/NF reference systems the effects of these two types of patches are mixed to an effective, attractive patchiness. Here, a systematic study that separates these effects would be intriguing, however, no law of corresponding states that would enable a investigation similar to this one has been observed for IPCs so far.57

5.3 Critical density ρc

Also the critical density ρc is sensitive to the number of patches within the KF model, but in contrast to image file: d4sm00867g-t65.tif the range of critical densities for long-range isotropic fluids overlaps with the corresponding range for KF fluids with 4 and 5 patches (see Fig. 10). The corresponding quantity in our model systems is the critical density (Table 1) evaluated in units of 1/σeff3, using the various definitions from Sections 5.1 and 5.2 (Fig. 6(b) and Table 2), see the symbols in Fig. 10. The analysis confirms the conclusions from Sections 5.1 and 5.2 from the analysis of image file: d4sm00867g-t66.tif, especially regarding the difference between σNFeff on the one hand and the alternative definitions (ii)–(iv) from Section 5.2. If the latter are used, the critical densities for the CM = 0 systems are in the range of isotropic short-range square well fluids and decrease slightly with increasing the charge. Only for σdryeff, the highly charged CM = 2 systems appear slightly patchy (5 KF patches), however, overlapping with the range of isotropic long-range square well fluids. Overall, this is comparable to the slight decrease of image file: d4sm00867g-t67.tif with increasing charge seen in Fig. 9(a) for the alternative definitions of σeff. The relative decrease of ρc (CM = 2) vs. ρc (CM = 0) of about 10% roughly corresponds to the relative decrease in ρc of KF 5-patch to 4-patch models (and would constitute an estimate independent of σeff). If however method (i) is used for the anisotropic systems, the critical densities for the CM = 0 systems are compatible with a 3-patch KF system and the charges tend to increase the critical density towards the one of a 4-patch KF system. This effect of charge is similar to the increase of image file: d4sm00867g-t68.tif for this definition of the effective diameter, see Fig. 9(a), the only difference being that according to the critical density analysis the CG model has an overall more pronounced patchy nature. The critical densities of the isotropic systems are in the range of isotropic long-range fluids (in agreement with the long-range character of Ueff in Fig. 5).
image file: d4sm00867g-f10.tif
Fig. 10 Critical densities ρc (Table 1) in units of 1/σeff3, where σeff is taken as σNF(,iso)eff, σrepeff, σdryeff or σcpeff (see Section 4.3 and Fig. 6(b)). The shaded areas correspond to the critical densities for KF particles reported in ref. 10 and to short- (δ ≤ 0.25) and long-range (δ ≥ 0.25) hard sphere square well systems reported in ref. 9.

5.4 Relation to experiment

Scattering experiments on BSA solutions in monovalent salt indicated stability of the BSA solution.28 Available experimental data for B2 in sodium chloride solution at pH = 7.4 and ambient temperatures were fitted in ref. 23 using a 20-bead model with DLVO forces, similar to ours. The resulting B2 values (in our unit σ3) in the range −5…25 are considerably larger than the B2 values at the critical point for our models (which are in the range −28…−37, see Fig. 6(a)). Since the fit of ref. 23 was for high salt concentrations and therefore effective charge screening, one should compare to the maximal value for the chargeless model (CM = 0), B2 ≈ −28. Thus one would conclude that a hypothetical critical temperature in a sodium chloride solution is way below the freezing point of water which is outside of physical interest.

In contrast to sodium chloride solutions, liquid–liquid phase separations have been observed in BSA solutions with trivalent cations.24 The necessary attraction has been proposed to be due to cation bridges at specific places on the BSA surface, i.e. patchy attractions in essence.25 The phase separation occurs near the point of zero charge. Measurements using small angle X-ray scattering (SAXS) of B2/B2,ex40,58–60 have indicated though that image file: d4sm00867g-t69.tif, which is of the order of image file: d4sm00867g-t70.tif and thus one would ascribe no patchiness to BSA. However, in these works the effective hard sphere diameter σeff used in B2,ex was fitted from an ellipsoidal form factor and gave σeff ≈ 9.2 nm, which is larger than our σeff by a factor 1.5…2.6. This particular ellipsoidal form factor fit is probably due to clusters. Note that a similar form factor fit in the SAXS measurements of ref. 28 (with monovalent salt, less attraction and thus less propensity for cluster formation) gave σeff ≈ 6.7 nm which is much closer to our σrepeff ≈ 6.3 nm which can be considered similar in spirit. Thus, using our range of σeff, image file: d4sm00867g-t71.tif in the BSA-trivalent cation system is estimated to be around −4.6…−32 which indicates much more patchiness and also a larger patchiness than found in our model (ranging between −1.4…−5.1). Our model has not been designed with very specific bridging attractions, rather we have put the effect of the attraction into an increase of AH and the attraction between positive and negative charges at the BSA surface. We would expect that image file: d4sm00867g-t72.tif attains larger negative values if specific bridging sites are included in our model.

A somewhat different perspective on the validity of the Noro–Frenkel rule for protein solutions has been given in ref. 9, discussing the phase diagrams of globular proteins, mainly of lysozyme. The determination of a suitable effective diameter and thus of image file: d4sm00867g-t73.tif uses the liquid–liquid binodal, which is in effect similar to method (iv) of Section 4.3. Specifically, the argument in ref. 9 concerning the validity of the Noro–Frenkel rule proceeds as follows: from an extended body of simulation data for binodals of isotropic fluids a “Noro–Frenkel band” of binodals is defined in the image file: d4sm00867g-t74.tif plane (packing fraction vs. reduced second virial coefficient). Considering now a protein solution, an effective diameter is defined by fitting the gas part of its binodal to the binodal of a square well fluid with reduced well width δ = 0.05 (this binodal is located approximately in the middle of the “Noro–Frenkel band”). As a result, the investigated lysozyme binodals nicely fall into this band and image file: d4sm00867g-t75.tif. Thus, ref. 9 argues, that the Noro–Frenkel rule is valid and would presumably also hold for other globular proteins, i.e. from critical liquid–liquid demixing no patchiness can be inferred. As a result, a such defined diameter is considerably larger than values derived from an effective, isotropic DLVO potential (as also seen in ref. 9). This is in agreement with our findings for σcpeff ≈ 2.0 (whose derivation is in effect similar to the construction of ref. 9), which is larger than σNF(,iso)eff ≈ 1.2…1.6 (from the orientation-averaged effective potential). Consequently, the critical points of all systems with varying CM are well within the region where isotropic fluids are critical, see Fig. 7. If using a finer distinction of image file: d4sm00867g-t76.tif values, that discriminates short and long-range hard sphere square well systems, the CM = 2 system clearly corresponds to longer ranged wells (δ ≥ 0.25), whereas the CM = 0, 1 systems fall in an intermediate region where δ ≈ 0.25.

6 Conclusion

We have investigated the phase diagram and critical points of an anisotropic coarse-grained model for the globular protein BSA. Besides the shape anisotropy through the spatial arrangement of the beads, we investigated also the dependence of the critical point on the charge asymmetry, realized by different bead charges at the isoelectric point. We have attempted to quantify the anisotropy effects on the critical point through an “effective patchiness” which is built on the suggestion by Foffi and Sciortino10 for an extended Noro–Frenkel rule8 for patchy particles. A difficulty arises from the fact that in the analysis of Foffi and Sciortino the repulsive core of the patchy particles are isotropic hard spheres whereas in our case the repulsive core is soft and anisotropic.

We have argued that the analysis of the effective, orientation-averaged pair potential Ueff provides a useful link: first, it allows to define an effective hard sphere diameter σNFeff and a repulsive part B2,ex of the second virial coefficient and secondly, the ratio of critical temperatures of an isotropic model interacting with Ueff and of the fully anisotropic model contains a definite patch dependence. The reduced second virial coefficient at the critical point image file: d4sm00867g-t77.tif points to an effective patchiness of 4 Kern–Frenkel patches image file: d4sm00867g-t78.tif whereas the ratio of critical temperatures points to a 3-patch model for CM = 2, see Fig. 9. Through computing the ratio r of B2 (eqn (11)), the dependence on B2,ex drops out and for our BSA model we see that the CM = 0 and CM = 2 charges lead to an effective anisotropy of approximately 4 (CM = 2) to 5 (CM = 0) Kern–Frenkel patches.

Through its dependence on B2,ex, the image file: d4sm00867g-t79.tif-criterion does not appear to be unique. We have investigated two other sensible criteria for defining B2,ex and σeff directly related to the molecular properties (through the repulsive part of the bead–bead potentials and the dry protein volume, yielding σrepeff and σdryeff) and obtain larger effective diameters and consequently also values image file: d4sm00867g-t80.tif (and larger) which are mostly compatible with the Noro–Frenkel rule for isotropic particles, except for σeff = σdryeff and CM = 2, see Fig. 9(a).

The analysis presented here shows that care is needed in the definition of an effective hard sphere diameter when working with a image file: d4sm00867g-t81.tif-criterion. This has become evident when (as a third alternative) applying an analysis similar to the one by Platten et al.9 on the phase behavior of lysozyme, yielding σcpeff. Here an effective hard sphere diameter was defined through matching the critical points of the proteins to a square well system, which again resulted in an approximate validity of the Noro–Frenkel rule.

It appears difficult to straightforwardly apply rules for isotropic systems with reasonably well defined excluded volume to anisotropic systems. Our work is an example how much the results for image file: d4sm00867g-t82.tif can vary when applying different methods which all appear reasonable. In this respect, computing the isotropic potential Ueff offers two advantages. First, the definition of a hard sphere diameter via Barker–Henderson integration (as applied for soft isotropic systems in ref. 8) and second the possibility of using the ratio of critical temperatures between the anisotropic and the effective isotropic system as a variable to infer patchiness.

Author contributions

Jens Weimar and Frank Hirschmann: data curation, formal analysis, investigation, methodology, software, visualization, writing – original draft, writing – review & editing. Martin Oettel: conceptualization, funding acquisition, methodology, project administration, resources, supervision, validation, writing – original draft, writing – review & editing.

Data availability

The data supporting our results is included in Table 2. The code and algorithms yielding our results are standard computational techniques (multi-dimensional integration, grand canonical Monte Carlo simulation and histogram reweighting) which are described in the references provided in our manuscript.
Table 2 Summary of critical parameters and effective molecule diameters σeff (in units of σ) at criticality calculated via the different methods (i)–(iv) presented in Section 4.3 for the anisotropic and isotropic systems. We estimate the total relative error for σeff and image file: d4sm00867g-t83.tif to be ±0.7% and ±2.6% (respectively) for methods (i) and (ii), ±0.9% and ±3.4% for method (iii) and ±0.4% and ±1.7% for method (iv)
System CM Anisotropic systems Isotropic systems
T c A H,c (i): σNFeff (ii): σrepeff (iii): σdryeff (iv): σcpeff T isoc

image file: d4sm00867g-t84.tif

(i): σNF,eff

image file: d4sm00867g-t85.tif

σ NFeff

image file: d4sm00867g-t86.tif

σ repeff

image file: d4sm00867g-t87.tif

σ dryeff σ beadeff

image file: d4sm00867g-t88.tif

σ cpeff

image file: d4sm00867g-t89.tif

σ NF,isoeff
I 0 0.405 4 −5.10 1.39 −1.41 2.14 −2.01 1.90 1.05 −1.70 2.01 1.12 0.36 −1.66 1.23
II 0 1 8.59 −5.09 1.37 −1.41 2.10 −2.09 1.85 1.02 −1.69 1.98 2.95 0.34 −1.61 1.19
III 1 0.411 4 −4.57 1.44 −1.40 2.14 −2.00 1.90 1.04 −1.69 2.01 1.02 0.40 −1.67 1.28
IV 1 1 8.46 −4.54 1.42 −1.41 2.10 −2.08 1.85 1.02 −1.67 1.98 2.61 0.38 −1.62 1.24
V 2 0.470 4 −4.36 1.57 −1.75 2.13 −2.51 1.89 1.04 −1.93 2.06 1.00 0.47 −1.76 1.43
VI 2 1 7.44 −4.50 1.56 −1.83 2.10 −2.71 1.85 1.02 −1.95 2.06 2.25 0.44 −1.73 1.40


Conflicts of interest

There are no conflicts to declare.

Acknowledgements

The authors acknowledge support by the state of Baden-Württemberg through bwHPC and thank Dr Fajun Zhang for insightful discussions on this manuscript and for providing experimental structure factors.

Notes and references

  1. P.-H. Peng, K.-W. Hsu and K.-J. Wu, Am. J. Cancer Res., 2021, 11, 3766–3776 CAS.
  2. A. Pande, J. Pande, N. Asherie, A. Lomakin, O. Ogun, J. King and G. B. Benedek, Proc. Natl. Acad. Sci. U. S. A., 2001, 98, 6116–6120 CrossRef CAS PubMed.
  3. X. Zhang, Y. Lin, N. A. Eschmann, H. Zhou, J. N. Rauch, I. Hernandez, E. Guzman, K. S. Kosik and S. Han, PLoS Biol., 2017, 15, 1–28 Search PubMed.
  4. P. R. ten Wolde and D. Frenkel, Science, 1997, 277, 1975–1978 CrossRef CAS PubMed.
  5. F. Zhang, G. Zocher, A. Sauter, T. Stehle and F. Schreiber, J. Appl. Crystallogr., 2011, 44, 755–762 CrossRef CAS.
  6. A. Stradner and P. Schurtenberger, Soft Matter, 2020, 16, 307–323 RSC.
  7. G. A. Vliegenthart and H. N. W. Lekkerkerker, J. Chem. Phys., 2000, 112, 5364–5369 CrossRef CAS.
  8. M. G. Noro and D. Frenkel, J. Chem. Phys., 2000, 113, 2941–2944 CrossRef CAS.
  9. F. Platten, N. E. Valadez-Pérez, R. Castañeda-Priego and S. U. Egelhaaf, J. Chem. Phys., 2015, 142, 174905 CrossRef PubMed.
  10. G. Foffi and F. Sciortino, J. Phys. Chem. B, 2007, 111, 9702–9705 CrossRef CAS PubMed.
  11. H. Liu, S. K. Kumar, F. Sciortino and G. T. Evans, J. Chem. Phys., 2009, 130, 044902 CrossRef PubMed.
  12. M. A. Blanco, mAbs, 2022, 14, 2044744 CrossRef PubMed.
  13. J. Jin, A. J. Pak, A. E. P. Durumeric, T. D. Loose and G. A. Voth, J. Chem. Theory Comput., 2022, 18, 5759–5791 CrossRef CAS PubMed.
  14. S. Kmiecik, D. Gront, M. Kolinski, L. Wieteska, A. E. Dawid and A. Kolinski, Chem. Rev., 2016, 116, 7898–7936 CrossRef CAS PubMed.
  15. S. J. Marrink, H. J. Risselada, S. Yefimov, D. P. Tieleman and A. H. de Vries, J. Phys. Chem. B, 2007, 111, 7812–7824 CrossRef CAS PubMed.
  16. E. Fagerberg and M. Skepö, J. Chem. Inf. Model., 2023, 63, 4079–4087 CrossRef CAS PubMed.
  17. L. Monticelli, J. Chem. Theory Comput., 2012, 8, 1370–1378 CrossRef CAS PubMed.
  18. A. Chaudhri, I. E. Zarraga, S. Yadav, T. W. Patapoff, S. J. Shire and G. A. Voth, J. Phys. Chem. B, 2013, 117, 1269–1279 CrossRef CAS PubMed.
  19. G. Wang, Z. Varga, J. Hofmann, I. E. Zarraga and J. W. Swan, J. Phys. Chem. B, 2018, 122, 2867–2880 CrossRef CAS PubMed.
  20. S. Izadi, T. W. Patapoff and B. T. Walters, Biophys. J., 2020, 118, 2741–2754 CrossRef CAS PubMed.
  21. H. Shahfar, J. K. Forder and C. J. Roberts, J. Phys. Chem. B, 2021, 125, 3574–3588 CrossRef CAS PubMed.
  22. A. Grünberger, P.-K. Lai, M. A. Blanco and C. J. Roberts, J. Phys. Chem. B, 2013, 117, 763–770 CrossRef PubMed.
  23. S. Pusara, P. Yamin, W. Wenzel, M. Krstić and M. Kozlowska, Phys. Chem. Chem. Phys., 2021, 23, 12780–12794 RSC.
  24. F. Zhang, M. W. A. Skoda, R. M. J. Jacobs, S. Zorn, R. A. Martin, C. M. Martin, G. F. Clark, S. Weggler, A. Hildebrandt, O. Kohlbacher and F. Schreiber, Phys. Rev. Lett., 2008, 101, 148101 CrossRef CAS PubMed.
  25. F. Roosen-Runge, F. Zhang, F. Schreiber and R. Roth, Sci. Rep., 2014, 4, 7016 CrossRef CAS PubMed.
  26. M. R. Fries, D. Stopper, M. W. A. Skoda, M. Blum, C. Kertzscher, A. Hinderhofer, F. Zhang, R. M. J. Jacobs, R. Roth and F. Schreiber, Sci. Rep., 2020, 10, 10349 CrossRef CAS PubMed.
  27. A. K. Gaigalas, J. B. Hubbard, M. McCurley and S. Woo, J. Phys. Chem., 1992, 96, 2355–2359 CrossRef CAS.
  28. F. Zhang, M. W. A. Skoda, R. M. J. Jacobs, R. A. Martin, C. M. Martin and F. Schreiber, J. Phys. Chem. B, 2007, 111, 251–259 CrossRef CAS PubMed.
  29. A. Bujacz, Acta Crystallogr., Sect. D: Biol. Crystallogr., 2012, 68, 1278–1289 CrossRef CAS PubMed.
  30. F. Ameseder, R. Biehl, O. Holderer, D. Richter and A. M. Stadler, Phys. Chem. Chem. Phys., 2019, 21, 18477–18485 RSC.
  31. C. Beck, M. Grimaldo, M. K. Braun, L. Bühl, O. Matsarskaia, N. H. Jalarvo, F. Zhang, F. Roosen-Runge, F. Schreiber and T. Seydel, Soft Matter, 2021, 17, 8506–8516 RSC.
  32. F. Hirschmann, H. Lopez, F. Roosen-Runge, T. Seydel, F. Schreiber and M. Oettel, J. Chem. Phys., 2023, 158, 084112 CrossRef CAS PubMed.
  33. M. J. Abraham, T. Murtola, R. Schulz, S. Páll, J. C. Smith, B. Hess and E. Lindahl, SoftwareX, 2015, 1–2, 19–25 CrossRef.
  34. W. Humphrey, A. Dalke and K. Schulten, J. Mol. Graphics, 1996, 14, 33–38 CrossRef CAS PubMed.
  35. H. Hamaker, Physica, 1937, 4, 1058–1072 CrossRef CAS.
  36. E. Jurrus, D. Engel, K. Star, K. Monson, J. Brandi, L. E. Felberg, D. H. Brookes, L. Wilson, J. Chen, K. Liles, M. Chun, P. Li, D. W. Gohara, T. Dolinsky, R. Konecny, D. R. Koes, J. E. Nielsen, T. Head-Gordon, W. Geng, R. Krasny, G.-W. Wei, M. J. Holst, J. A. McCammon and N. A. Baker, Protein Sci., 2018, 27, 112–128 CrossRef CAS PubMed.
  37. C. R. Søndergaard, M. H. M. Olsson, M. Rostkowski and J. H. Jensen, J. Chem. Theory Comput., 2011, 7, 2284–2295 CrossRef PubMed.
  38. A. Salis, M. Boström, L. Medda, F. Cugia, B. Barse, D. F. Parsons, B. W. Ninham and M. Monduzzi, Langmuir, 2011, 27, 11597–11604 CrossRef CAS PubMed.
  39. F. Zhang, S. Weggler, M. J. Ziller, L. Ianeselli, B. S. Heck, A. Hildebrandt, O. Kohlbacher, M. W. A. Skoda, R. M. J. Jacobs and F. Schreiber, Proteins: Struct., Funct., Bioinf., 2010, 78, 3450–3457 CrossRef CAS PubMed.
  40. O. Matsarskaia, F. Roosen-Runge, G. Lotze, J. Möller, A. Mariani, F. Zhang and F. Schreiber, Phys. Chem. Chem. Phys., 2018, 20, 27214–27225 RSC.
  41. F. Zhang, personal communication, 2024.
  42. N. Gnan, F. Sciortino and E. Zaccarelli, in Patchy Particle Models to Understand Protein Phase Behavior, ed. J. J. McManus, Springer New York, New York, NY, 2019, pp. 187–208 Search PubMed.
  43. K. Shoemake, SIGGRAPH Comput. Graph., 1985, 19, 245–254 CrossRef.
  44. A. M. Ferrenberg and R. H. Swendsen, Phys. Rev. Lett., 1988, 61, 2635–2638 CrossRef CAS PubMed.
  45. A. M. Ferrenberg and R. H. Swendsen, Phys. Rev. Lett., 1989, 63, 1195–1198 CrossRef CAS PubMed.
  46. R. H. Swendsen, Phys. A, 1993, 194, 53–62 CrossRef.
  47. N. B. Wilding, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1995, 52, 602–611 CrossRef CAS PubMed.
  48. P. Virnau and M. Müller, J. Chem. Phys., 2004, 120, 10925–10930 CrossRef CAS PubMed.
  49. M. A. Blanco and V. K. Shen, J. Chem. Phys., 2016, 145, 155102 CrossRef PubMed.
  50. J. K. G. Dhont, An Introduction to Dynamics of Colloids, Elsevier, 1996 Search PubMed.
  51. J. A. Barker and D. Henderson, J. Chem. Phys., 1967, 47, 4714–4721 CrossRef CAS.
  52. D. Gazzillo and D. Pini, J. Chem. Phys., 2013, 139, 164501 CrossRef PubMed.
  53. N. Kern and D. Frenkel, J. Chem. Phys., 2003, 118, 9882–9889 CrossRef CAS.
  54. E. Bianchi, G. Kahl and C. N. Likos, Soft Matter, 2011, 7, 8313–8323 RSC.
  55. S. Kohara, Y. Matsuzawa and Y. Kuroda, Chem. Phys. Lett., 2022, 802, 139767 CrossRef CAS.
  56. A. Lošdorfer Božič, Soft Matter, 2018, 14, 1149–1161 RSC.
  57. D. Notarmuzi and E. Bianchi, Soft Matter, 2024, 20, 7601–7614 RSC.
  58. M. K. Braun, M. Wolf, O. Matsarskaia, S. Da Vela, F. Roosen-Runge, M. Sztucki, R. Roth, F. Zhang and F. Schreiber, J. Phys. Chem. B, 2017, 121, 1731–1739 CrossRef CAS PubMed.
  59. F. Surfaro, R. Maier, K.-F. Pastryk, F. Zhang, F. Schreiber and R. Roth, J. Chem. Phys., 2023, 158, 164902 CrossRef CAS PubMed.
  60. M. K. Braun, A. Sauter, O. Matsarskaia, M. Wolf, F. Roosen-Runge, M. Sztucki, R. Roth, F. Zhang and F. Schreiber, J. Phys. Chem. B, 2018, 122, 11978–11985 CrossRef CAS PubMed.

This journal is © The Royal Society of Chemistry 2024
Click here to see how this site uses Cookies. View our privacy policy here.