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
First published on 8th October 2024
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 is frequently invoked which is based on the behavior of short-range attractive fluids (Noro–Frenkel rule, ). However, if anisotropic models for the protein–protein interaction are considered, e.g. the Kern–Frenkel (KF) patchy particle model, the value of the 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 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 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.
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. for three patches, for four patches etc.10 Thus one would expect also a deviation of 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 vs. patch number.11 This implies that using universal 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 -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 ). Furthermore the relationship of our results to experiments is discussed, before we conclude in Section 6.
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) |
Volume exclusions are modelled by the repulsive part of the Lennard-Jones potential Urep, which is cut-off at its minimum:
(2) |
Dispersion interactions are modelled by the van der Waals potential UvdW for two equally sized spheres with:35
(3) |
Electrostatic interactions are approximated via a screened Yukawa potential Uel:
(4) |
Z → CMZ, | (5) |
The last term Ucont in eqn (1) compensates the usage of cutoffs by making the total potential continuous. It reads:
(6) |
For illustration, the complete interaction potential specified by eqn (1) is depicted in Fig. 2 for some typical parameters.
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.
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.
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 = N − sE) 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.
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 = N − sE 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.
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).
(7) |
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.
(8) |
For obtaining reduced second virial coefficients 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 .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.
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. |
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 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%.
Fig. 7 Reduced second virial coefficient 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 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.
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.
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) |
(10) |
Fig. 9 (a) Reduced second virial coefficients = 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 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 for all critical BSA systems together with regions in which 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 in Fig. 6(b) and 9(a) (respectively) and also in Table 2. For all investigated isotropic systems, which is different from . 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: 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 , as well as 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 are around 0.34…0.47 with a moderate spread and all values of 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 (similar to ).
Comparing this to the KF results, these values correspond to 4 patches (from ) or 3 patches for the CM = 2 systems (from ), which points to a considerable effective anisotropy in the BSA model.
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 is reduced in magnitude and they all exhibit a stronger charge dependence than before: for method (ii) (with σeff = σrepeff) we find for the systems with CM = 0 and 1, and 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: (CM = 0, 1) and ≈−2.5…−2.7 (CM = 2). Method (iv) yields approximately values that fall in the middle of (ii) and (iii): (CM = 0, 1) and (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 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
(11) |
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
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. |
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 , which is of the order of 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, 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 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 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 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 . 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 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.
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 points to an effective patchiness of 4 Kern–Frenkel patches 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 -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 (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 -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 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.
System | CM | Anisotropic systems | Isotropic systems | |||||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
T c | A H,c | (i): σNFeff | (ii): σrepeff | (iii): σdryeff | (iv): σcpeff | T isoc | (i): σNF,eff | |||||||||
σ NFeff | σ repeff | σ dryeff | σ beadeff | σ cpeff | σ 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 |
This journal is © The Royal Society of Chemistry 2024 |