Using schematic models to understand the microscopic basis for inverted solubility in $\gamma$D-crystallin

TL;DR

This study investigates the microscopic mechanisms behind inverted solubility in human γD-crystallin, finding that hydrophobicity changes are unlikely the cause and that temperature-dependent interactions require fine-tuning to produce solubility inversion.

Contribution

The paper combines experimental surface hydrophobicity measurements with schematic modeling to assess the role of hydrophobic effects and temperature-dependent interactions in inverted protein solubility.

Findings

Hydrophobicity does not significantly increase with mutation.

Hydrophobic effects alone may not explain solubility inversion.

Temperature-dependent interactions have minimal impact on critical properties.

Abstract

Inverted solubility--a crystal melting upon cooling--is observed in a handful of proteins, such as carbomonoxy hemoglobin and $\gamma$D-crystallin. In human $\gamma$D-crystallin, the phenomenon is associated with the mutation of the 23$^\mathrm{rd}$ residue, a proline, to a threonine, serine or valine. One proposed microscopic mechanism for this effect entails an increase in hydrophobicity upon mutagenesis. Recent crystal structures of a double mutant that includes the P23T mutation allows for a more careful investigation of this proposal. Here, we first measure the surface hydrophobicity of various mutant structures of this protein and determine that it does not discernibly increase upon the mutating the 23$^\mathrm{rd}$ residue. We then investigate the solubility inversion regime with a schematic patchy particle model that includes one of three models for temperature-dependent patch energies: two of the hydrophobic effect, and a more generic description. We conclude that while solubility inversion due to the hydrophobic effect may be possible, microscopic evidence to support it in $\gamma$D-crystallin is weak. More generally, we find that solubility inversion requires a fine balance between patch strengths and the temperature-dependent contribution, which may explain why inverted solubility is not commonly observed in proteins. In any event, we also find that the temperature-dependent interaction has only a negligible impact on the critical properties of the $\gamma$D-crystallin, in line with previous experimental observations.

Tables

Table 1. Table 1: Patch parameters obtained from all-atom molecular dynamics simulations and slightly modified to adjust the crystal geometry

Contact	$\cos {\delta }_{\alpha }$	$\cos {\delta }_{\beta }$	$\Delta {\varphi }_{\alpha \beta }$ (rad)	$\epsilon$ ($k_{B}T$)	$\lambda$ ($\sigma$)
c1	0.956	0.963	0.07	21.5	1.025
c2	0.935	0.935	0.21	6.3	1.1
c3	0.992	0.994	0.14	10.3	1.06
c4	0.997	0.997	0.06	15.8	1.025
c5	0.947	0.966	0.13	6.0	1.037
cI	0.972	0.981	0.21	6.3	1.1
cII	0.96	0.96	0.3	8.7	1.11
cIII	0.95	0.95	0.2	10.1	1.1
cIV	0.96	0.937	0.3	4.3	1.14
cV	0.95	0.95	0.21	18.5	1.153
cVI	0.898	0.927	0.3	8.6	1.11
cVII	0.95	0.95	0.21	18.5	1.53

Table 2. Table 2: Vectors defining each patch location and dihedral angle, ϕ α ​ β subscript italic-ϕ 𝛼 𝛽 \phi_{\alpha\beta} . Note that because of geometric constraints of the DBN structure, contacts cII, cIV, and cVI are collocated.

Contact	${\hat{e}}_1$	${\hat{e}}_2$	${\hat{e}}_3$	${\varphi }_{\alpha \beta }$
c1	0.6614	0.6543	0.3667	-2.2703
	0.6473	-0.2853	0.7068
c2/cI	-0.9854	0.0471	0.1101	0.0
	0.9854	-0.0471	-0.1101
c3	-0.0141	-0.9402	0.3403	0.0
	0.0141	0.9402	-0.3403
c4	0.3231	0.1911	-0.9269	-2.1288
	0.3089	-0.7486	-0.5867
c5	-0.6474	0.2851	-0.7068	2.2698
	-0.6613	-0.6543	-0.3669
cII*	-0.9494	0.3101	-0.0499	-2.4946
	0.0011	0.4334	0.9012
cIII	-0.6455	-0.5123	-0.5665	-1.80
	-0.7102	-0.5524	0.4365
cIV*	-0.9494	0.3101	-0.0499	-2.4946
	0.0011	0.4334	0.9012
cV/VII	0.3162	-0.8026	-0.5058	-2.1529
	0.4927	0.8531	0.1716
cVI*	-0.9494	0.3101	-0.0499	-2.4946
	0.0011	0.4334	0.9012

Equations

$$H_{\zeta }=\frac{\sum _{i\in S}{f}_{\zeta }(i)A(i)}{\sum _{i\in S}A(i)},$$

$$u(r_{ij},{\Omega }_i,{\Omega }_j)=u_{\mathrm{HS}}(r_{ij})+\sum _{a,b}^n{u}_{ab}(r_{ij},{\Omega }_i,{\Omega }_j),$$

$$u_{ab}=v_{ab}(r_{ij})f_{ab}({\Omega }_i,{\Omega }_j).$$

$$\displaystyle v_{ab}(r_{ij})=\left\{\begin{array}[]{ll}-\varepsilon_{ab}(T),&\sigma<r_{ij}<\lambda_{a}+\lambda_{b}\\ 0,&\text{otherwise}\end{array}\right.,$$

$$\displaystyle f_{ab}=\left\{\begin{array}[]{ll}1,&\theta_{a,ij}\leq\delta_{a}\;\mathrm{and}\;\theta_{b,ij}\leq\delta_{b}\\ 0,&\text{otherwise}\end{array}\right.\times\left\{\begin{array}[]{ll}1,&\psi_{ij}\in[\varphi_{ab}-\Delta\varphi_{ab},\varphi_{ab}+\Delta\varphi_{ab}]\\ 0,&\text{otherwise}\end{array}\right.$$

$$E_s=\frac{E_{\mathrm{os}}+E_{\mathrm{ds}}{e}^{-\beta (E_{\mathrm{ds}}-E_{\mathrm{os}})}}{1+e^{-\beta (E_{\mathrm{ds}}-E_{\mathrm{os}})}}$$

$$E_b=\frac{E_{\mathrm{ob}}+E_{\mathrm{db}}{e}^{-\beta (E_{\mathrm{db}}-E_{\mathrm{ob}})}}{1+e^{-\beta (E_{\mathrm{db}}-E_{\mathrm{ob}})}}$$

$$\displaystyle S_{s}/k_{B}=\log\Big{(}\frac{q_{\mathrm{os}}+q_{\mathrm{ds}}e^{-\beta(E_{\mathrm{ds}}-E_{\mathrm{os}})}}{1+e^{-\beta(E_{\mathrm{ds}}-E_{\mathrm{os}})}}\Big{)}$$

$$\displaystyle S_{b}/k_{B}=\log\Big{(}\frac{q_{\mathrm{ob}}+q_{\mathrm{db}}e^{-\beta(E_{\mathrm{db}}-E_{\mathrm{ob}})}}{1+e^{-\beta(E_{\mathrm{db}}-E_{\mathrm{ob}})}}\Big{)}.$$

$${\epsilon }_4^{\prime }={\epsilon }_4+n_w\Delta \epsilon (\beta ),$$

$$-{\epsilon }_4^{\prime }=-{\epsilon }_4-2{\epsilon }_w+\frac{2}{\beta }\Delta s_w,$$

$$\displaystyle\varepsilon_{4}^{\prime}(T)=\frac{\varepsilon_{4}}{2}\Bigg{[}1+\tanh\Bigg{(}\frac{T-T_{a}}{\tau}\Bigg{)}\Bigg{]},$$

$$\frac{\beta p}{\rho }=1+B_2\rho ,$$

$$\beta {\mu }_f=\beta {\mu }^{\mathrm{id}}+2B_2\rho =\log {\Lambda }^3\rho +2B_2\rho ,$$

$$\rho =\frac{-1+\sqrt{1+4B_2\beta p}}{2B_2}.$$

$$\beta {\mu }_c(\beta ,p)={\beta }_0{\mu }_c({\beta }_0,p)+\underset{{\beta }_0}{\overset{\beta }{\int }}\frac{⟨H({\beta }^{\prime })⟩}{N}d{\beta }^{\prime }+\underset{{\beta }_0}{\overset{\beta }{\int }}{\beta }^{\prime }\frac{⟨dU/d{\beta }^{\prime }⟩}{N}d{\beta }^{\prime },$$

$$\underset{{\beta }_0}{\overset{\beta }{\int }}\frac{⟨H({\beta }^{\prime })⟩}{N}d{\beta }^{\prime }\approx \frac{1}{N}\underset{{\beta }_0}{\overset{\beta }{\int }}\frac{⟨U({\beta }^{\prime })⟩}{N}d{\beta }^{\prime }+\frac{p}{\rho }(\beta -{\beta }_0).$$

$$\beta {\mu }_c(\beta ,p)\approx {\beta }_0{\mu }_c({\beta }_0,p)+\frac{p}{\rho }(\beta -{\beta }_0)+\xi (\beta ),$$

$$\displaystyle\xi(\beta)=\int\limits_{\beta_{0}}^{\beta}\Big{[}-\varepsilon_{\mathrm{tot}}+\Delta\varepsilon(1)n_{w}-2n_{w}\varepsilon_{w}+2n_{w}(\frac{d\Delta s_{w}}{d\beta^{\prime}}-\beta^{\prime}\frac{d\varepsilon_{w}}{d\beta^{\prime}})\Big{]}d\beta^{\prime}$$

$$\displaystyle=(\Delta\varepsilon(1)n_{w}-\varepsilon_{\mathrm{tot}})(\beta-\beta_{0})+2n_{w}\Big{[}(\Delta s_{w}(\beta)-\Delta s_{w}(\beta_{0}))-(\beta\varepsilon_{w}(\beta)-\beta_{0}\varepsilon_{w}(\beta_{0}))\Big{]}.$$

$$\displaystyle\frac{\beta\mu_{c}}{d\beta}\approx(\Delta\varepsilon(1)n_{w}-\varepsilon_{\mathrm{tot}})+2n_{w}\Big{[}\frac{d\Delta s_{w}}{d\beta}-\varepsilon_{w}-\beta\frac{d\varepsilon_{w}}{d\beta}\Big{]},$$

$$\Gamma (\beta )\equiv \frac{d\Delta s_w}{d\beta }-{\epsilon }_w-\beta \frac{d{\epsilon }_w}{d\beta }=\frac{{\epsilon }_{\mathrm{tot}}}{2n_w}-\frac{\Delta \epsilon (1)}{2}.$$

$$\displaystyle U(\beta)/N=-\varepsilon_{\mathrm{tot}}-2\Delta s_{w}\Big{(}\frac{1}{\beta_{\mathrm{ref}}}-\frac{1}{\beta}\Big{)}.$$

$$\displaystyle\beta\mu_{c}(\beta,p)=\beta\mu_{c}(\beta_{0},p)+\Big{(}-\varepsilon_{\mathrm{tot}}-\frac{2\Delta s_{w}}{\beta_{\mathrm{ref}}}-\frac{p}{\rho}\Big{)}\Big{(}\beta-\beta_{0}\Big{)}.$$

$$\beta A^{\prime }=\beta A-\beta U_0/N+\beta U_0^{\prime }/N,$$

$$n_w=A_4{\rho }_w\underset{3\text{A˚}}{\overset{4.5\text{A˚}}{\int }}{g}_C(r)dr,$$

$$\displaystyle B_{2,\alpha\beta}=-\frac{1}{2}4\pi\int\Big{(}e^{-\beta U(r)}-1\Big{)}r^{2}dr,$$

$$\displaystyle\lambda_{\alpha\beta}=\Bigg{(}\frac{3\int(e^{-\beta U_{\mathrm{PMF}}(r)}-1)r^{2}dr}{e^{\beta\varepsilon_{\alpha\beta}}-1}+1\Bigg{)}.$$

$$B_2=-\frac{1}{2V}\int d\mathbf{r}\int d\chi \int d{\mathbf{\Omega }}_1\int d{\mathbf{\Omega }}_2(e^{-\beta u(|\mathbf{r}|,{\mathbf{\Omega }}_1,{\mathbf{\Omega }}_2)}-1),$$

$$\displaystyle\frac{B_{2}}{B_{\mathrm{HS}}}=1-\sum(\lambda_{ij}^{3}-1)\sin^{2}\Big{(}\frac{\delta_{i}}{2}\Big{)}\sin^{2}\Big{(}\frac{\delta_{j}}{2}\Big{)}(e^{\beta\varepsilon_{ij}}-1)(2\Delta\phi_{ij})$$

$$\displaystyle-(\lambda^{3}_{IV}-\lambda^{3}_{II})(e^{\beta\varepsilon_{IV}}-1)\sin^{2}\Big{(}\frac{\delta_{IV,\alpha}}{2}\Big{)}\sin^{2}\Big{(}\frac{\delta_{IV,\beta}}{2}\Big{)}(2\Delta\chi)$$

$$\displaystyle-(\lambda_{II}^{3}-1)(e^{\beta(\varepsilon_{II}+\varepsilon_{IV}+\varepsilon_{VI})}-1)\sin^{2}\Big{(}\frac{\delta_{II,\alpha}}{2}\Big{)}\sin^{2}\Big{(}\frac{\delta_{II,\beta}}{2}\Big{)}(2\Delta\chi)$$

$$\displaystyle-(\lambda_{VI}^{3}-1)(e^{\beta(\varepsilon_{IV}+\varepsilon_{VI})})\sin^{2}\Big{(}\frac{\delta_{IV,\alpha}}{2}\Big{)}\sin\Big{(}\frac{\delta_{IV,\beta}+\delta_{IV,\alpha}}{2}\Big{)}\sin\Big{(}\frac{\delta_{IV,\beta}-\delta_{IV,\alpha}}{2}\Big{)}(2\Delta\chi)$$

Full text

1 Introduction

Proteins can self-organize into a rich variety of superstructures 1, such as crystals 2, virus capsids 3, disease-forming aggregates 4, and biomaterials 5. A key challenge is to understand how microscopic features of solvated proteins can give rise to such complex phase diagrams, and eventually to design systems that reliably assemble 6, 7, 8, 5, 9, 10. In this context, coarse-grained models are especially valuable, because they help both pinpoint and abstract the microscopic features that can reproduce the experimentally observed behavior. (Since simulating protein self-assembly typically requires hundreds to thousands of protein copies, which are themselves comprised of thousands of atoms, such models are also a computational necessity 11, 12, 13.) For example, even simple models of short-ranged 14, 15, anisotropic pair interactions can largely explain the phase behavior of globular proteins 16, 17, 1. Understanding the assembly of some systems, however, requires coarse-grained models with additional features, such as shape anisotropy for viral capsid and amyloid forming proteins 3, 1. Capturing certain features of protein crystallization, which is key to protein structure determination by diffraction methods 18, also requires enhanced patchy particle models 2.

Systems that exhibit atypical phase behaviors are essential test of our understanding of the physico-chemical processes that underlie protein assembly. One such phenomenon is inverted crystal solubility, i.e., the decrease of solubility with increasing temperature. This phenomenon is observed in a handful of proteins, such as some single mutants of $\gamma$D-crystallin 19, 20, and the wild type carbomonoxy-hemoglobin C 21. (The temperature invariant solubility of apoferritin is a limit case 22.) Thermodynamically, inverted solubility suggests that as temperature increases, the Gibbs free energy of crystallization decreases, and hence that the crystal becomes increasingly more stable than the fluid. The phenomenon is often attributed to a large and positive entropy gain upon crystallization. Crystal formation is then possible despite the enthalpy of crystallization being non-negative 21, 23, 24. Because the solute contribution to the change in entropy is typically negative, the solvent contribution is considered to be the key microscopic determinant 24, 23, 21.

The association of inverted solubility in proteins with the hydrophobic effect comes from our understanding of the aqueous solvation of hydrocarbons, which presents an analogous solubility regime 25. A minimal model for this effect was proposed by Lee and Graziano 26, who refined Muller’s two state description of water as having either active or broken hydrogen bonds 27. The final model formulation considers water as being in one of four states: disordered shell (ds), ordered shell (os), disordered bulk (db), and ordered bulk (ob). Shiryayev et al. then used this model to estimate the phase diagram of model globular proteins with isotropic interactions assumed to be driven exclusively by hydrophobic interactions 23. Although the resulting phase behavior does present an inverted solubility regime, it is unclear whether it persists for more realistic protein models, with directional interactions and a complex surface mosaic of hydrophilic and hydrophobic interactions with the solvent. In other words, while the hydrophobic scenario for solubility inversion in proteins is thermodynamically sound, microscopic evidence for it remains limited. The generality of the underlying physical arguments is also seemingly incompatible with the relatively rare occurrence of inverted solubility in experiments.

Here, we examine this microscopic scenario in the context of a double mutant (R36S+P23T) of the human $\gamma$D-crystallin, which has been shown to form two competing crystals: a normal solubility structure (DBN, PDB 28 ID: 6ETA) and an inverted solubility structure (DBI, PDB ID: 6ETC) 29, 30. The solubility inversion is here most likely associated with the mutation in the 23rd residue because the single mutants P23T, P23S, and P23V also exhibit inverted solubility 31. In earlier work, we have parameterized a patchy model for this system and have obtained a solubility inversion regime by completely deactivating the patch containing the 23rd residue at low temperatures 30. Interestingly, the DBI crystal does not present any obvious structural feature, other than the formation of a hydrogen bond through the 23rd residue. Here, we test three different temperature-dependent interaction potentials: the generic model we previously considered, and two that explicitly model the hydrophobic scenario. We use these models to test the hydrophobic scenario as well as the robustness of the inverted solubility regime with respect to model parameters. We thus attempt to elucidate why inverted solubility is not more commonly observed. We further explore the relationship between the liquid-liquid critical point and the solubility curve which has been studied experimentally for some of these systems 31. The plan for the rest of this paper is as follows. We first consider whether an increase in surface hydrophobicity can be discerned upon introducing the solubility inverting mutations (Sec. 2). We then introduce a patchy protein model for these proteins (Sec. 3.1) along with the different temperature-dependent patch models (Sec. 3.2), and the methods used to determine solubility lines (Sec. 3.3). Sections 4 and 5 provide a detailed analysis of these patchy models, and we conclude with proposals for further discerning experiments in Section 6.

2 Surface Hydrophobicity

As a first consideration of the reasonableness of the hydrophobicity scenario, we first evaluate the surface hydrophobicity of various human $\gamma$D-crystallin crystal structures. Were the P23T mutation to consistently increase surface hydrophobicity, it would serve as strong evidence for the decrease in protein solubility upon mutagenesis to be driven by the hydrophobic effect. By studying the relative binding propensity of two dyes known to bind hydrophobic surfaces, Pande et al. have inferred that P23T, P23S, and P23V mutants of human $\gamma$D-crystallin do present a higher surface hydrophobicity than the wild type (WT) protein. 20 In order to test the robustness of this interpretation, we here consider different scales that quantify hydrophobicity at the amino acid level. More specifically, we compute an average of hydrophobicity indices of solvent-exposed residues 32 weighted by their solvent accessible surface area (SASA) 33, for five different scales: grand average of hydropathy (GRAVY) 34, as well as the scales of Wimley and White (ww) 35, Hessa et al. (hh) 36, Moon and Fleming (mf) 37, and Zhao and London (also known as transmembrane tendency, tt) 38. Each of these scales assigns a hydrophobicity index to each residue type; all but hh and mf assign positive values to hydrophobic residues.

We compute hydrophobicity for three sets, $S$, of amino acids: (i) the entire protein surface, (ii) the surface of its N-terminus, i.e., the first 82 residues (including the solubility inverting 23rd residue), and (iii) the surface residues in the DBI contact that includes the 23rd residue (Patch 4 as per Sec. 3) 30. The hydrophobicity, $H_{\zeta}$, for a given scale $\zeta$ is then obtained as

[TABLE]

where $f_{\zeta}(i)$ is the hydrophobicity index for residue $i$, and $A(i)$ is its SASA. We specifically consider: WT (PDB ID: 1HK0 39), the P23T single mutant (PDB ID: 4JGF 40), the R36S single mutant (PDB ID: 2G98 41), the R58H single mutant (PDB ID: 1H4A 39), DBI (PDB ID: 6ETC 30), and DBN (PDB ID: 6ETA 30). Of these, only WT, R36S, and R58H do not have a mutation at the 23rd residue. Note that missing residues are completed using Modeller 42 within Chimera 43, and all crystal water molecules are removed. In order to estimate the error on these measured hydrophobicities, 100 structures per crystal are created by perturbing each particle coordinate by a random number selected from a Gaussian distribution with standard deviation corresponding to the coordinate error specified in the PDB file. Two assumptions are made in estimating these error bars. First, the coordinate error reported in the PDB entry is assumed to be distributed uniformly and isotropically across all protein atoms. In reality, certain domains or residues in proteins are more mobile and thus harder to resolve by X-ray diffraction than others, but residue-level information is not available. This assumption thus overestimates the error in more localized parts of the protein and underestimates the error in more mobile parts. Second, the refined structures do not precisely capture the actual protein structure, as suggested by $R_{\mathrm{free}}$ values ranging from 0.174 to as high as 0.284, hence possibly creating artificial hydrophobicity differences between different mutants, or, conversely, underestimating them.

The resulting hydrophobicity estimates are shown in Fig. 1. We first compare the DBN and DBI structures, which are obtained from the same double mutant, R36S+P23T, and which are structurally very similar 30. As expected, nearly all measurements for DBN and DBI overlap within their $95\%$ confidence intervals. The only exception is the hydrophobicity of Patch 4 measured by the mf scale. This could be because mf uniquely classifies prolines as hydrophobic. This discrepancy could then amplify the minute difference in surface exposure of Patch 4 prolines between DBI and DBN.

Overall, the N-terminus is the most hydrophobic region in nearly all scales and for all structures. However, other observations are not consistent across scales. A number of nonmonotonicities can indeed be observed across different hydrophobicity scales. For instance, Patch 4 is more hydrophobic in DBI than in WT in the GRAVY, hh, and mf scales, but the ww and tt scales present no discernible difference. Similarly, Patch 4 is more hydrophobic in R36S than WT for the mf scale, but the reverse is true for hh. These discrepancies reflect the different ordering of residues in the various hydrophobicity scales. For instance, GRAVY, which is calculated from experimental measurements of transfer free energies from water to water vapor, tends to assign aromatic side chains lower hydrophobicities than the other four scales, which instead consider the tendency of residues to transfer from water to within a lipid bilayer, a measurement prone to more experimental uncertainty 35.

Interestingly, the N-terminus of the P23T mutant is the least hydrophobic structure for the GRAVY and mf scales. This trend, however, disappears when only Patch 4 residues are considered. Patch 4, which controls solubility inversion, is actually less hydrophobic than the overall N-terminus or the entire protein, except on the mf scale. Only for this last scale is Patch 4 clearly more hydrophobic. A similar inconsistency exists for Patch 4 of DBI, which is more hydrophobic than the other proteins for GRAVY and mf, but for these two scales P23T and DBN are not discernibly more hydrophobic than the structures without the mutation in the 23rd residue.

In summary, in none of the hydrophobicity scales do the structures with the (solubility-inverting) P23T mutation have a statistically and consistently higher hydrophobicity than those without. P23T mutations even result in lower hydrophobicity estimates on some scales. While our measurements are subject to errors from the crystal structure accuracy, as well as the imperfections of the hydrophobicity scales themselves, a microscopic change to the protein surface that could underlie inversion of solubility remains elusive from this viewpoint. The question that then arises is whether other measurements could be more revealing. Because hydrophobicity scales are but an indirect measure of protein-water interactions (and thus protein-protein interactions), more direct measurement could be microscopically more revealing. In that context, we next consider the hydrophobicity scenario from a thermodynamic perspective.

3 Solubility Lines from Patchy Models

Because a clear enhancement of hydrophobicity cannot be detected directly in mutants with inverted solubility, we next consider the thermodynamics of patchy models that incorporate various temperature-dependent patch energies. A schematic model of the double mutant of human $\gamma$D-crystallin was previously studied in Ref. 30, but it is here modified to consider the hydrophobic scenario and then perturbed to evaluate the robustness of its inverted solubility regime.

3.1 Patchy Model

The schematic model consists of hard particles with attractive patches

[TABLE]

where $r_{ij}$ is the distance between particles $i$ and $j$, $\Omega$ denotes the particle orientation, and $u_{\mathrm{HS}}(r_{ij})$ is the hard sphere potential for particles of diameter $\sigma$. The sum runs over all patch pairs, with $n$ the total number of patches. The second contribution, $u_{ab}$, is further broken down into radial and angular parts

[TABLE]

The radial part, $v_{ab}$, is a square-well interaction

[TABLE]

with interaction ranges $\lambda_{a}$ and $\lambda_{b}$ of patches $a$ and $b$, respectively, and with either constant or temperature-dependent patch energy $-\varepsilon_{ab}(T)$. The orientational part

[TABLE]

contains two contributions. The first ensures that the relative particle orientation enables them to interact with $\delta_{a}$ and $\delta_{b}$ the angular width for patches $a$ and $b$, respectively (Fig. 2a). The second limits the range $\varphi_{ab}\pm\Delta\varphi_{ab}$ of dihedral angles $\psi_{ij}$ allowed for each pair (Fig. 2b), with $\theta_{a,ij}$ the angle between the vector defining the location of patch $a$ and the vector that connects the centers of particles $i$ and $j$, and $\theta_{b,ij}$ similarly for patch $b$.

This model is parameterized such that each patch corresponds to a crystal contact in either the DBI or DBN crystal structure. This choice assumes that these surface patches are most chemically relevant for crystal formation, which is reasonable for such a small protein and is consistent with earlier studies of protein crystallization 44. We then obtain five patches for DBI, labeled with Arabic numerals, and five patches for DBN, labeled with Roman numerals. Because Patch 4 of DBI contains the 23rd residue, which is associated with the inverted solubility regime, this patch is taken to be temperature dependent (see Sec. 3.2); other patches have a constant energy. Patch energies and interaction ranges were previously extracted from all-atom molecular dynamics simulations 45, using umbrella sampling 46. The resulting effective single-component system model thus coarse-grains the role of solvent and ions in the crystallization cocktail.

(The resulting patchy particle model is sketched in Fig. 3, and the geometry details are given in the supplementary information.) In what follows, unless otherwise specified, energies are reported in units of $k_{\mathrm{B}}T_{\mathrm{ref}}$, where $T_{\mathrm{ref}}=277$K is the temperature at which DBN was crystallized experimentally, and distances are reported in units of the particle diameter $\sigma$, which here is taken to be 2.54 nm.

It is important to highlight that this protocol presents a number of limitations, including inaccuracies of the protein force field 47 and of the water model 48, inefficient sampling, as well as the crudeness of representing potentials of mean force as square well interactions and proteins as spheres. In addition, determining the potential of mean force for each crystal contact is a computationally challenging task, and the 10 ns sampling of each umbrella window likely incompletely explores some of the protein conformational changes, such as loop motion 49. On the whole, this approach likely yields estimates of protein-protein interactions that are at best within 10 to 50% of the association free energy. If sufficiently robust, the resulting phase diagram should therefore be qualitatively, although not quantitatively captured.

3.2 Inverted Solubility Models

In order to represent the microscopic origin of the inverted solubility, we consider three models for the temperature-dependence of Patch 4: the MLG model, the Wentzel-Gunton model, and the temperature-(de)activated patchy model. Note that the parameters explicitly defined in these models are provided here, while the free parameters are left for discussion in Sec. 4

MLG model– In this model, each of four water states is assigned an energy, $E$, and a (relative) degeneracy, $q$. Degeneracies are ordered $q_{\mathrm{ds}}>q_{\mathrm{db}}>q_{\mathrm{ob}}>q_{\mathrm{os}}$. The last inequality follows from the hydration shell allowing the formation of only hydrogen bonds between water molecules and not to the hydrophobic solute. The higher degeneracy of the disordered shell compared to the disordered bulk follows from the additional orientational constraints in the former compared to the latter. Because only relative information about the degeneracies is needed, the estimates of Ref. 50 here suffice: $q_{\mathrm{ob}}=1.5$, $q_{\mathrm{db}}=30$, $q_{\mathrm{os}}=1$, and $q_{\mathrm{ds}}=48$. The energies are ordered $E_{\mathrm{ds}}>E_{\mathrm{db}}>E_{\mathrm{ob}}>E_{\mathrm{os}}$. The ordered shell is expected to have a lower energy than the ordered bulk state, because hydrogen bonds that form via tangentially oriented water molecules tend to be stronger than radially oriented ones. The disordered shell is expected to have a higher energy than the disordered bulk because replacing the solute with water molecules increases slightly the number of hydrogen bonds. Because energy values used by Ref. 25 are on an arbitrary scale, which is incompatible with the specific energy scale of our patchy model, we use instead $E$ and $q$ reported by Silverstein et al. for the Mercedes-Benz model of water 51, 50. Posing that the energy of the ordered bulk is about one hydrogen bond, $E_{\mathrm{ob}}=-5.82$ $k_{\mathrm{B}}T_{\mathrm{ref}}$ 52, the other three states have: $E_{\mathrm{db}}=-1.69$ $k_{\mathrm{B}}T_{\mathrm{ref}}$, $E_{\mathrm{os}}=-5.90$ $k_{\mathrm{B}}T_{\mathrm{ref}}$, and $E_{\mathrm{ds}}=-0.56$ $k_{\mathrm{B}}T_{\mathrm{ref}}$. The energy and entropy per water molecule in the shell and are then given as 23

[TABLE]

and

[TABLE]

The differences in energy and entropy per water molecule in the solvation shell of the protein and in the bulk are then simply $\varepsilon_{w}=E_{s}-E_{b}$ and $\Delta s_{w}=S_{s}-S_{b}$, respectively.

With this formulation the energy of Patch 4 is given by

[TABLE]

where we have defined $\Delta\varepsilon(\beta)=2(\varepsilon_{w}-\Delta s_{w}/\beta)$, and $n_{w}$ is the number of water molecules in the solvation shell around contact $i$. Note that because patch parameters are measured at $\beta_{\mathrm{ref}}=1$, parameters need to be tuned such that $\varepsilon_{4}^{\prime}(\beta=1)=\varepsilon_{4}$, and hence $\varepsilon_{4}^{\prime}=(\varepsilon_{4}-\Delta\varepsilon(1)n_{w})+\Delta\varepsilon(\beta)n_{w}$. Note also that the temperature scale for the MLG model cannot be changed arbitrarily by changing $\beta_{\mathrm{ref}}$, because its parameters already set the range of temperatures within which the hydrophobic effect changes the free energy of crystallization.

Wentzel-Gunton Model– Wentzel and Gunton proposed a simplified version of the MLG model in order to consider the phase behavior of anisotropic particles using Wertheim’s theory 53, 54, 55, 24. This simple model assigns a linear temperature dependence for the patch energies

[TABLE]

where $-\varepsilon_{w}$ and $-\Delta s_{w}$ are free parameters that account for the change in energy and in entropy, respectively, due to the displacement of water upon contact association. Patch energies should equal those of the original model at $\beta=\beta_{\mathrm{ref}}$, at which the model is parameterized. By contrast to the MLG model, this choice here suffices to set the overall temperature scale, because $\varepsilon_{w}$ and $\Delta s_{w}$ are arbitrary. Fixing $\varepsilon_{w}$, such that $\varepsilon_{4}^{\prime}(\beta=\beta_{\mathrm{ref}})=\varepsilon_{4}$, thus results in $\varepsilon_{4}^{\prime}=\varepsilon_{4}+2\Delta s_{w}(\frac{1}{\beta_{\mathrm{ref}}}-\frac{1}{\beta})$.

Temperature-(de)activated Patchy Model– de Las Heras and de Gama 56 proposed a model for patch (de)activation with temperature inspired by DNA-grafted colloids, which lose their attractive patches above the DNA melting temperature 57, Although this model does not correspond to a specific microscopic scenario in proteins, it can nevertheless be construed as a simple and elegant way to describe patch (de)activation. The temperature dependence of the interaction is then

[TABLE]

where $T_{a}$ is the deactivation temperature, $\tau$ controls the sharpness of that deactivation. For this model, Patch 4 is deactivated below $T_{a}$.

3.3 Crystal Solubility Determination

Solubility lines are determined by first calculating the fluid and crystal chemical potentials, and then identifying the coexistence points as the intersection of these curves at fixed temperature and pressure. For both DBI and DBN, experimental solubilities correspond to protein volume fractions of $\phi=10^{-3}$ or lower 29. At such low-densities simple local Monte Carlo (or molecular dynamics) sampling of the fluid phase is computationally inefficient, because transport is relatively slow. While this problem can be alleviated with advanced sampling methods such as aggregation volume bias Monte Carlo 58 and event chain Monte Carlo 59, we here instead estimate the fluid properties from the second virial coefficient, $B_{2}$, which is calculated as in Ref. 60 (see SI 61). Because the patch energies are high, $B_{2}$ can become very large and negative at low temperatures, but the protein density remains sufficiently low for $|B_{2}\rho|\ll 1$ in the regime of interest. In order to confirm that the third virial coefficient, $B_{3}$, can then be neglected, we bound its value as follows. Because the patch properties are such that triply-bonded triplets cannot form, the dominant contribution to $B_{3}$ comes from doubly-bonded triplets. This term scales as $B_{2}^{2}$, hence $|B_{3}|\lesssim|B_{2}|^{2}\ll 1$ in the regime of interest, and its contribution to $\beta\mu_{f}$ is negligible. Higher-order corrections are similarly expected to be negligible, thus justifying this theoretical expediency.

The fluid equation of state and chemical potential, $\mu_{f}$, can then be written as

[TABLE]

[TABLE]

where $\beta\mu^{\mathrm{id}}=\log\Lambda^{3}\rho$ is the chemical potential of the ideal gas, and the thermal de Broglie wavelength $\Lambda$ is set to unity, without loss of generality. With this formulation, we have

[TABLE]

If $B_{2}$ is positive, $\partial\beta\mu_{f}/\partial\rho$ is also positive. If $B_{2}$ is negative, $\partial\beta\mu_{f}/\partial\rho>0$ for $\rho<-1/2B_{2}$, which is always true. Thus, $\beta\mu_{f}$ decreases with decreasing pressure.

The crystal free energy at a given pressure and temperature is calculated using numerical simulations (see SI for simulation details). and the Frenkel-Ladd method 62, which involves thermodynamically integrating from an ideal Einstein crystal. From this reference free energy, thermodynamic integration along an isobar provides the free energy at different temperatures,

[TABLE]

where $\langle H\rangle=p\langle V\rangle+\langle U\rangle$ is the enthalpy and $\langle\cdot\rangle$ denotes thermal averaging. Because of the highly constrained geometry of the patchy models, both crystals are almost incompressible. As a result, $\langle V\rangle$ is essentially independent of temperature. To high accuracy, we can thus write

[TABLE]

Note that at sufficiently low pressures, the second term is also negligible.

We further approximate that all the crystal bonds are active, and hence $\langle U(\beta)\rangle\approx U_{0}(\beta)$, where $U_{0}(\beta)$ is the ground state energy, and $\langle dU/d\beta\rangle\approx dU/d\beta$. While this approximation is generally quite good, it is overly crude in the patch deactivation regime, where the patch energy decreases rapidly around $\beta_{a}$, and vanishes when temperature is reduced further. As a result, $\langle dU/d\beta\rangle\ll dU/d\beta$, which can result in a significant correction to $\beta\mu_{c}$ (see Fig. 4a). In the Wentzel-Gunton model, the patch similarly becomes non-attractive for $\beta>\beta_{\mathrm{ref}}$, and upon further lowering the temperature, it eventually becomes repulsive. The topology of the DBI crystal then changes and the energy of the crystal once again becomes temperature-independent, which leads to a bending of the evolution of the chemical potential with temperature (Fig. 4b). In both cases, however, the DBI solubility curve is unaffected, because these changes occur in a region where DBI is metastable with respect to DBN.

Under this approximation, the chemical potential of the crystal for the MLG model can be written as

[TABLE]

where

[TABLE]

We thus have that

[TABLE]

which has a minimum when

[TABLE]

As already noted, $\beta\mu_{f}$ decreases with decreasing pressure, and because by thermodynamic stability so does $\rho$, an inverted solubility regime is only obtained when the slope of $\beta\mu_{c}$ with respect to $\beta$ is positive. For $\Gamma(\beta)>\varepsilon_{\mathrm{tot}}/(2n_{w})-\Delta\varepsilon(1)/2$, the slope of $\beta\mu_{c}$ is positive, hence inverted solubility is observed.

For the Wentzel-Gunton model, the change in $\beta\mu_{c}$ with temperature can be similarly estimated. We can write the energy per particle in the crystal as

[TABLE]

and hence, following Eq. (23),

[TABLE]

The slope of $\beta\mu_{c}$ with respect to $\beta$ is positive when $-\varepsilon_{\mathrm{tot}}-p/\rho>2\Delta s_{w}/\beta_{\mathrm{ref}}$, thus resulting in inverted solubility.

Although it is not possible to write a compact expression for $\beta\mu_{c}$ for the temperature-(de)activated patchy model–the associated integrals need to be evaluated numerically– the phenomenology is similar. The solubility is inverted in the region where $\beta\mu_{c}$ has a positive slope, i.e., around $T_{a}$. The associated jump in $\beta\mu_{c}$ due to patch deactivation with increasing $\beta$ can be seen in Fig. 4a.

If patch energies are modified by either randomly perturbing them or by scaling them by a constant factor, the free energy of this altered model is estimated from the original model, assuming that the crystal free energy can be expressed as

[TABLE]

where $A^{\prime}$ is the Helmholtz free energy and $U_{0}^{\prime}$ is the ground state crystal energy for altered model. This treatment amounts to neglecting the change in crystal entropy upon weakening or strengthening the patches, which is a small contribution. We also verify that the crystal remains stable at the temperatures of interest.

The approximations described above allow the expedited consideration of coexistence points that constitute the solubility curves by generating $\beta\mu_{f}$ and $\beta\mu_{c}$ as functions of $\beta$ at various pressures.

4 Inverted Solubility from Hydrophobicity Models

In order for the microscopic hydrophobicity models described above to give rise to solubility inversion, a sufficient number of water molecules need to be involved. In this section we first consider physical bounds on that number, and then consider how the corresponding crystal solubility lines are affected.

4.1 Effect of Parameters on Solubility Lines

The key free parameter in hydrophobicity models is the number of water molecules solvating the hydrophobic patch, $n_{w}$. We first estimate the number of water molecules possibly available around Patch 4, by calculating the SASA for the participating residues 30 and then compute

[TABLE]

where $A_{4}$ is the solvent accessible surface area of Patch 4, $\rho_{w}=3.3\times 10^{-2}\text{\AA}^{-3}$ is the number density of water in the bulk at room temperature, and $g_{C}(r)$ is the radial distribution function of water around carbon atoms determined in Ref. 48. This estimate thus assumes that (i) the solvent has a radius of $1.4\text{\AA}$ (the SASA definition), (ii) the average van der Waals radii of protein heavy atoms is $\sim 1.6\text{\AA}$, and (iii) the first solvation shell ends with the first peak of $g(r)$ at $4.5\text{\AA}$. We further assume that the measured surface is flat, which is here but a small correction. If we furthermore assume that all residues contributing to Patch 4 are hydrophobic, then $n_{w}=133-140$ for all six protein structures. However, because Patch 4 contains only a handful of hydrophobic residues a more realistic estimate should decrease $A_{4}$. Assuming that a residue is hydrophobic if it is labeled as such in any of the hydrophobicity scales considered in Sec. 2 gives instead $n_{w}=43-48$. Because the hydrophobic residues within Patch 4 do not constitute a contiguous area, the configuration of contributing water molecules solvating them will additionally be affected by the nearby hydrophilic surface residues. This estimate should thus be treated as an upper bound. Note that the P23T mutation does not seem to be associated with a systematic change in $A_{4}$, and thus $n_{w}$.

We compare this result with the number of water molecules needed for Patch 4 to have its measured bond strength. In particular, if we attribute the entire Patch 4 energy to the change in free energy upon moving solvating water molecules to the bulk, then the MLG model gives $\varepsilon_{4}=n_{w}\Delta\varepsilon(1)$, and thus $n_{w}\sim 23$. Because multiple hydrogen bonds also contribute, however, this number should also be treated as an upper bound, that is consistent with yet tighter than the above bound.

We can now contrast these bounds with the minimum number of water molecules, $n_{w}^{*}$, that need to be displaced to invert solubility. For the MLG model, we use Eq. (29) and the sum of DBI patch energies, $\varepsilon_{\mathrm{tot}}=60$, to estimate $n_{w}^{*}$; it must be such that $\xi(\beta_{\mathrm{min}})$ is a minimum, i.e., $\Gamma(\beta)>\varepsilon_{\mathrm{tot}}/(2n_{w})-\Delta\varepsilon(1)/n_{w}$. In other words, the solubility is inverted if $\beta>\beta_{\mathrm{min}}$. The numerical solution in Fig. 5 shows that $n_{w}^{*}\gtrsim 71$. The corresponding change in $\xi(\beta)$ is given in the inset. It should be noted, however, that $n_{w}^{*}$ depends on the MLG model parameters. For the multiplicities of Shiryayev et al., for instance, inverted solubility is possible with a mere $n_{w}^{*}\sim 18$. These degeneracies, however, are fairly arbitrary 25. Choosing $q_{\mathrm{ob}}=10$ and $q_{\mathrm{os}}=1$, in particular, seems unphysical, because the ordered bulk degeneracy is unlikely to be an order of magnitude larger than that of the ordered shell. We thus expect $n_{w}^{*}\gtrsim 71$ to be physically more reasonable.

However, because $\varepsilon_{\mathrm{tot}}=60$ results in room temperature solubilities that are orders of magnitude lower than their experimental counterpart, and in light of the various sources error in patch energy determination (Sec. 3.1), Ref. 30 proposed a reasonable correction would be to halve patch energies. For $\varepsilon_{\mathrm{tot}}=30$, we have $n_{w}^{*}\sim 35$ (Fig. 6), which is less than the 40 or so water molecules solvating hydrophobic residues in Patch 4, but more than the energetic estimate.

In light of the many estimates involved in the above analysis, the hydrophobic effect as a cause of inverted solubility, although weakly supported, cannot be eliminated outright. Even if the hydrophobicity model parameters are kept constant, a possible resolution could be for Patch 4 to be stronger than estimated and the other patches weaker. The hydrophobicity scenario, however, does severely constrain the patch model parameters.

4.2 Solubility Lines for Models of Hydrophobicity

In the previous subsection we determined that the hydrophobicity scenario for inverting solubility requires a fine balance between the protein-protein patch energies, the size of the hydrophobic patch, and the number of water molecules solvating it. While this rare confluence of factors could explain why inverted solubility is not common among proteins, it is natural to wonder whether the presence of weak hydrophobic patches, which are ubiquitous in proteins, affects solubility lines without engendering a regime of inverted solubility. In this section, we study the Wentzel-Gunton model in order to examine this possibility.

The solubility lines for the Wentzel-Gunton model in Fig. 7 are specifically obtained for $\beta_{\mathrm{ref}}=0.5$, to match the experimental solubility as in Ref. 30, but our observations are qualitatively independent of this choice. Setting $\Delta s_{w}=-10$, which is here akin to $n_{w}\sim 20$ (assuming that the temperature-dependent energy in the MLG model scales as $n_{w}$), results in normal solubility, but the steepness of the solubility curve changes markedly compared to $\Delta s_{w}=0$. Setting $\Delta s_{w}=-15$ ($n_{w}\sim 30$) results in the DBI solubility being almost independent of temperature and in DBN being more stable than DBI at $T<T_{\mathrm{tp}}\sim 1.7$. Further reducing $\Delta s_{w}$ results in an inverted solubility regime. The solubility curve then flattens below $T\sim 2$ and $\phi_{\mathrm{tp}}$ moves to higher packing fractions. These observations thus reveal that the presence of an inverted solubility regime is the limit case of a continuum of how hydrophobicity impacts the solubility line.

5 Solubility Lines for Temperature Deactivated Patches

Absent clear microscopic evidence for the hydrophobic effect, we finally consider a generic model for patch deactivation. The temperature-deactivated patchy model, which was used to successfully capture the inverted solubility of DBI 30, stabilizes the crystal with increasing temperature without referring to any specific microscopic mechanism. In this section, we first discuss the physical constraints on the model parameters and then consider how solubility lines change with model parameters, paying particular attention to the robustness of the inverted solubility regime. We also estimate the binodal and the critical temperature, which has been experimentally for some human $\gamma$D-crystallin variants 31.

5.1 Parameter Estimates

Despite the absence of an explicit microscopic interpretation for the (de)activation model, one can still place some reasonably solid physical constraints its tuning parameters. First, the (de)activation temperature $T_{a}$ must lie in the vicinity of the triple point, and thus $T_{a}\sim T_{\mathrm{tp}}$. For our model, the choice $T_{a}=1.9$ ensures that the deactivation of Patch 4 makes DBI metastable with respect to DBN for $T<T_{\mathrm{tp}}$. Second, $\tau$, which sets the temperature range over which (de)activation takes place, ought to capture the degree of cooperativity of the underlying microscopic effect. It therefore cannot be arbitrarily small, as it would be at a thermodynamic phase transition. Denaturing a protein, for instance, takes place over a few degrees, and any smaller scale rearrangement that involves tens to hundreds of atoms should spread over at least $\gtrsim 10K$. However, if this change takes place too gradually, say over $\gg 10K$, then inverted solubility cannot be observed for typical protein-protein interactions. We thus here consider a temperature range of $\sim 10K$, which corresponds to setting $\tau=0.05$.

We first investigate how varying the patch energies impacts the phase diagram, keeping $T_{a}=1.9$ and $\tau=0.05$ constant. As previously reported 30, the resulting phase diagram (Fig. 8a) exhibits a re-entrance regime bounded by the DBI solubility line, as well as a triple point between the fluid and the two crystal forms. The solubility lines that result from perturbing the patch energies by $5\%$ and $10\%$ are shown in Fig. 8a. Although the solubility lines then shift to substantially lower or higher densities, the existence of an inverted solubility regime is robust. The errors inherent to the overall parameterization of the model are therefore qualitatively benign.

We next investigate the robustness of the results with respect to the relative strength of the temperature-deactivated fourth patch, $\varepsilon_{4}$. This question is of interest for two main reasons: (i) the strength and robustness of solubility inversion depend sensitively on the strength of that patch; and (ii) the ordering of the single mutant solubilities directly correlates with their respective Patch 4 energies.

We first multiply $\varepsilon_{4}$ by $f\in\{0.9,1.0,1.5\}$, while keeping the other patch parameters constant (Fig. 8b). Increasing the strength of Patch 4 systematically decreases the solubility of DBI as well as $\phi_{\mathrm{tp}}$. Interestingly, the decrease in solubility with increasing $f$ is consistent with the experimental observations for the single mutants, P23T, P23S, and P23V 31. Because a stronger Patch 4 decreases the DBI solubility (R36S+P23T double mutant), assuming that the difference between crystals will arise due to Patch 4 only, we speculate that if two other double mutants, R36S+P23S and R36S+R23V, were crystallized with similar crystal contacts, then their inverted solubility would have a similar ordering.

We next change the energy of the temperature-dependent patch while keeping the total energy of DBI patches constant, i.e., $\varepsilon_{\mathrm{tot}}=f_{1}\varepsilon_{4}+f_{2}(\varepsilon_{1}+\varepsilon_{2}+\varepsilon_{3}+\varepsilon_{5})$ (Fig. 8c). (Because the second patch corresponds to a shared contact between DBI and DBN, the DBN solubility is then also slightly perturbed.) As in the first case, the inverted solubility regime vanishing upon markedly reducing the strength of Patch 4. The difference is that DBN is here metastable with respect to DBI within the probed temperature range, while DBI becomes metastable with respect to DBN otherwise. For $f_{1}=0.4$, DBN is still metastable with respect to DBI, but the narrow range of inverted solubility is then replaced by standard solubility at both lower and higher temperatures.

We finally investigate the robustness of the phenomenology with respect to changes in $\tau$. Decreasing $\tau$ corresponds to a faster temperature (de)activation of the patch, which flattens the inverted solubility region and results in $T_{\mathrm{tp}}\rightarrow T_{a}$ as $\tau\rightarrow 0$ (Fig. 9a). $\phi_{\mathrm{tp}}$ similarly gets pushed to higher packing fractions, suggesting that a protein solution prepared very near $T_{a}$ could reach remarkably high concentrations compared to solutions prepared at surrounding temperatures. However, as argued above, very small values of $\tau$ are not here physically meaningful. Conversely, increasing $\tau$ weakens this transition and eventually eliminates the solubility inversion regime. Increasing $\tau$ also renders the inverted solubility less robust with respect to varying patch energies. Interestingly, a specific choice of $\tau$, with a minor tweak to patch energies ($\tau=0.28$, $\varepsilon_{i}^{\prime}=1.1\varepsilon_{i}$), gives rise to a nearly vertical solubility curve (Fig. 9b), similar to the temperature-independent solubility of apoferritin 21.

5.2 Estimation of the Critical Temperature

Although various theoretical results suggest that a closed-loop binodal with multiple critical points is possible upon introducing temperature-dependent binding energies 23, 24, 56, no experimental evidence of such a closed-loop binodal is found for any human $\gamma$D-crystallin mutant. In addition, experiments find that the P23V mutation, which also inverts solubility, has a binodal indistinguishable of that of the wild type 31. We thus estimate the liquid-liquid binodal and the associated critical temperature, $T_{c}$, to determine if the deactivation of Patch 4 affects the solution properties. Here we use Wertheim’s perturbation theory 53, 54, which provides quantitatively good estimates of the binodals in patchy models (see SI for details 63).

Choosing $\tau=0.05$ and $T_{a}=1.9$, as above, results in a typical binodal with a single critical point at $T_{c}=1.85$ (Fig. 10). Hence, without altering patch energies, we do not observe any signature of a closed-loop binodal in our model (Fig. 10), which is consistent with experimental results 31. In order to determine how far our model is from exhibiting a closed-loop binodal, we systematically increase the Patch 4 energy. Only above $\varepsilon_{4}>36$, does a closed-loop binodal appear. See example in Fig. 10. If only the energy of Patch 4 is more than doubled and is disproportionately stronger than any other patch in the model does this phenomenon appear. This perturbation falls far outside of the error estimates for patch energies, which supports the robustness of our model prediction.

6 Conclusion

We have here attempted to rationalize the inverted solubility of certain mutants of $\gamma$D-crystallin based on microscopic models of protein-protein interactions and their temperature dependence. We have paid particular attention to the putative role of hydrophobicity. Estimating surface hydrophobicity using different scales did not reveal the presence of any pertinent surface feature, but microscopic models of hydrophobicity suggest that the amount of available surrounding water molecules might be close be sufficient. Although our analysis falls short of concluding whether hydrophobicity plays a determining role, the scenario nonetheless seems a bit far fetched. A more convincing determination would likely require for the water structure around the region of interest to be more specifically probed. Because numerical simulations of standard model of water are insufficiently sensitive to this feature 48 to be conclusive, simulations with more sophisticated water models 64, 65, and neutron diffraction or hydrogen-deuterium exchange experiments might be more productive avenues.

Even though the microscopic origin of inverted solubility in human $\gamma$D-crystallin thus still remains somewhat elusive, additional insight from the crystallization of other double mutants, such as R36S+P23V and R36S+P23S, might be helpful in identifying generic features that might have eluded the analysis thus far. Repeating the above structural and thermodynamical for these mutants could help tease out more subtle features.

{acknowledgement}

We wish to thank our colleague David Beratan for his continued support. It has been a pleasure to exchange fruits and science with him over the last decade. We acknowledge support from National Science Foundation Grant no. NSF DMR-1749374. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1548562, as well as the Duke Compute Cluster. We warmly thank Diana Fusco for stimulating discussions as well as Amir Khan for discussing and sharing experimental data ahead of publication. Data relevant to this work have been archived and can be accessed at https://doi.org/10.7924/XXXXXXX.

{suppinfo}

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/XXXX.

Supporting Information

Analyzing Crystal Structures

The residues identified to be involved in Patch 4 are:

LYS 2, TYR 16, GLU 17, SER 19, SER 20, ASP 21, HIS 22, PRO/THR 23, ASN 24, PRO 27, TYR 28, PRO 48, ASN 49, TYR 50, PRO 82, HIS 83, SER 84, SER 86, ARG 88, TYR 97, VAL 125, LEU 126, GLU 127, GLY 128, SER 129, MET 146, PRO 147, GLY 148, ASP 149, TYR 150, ARG 151, ARG 152, GLN 154, ASP 155.

Residues in bold are deemed hydrophobic in at least one hydrophobicity scale. We compare the solvent accessible areas (SASA) obtained from these selections, $A_{4}$ , across the available mutant crystal structures to evaluate if a trend in hydrophobicity that could explain the difference between normal and inverted solubility could be measured. If the P23T mutation were to increase the SASA, then maybe a tipping $n_{w}$ could be reached. The bare value of $A_{4}$ gives R36S $>$ DBI∗ $>$ R58H $>$ WT $>$ DBN∗ $>$ P23T∗, where the asterisk denotes structures that contain the solubility inverting mutation. If we consider only the hydrophobic residue contribution to Patch 4, we instead have R58H $>$ R36S $>$ DBI∗ $>$ DBN∗ $>$ WT $>$ P23T∗. In any case, no systematic change is detected, thus further weakening the case for hydrophobicity as key to solubility inversion.

Model Parameters and Geometry

Model parameters are obtained from all-atom molecular dynamics (MD) simulations using Gromacs 45 versions 5.1 and 2016. Patch energies are chosen as the depth of the well in the PMF, after zeroing its long distance energy value. The interaction range for each patch was then calculated by matching the second virial coefficient, $B_{2,\alpha\beta}$, of the square well potential to that of the PMF

[TABLE]

from which the interaction range is found to be

[TABLE]

The details of the all-atom simulations resulting in model parameters are given in the supporting information of Ref. 30. The resulting model parameters are summarized in Table 1. Contacts are labeled cY, where Y is a Roman numeral for DBN contacts, and an Arabic numeral for DBI contacts. Contact 2 (c2) of DBI and contact I (cI) of DBN are very similar in terms of geometry and the nature of their interactions and were hence merged. Contact V (cV) and VII (cVII) of DBN are also almost identical and were also merged.

The patch locations and dihedral angles are listed in Table 2. A number of modifications were made to the geometry of the actual protein crystal to be able to represent the proteins as spherical objects in an orthorhombic unit cell and are detailed in the SI of Ref. 30. Here we briefly quote the resulting geometry (Table 2), and note that we do not expect these topological changes to qualitatively affect our results, as it was previously shown that patch locations have only a negligible effect on the phase diagram topology 61. Note that the resulting patch geometry is such that dihedral constraints are necessary to properly separate the crystal ground states within the energy landscape.

Determining Solubility Lines

For the patchy model studied here, the second virial coefficient, $B_{2}$, is given as

[TABLE]

where the integrand is the Mayer function, and the integration is performed over particle positions, $\mathbf{r}$, orientations $\mathbf{\Omega}$, and dihedral $\chi$. Computing this integral for our patchy model gives

[TABLE]

where the sum is over all patch pairs excluding cII, cIV, and cVI, which are collocated. The terms after the sum specifically describe the contribution of these patches. Note that $2\Delta\chi=0.6$ is the width of range of dihedral angles for these three patches, and $\Delta\phi_{V}=1.0632$ because there are four valid dihedral angles for cV with overlapping ranges.

Specialized MC simulations and estimates using the second virial coefficient were used to determine solubility lines. MC simulations have 200,000 to 1,000,000 MC sweeps, each of which consists of $N$ displacement and $N$ rotational moves for constant $NVT$ simulations, and two additional volume moves for $NPT$ simulations. The amplitude of displacements is chosen such that an acceptance rate of about 40% is obtained. System sizes are chosen to be similar from one phase to the other, while respecting the crystal symmetry, here, $N_{\mathrm{DBI}}$=432 and $N_{\mathrm{DBN}}$=512. For convenience, these simulations report distances in units of the particle diameter, $\sigma=2.54$nm, determined as the shortest center of mass distance between two proteins in contact, and energies in units of $k_{\mathrm{B}}T_{\mathrm{ref}}$ with $T_{\mathrm{ref}}=277$K.

Estimating the Critical Temperature with Wertheim’s Perturbation Theory

Because Monte Carlo sampling gets increasingly inefficient as temperature is decreases, the critical temperature and the liquid-liquid coexistence regimes are determined using Wertheim’s perturbation theory. The fluid free energy is then approximated as

[TABLE]

where $A_{\mathrm{HS}}$ is the free energy of the hard sphere fluid and $A_{\mathrm{bond}}$ is the bonding contribution

[TABLE]

where the sum runs over all patches, $X_{a}$ is the bonding probability of patch $a$, and $M$ is the number of patches. Bonding probability is calculated as

[TABLE]

where $K_{a}$ quantifies the interaction strength of patch $a$ and is calculated as in Ref. 63,

[TABLE]

where $g_{\mathrm{ref}}^{c}$ is the contact value of the radial distribution function of a reference fluid, which we take to be the hard sphere fluid for simplicity, and $\varepsilon_{a}$, $\delta_{a}$, $\delta_{a}^{\prime}$, $\lambda_{a}$, and $\Delta\chi_{a}$ are patch parameters. The $(\Delta\chi_{a}/\pi)$ term takes into account the reduced bonding volume due to the dihedral constraint. The critical temperature, $T_{\mathrm{c}}$ is then calculated by setting

[TABLE]

The binodal is then traced using a Maxwell construction at various temperatures.

Note that for this calculation patches cII, cIV, and cVI, which are collocated, are treated by using $\lambda_{II}$, $\varepsilon=\varepsilon_{II}+\varepsilon_{IV}+\varepsilon_{VI}$, and $\delta_{II}=\delta_{II}^{\prime}$ as the patch parameters.