Packing in Protein Cores

TL;DR

This paper investigates the packing density and structural organization of amino acids in protein cores, revealing a lower packing fraction than traditionally thought and comparing it to jammed particle packings, with implications for protein stability and design.

Contribution

It demonstrates the importance of calibrated atom sizes and explicit hydrogen atoms in modeling protein core packing, providing new insights into amino acid arrangements and packing density.

Findings

Protein cores have a packing fraction of approximately 0.56.

Amino acids pack as densely as disordered jammed particles with similar shapes.

Proper modeling requires explicit hydrogen atoms and calibrated atom sizes.

Abstract

Proteins are biological polymers that underlie all cellular functions. The first high-resolution protein structures were determined by x-ray crystallography in the 1960s. Since then, there has been continued interest in understanding and predicting protein structure and stability. It is well-established that a large contribution to protein stability originates from the sequestration from solvent of hydrophobic residues in the protein core. How are such hydrophobic residues arranged in the core? And how can one best model the packing of these residues? Here we show that to properly model the packing of residues in protein cores it is essential that amino acids are represented by appropriately calibrated atom sizes, and that hydrogen atoms are explicitly included. We show that protein cores possess a packing fraction of $\phi \approx 0.56$, which is significantly less than the typically quoted value of 0.74 obtained using the extended atom representation. We also compare the results for the packing of amino acids in protein cores to results obtained for jammed packings from disrete element simulations composed of spheres, elongated particles, and particles with bumpy surfaces. We show that amino acids in protein cores pack as densely as disordered jammed packings of particles with similar values for the aspect ratio and bumpiness as found for amino acids. Knowing the structural properties of protein cores is of both fundamental and practical importance. Practically, it enables the assessment of changes in the structure and stability of proteins arising from amino acid mutations (such as those identified as a result of the massive human genome sequencing efforts) and the design of new folded, stable proteins and protein-protein interactions with tunable specificity and affinity.

Tables

Table 1. Table 1: Atomic radii used in the hard-sphere model and three other studies (one using explicit hydrogens  [ 35 ] and two others using the extended atom model  [ 22 , 23 ] ). All values are given in Å.

Atom Type	Hard-sphere Model	Word 1999 [35]	Richards 1974 [22]	Liang 2001 [23]
C${}_{s{p}^3}$	1.5	1.65	2.0	1.88
Caromatic	1.5	1.65	1.7	1.61 - 1.76
CO	1.3	1.65	1.7	1.61 - 1.76
N	1.3	1.55	1.7	1.64
O	1.4	1.4	1.4 - 1.6	1.42 - 1.46
S	1.75	1.8	1.8	1.77
H	1.1	1.17	N/A	N/A

Table 2. Table 2: Volumes for the 11 11 11 residues that occur most frequently in protein cores using the explicit hydrogen representation.

Residue	Volume (Å3)
Ala	48.8
Cys	64.3
Gly	35.6
Ile	88.1
Leu	88.1
Met	92.7
Phe	100.7
Thr	69.0
Trp	121.9
Tyr	107.5
Val	75.0

Equations

$${\varphi }_r=\frac{{\sum }_i{\mathrm{V}}_i}{{\sum }_i{\mathrm{V}}_i^v},$$

$$U_{\mathrm{RLJ}}(r_{ij})=\frac{ϵ}{72}{\left[1-{\left(\frac{{\sigma }_{ij}}{r_{ij}}\right)}^6\right]}^2\Theta ({\sigma }_{ij}-r_{ij}),$$

$$\Delta \chi =\sqrt{({\chi }_1^{xtal}-{\chi }_1^{HS}{)}^2+\dots +({\chi }_n^{xtal}-{\chi }_n^{HS}{)}^2}.$$

$$B=\sqrt{\text{bigint}d\hat{u}\frac{{\left(\vec{R}(\hat{u})-\vec{\mathcal{R}}(\hat{u})\right)}^2}{{\mathcal{R}}^2(\hat{u})}},$$

Full text

Keywords: proteins, random close packing, jamming

1 Introduction

Proteins are biological polymers that play important roles in cellular processes ranging from the purely structural to the actively catalytic. Proteins are linear chains of different combinations of the $20$ naturally occurring amino acid residues with variable chain lengths from tens to tens of thousands. A key feature that distinguishes proteins from other polymers is that each folds into a unique three-dimensional structure. Proteins typically fold spontaneously in aqueous solution at room temperature. The amino acid sequence is the only information required to specify a protein’s unique structure [1, 2].

The amino acids can be grouped into two main categories: hydrophobic and hydrophilic. Hydrophobic residues form the solvent-inaccessible core of a protein and hydrophilic residues, both polar and charged, are on the solvent-accessible surface. As of 2017, the structures of more than 125,000 proteins have been determined, primarily by x-ray crystallography, with a median resolution of $\approx 2.5$ Å and deposited in the protein data bank (PDB) [3]. This large database of atomic coordinates provides a wealth of structural information that can be used to analyze the physical properties of proteins and to understand how proteins interact and carry out their functions [4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14].

Each amino acid is made up of the same backbone unit of four heavy (non-hydrogen) atoms, N-Cα-C-O, and different combinations of side chain atoms that branch from the Cα atom (Fig. 1). The repeating units are joined by a peptide bond between the carboxyl carbon (C) of a given amino acid and the nitrogen (N) of the next. All bond lengths and bond angles are specified by the same basic stereochemistry that defines the structures of small molecules [15, 16]. The three-dimensional structure that a protein adopts is specified by the amino acid dihedral angles. For each amino acid in the protein chain, there are two backbone dihedral angle degrees of freedom, $\phi$ and $\psi$, and $N_{s}$ side chain dihedral angle degrees of freedom, $\chi_{1},\ldots,\chi_{N_{s}}$. (See Fig. 1.) $N_{s}$ ranges from zero (for alanine and glycine) to five (for arginine). The third backbone dihedral angle is typically constrained to be $\omega=180^{\circ}$ or $0^{\circ}$. Repetition of certain backbone $\phi$ and $\psi$ values in a stretch of amino acids gives rise to specific secondary structures, such as $\alpha$-helices and $\beta$-sheets [17, 18]. All proteins are formed from different combinations of $\alpha$-helix, $\beta$-sheet, and ‘random coil’ structures. Interactions between different elements of secondary structure are stabilized by interactions between the side chains [19, 20, 21]. In addition, side chain interactions on the surfaces of proteins also specify how different proteins bind to each other and to other molecules [6].

A minimal physical model for an amino acid is a composite particle formed from connected spheres with stereochemical constraints (Fig. 2). As is clear from Fig. 2, amino acids are non-spherical objects with complex shapes. Thus, we can imagine proteins as strings of interconnected non-spherical objects that fold into compact three-dimensional structures. Many prior studies have argued that the cores of folded proteins are tightly packed. For example, several studies have measured the ratio between the volume of a core amino acid and its Voronoi volume to be greater than 70%, which suggests dense crystalline packing [22, 23]. In addition, experimental studies find that mutations in protein cores from small to large residues typically destabilize the protein, suggesting that there is very little empty space present to accommodate additional atoms [24, 25].

In this review, we summarize prior work on the structural properties of protein cores and provide strong evidence that although protein cores are densely packed, they are not as densely packed as crystalline solids. Instead, protein cores possess packing fractions of $\sim 0.56$ [14]. Even though this value is lower than that for crystalline solids (e.g. $0.74$ for face-centered-cubic crystals), protein cores are solid-like with very little free volume that would allow side chain motion. We also show that static packings of particles with complex, non-spherical shapes possess packing fractions below $0.6$, yet still display solid-like properties and that the amino acids in protein cores can be modeled as random, densely packed non-spherical objects. We then relate our computational studies of dense packing in protein cores to experimental studies of mutations that are able to alter the structure and stability of proteins.

2 Packing efficiency in protein cores

By determining the packing fraction of protein cores one can begin to understand their structural and mechanical properties. For example, the shear modulus (i.e. the material response to applied shear stress) typically increases monotonically with the packing fraction since the number of stress-bearing interatomic contacts increases with the packing fraction [26]. Thus, the rigidity of proteins is strongly correlated with the packing density [27, 28]. In addition, knowing the packing density is vital for predicting changes in stability from mutations to protein cores, many of which are disease-associated [29]. Accurate calculations of the packing density are also necessary to predict structure from sequence and to design new stable proteins [10, 30, 31].

One of the first studies of the packing density of protein cores was performed by Richards in 1974. At this time, only a few protein crystal structures were available. Richards focused on two proteins: lysozyme and ribonuclease S [22]. When a protein structure is obtained from x-ray crystallography, the resolution of the structure typically does not allow for the placement of the hydrogen atoms in the protein. In the past, researchers circumvented this problem by implementing an “extended atom” model, where the atomic radii of each heavy atom are increased by a factor that depends on the number of hydrogen atoms that are bound to it [22, 23, 34]. New computational techniques allow for the accurate placement of hydrogen atoms in a protein crystal structure [35, 36], which provides a more detailed “explicit hydrogen” model of proteins. Since hydrogen atoms comprise $\sim 50\%$ of the atoms in a protein, the extended atom approximation can have major effects on the accuracy of the structural model of the protein.

To accurately assess the packing fraction of proteins, one must calibrate and select proper atomic radii. In our recent work [14], we have chosen atomic radii that when used in a hard-sphere model of a dipeptide mimetic can reproduce the observed side chain dihedral angle distributions of non-polar amino acids in a database of high resolution crystal structures [14, 32, 33, 37, 38]. The values for the seven atomic radii are $C_{sp^{3}}$, $C_{\rm aromatic}$: $1.5$ Å; $C_{O}$: $1.3$ Å; $O$: $1.4$ Å; $N$: $1.3$ Å; $H$: $1.10$ Å; and $S$: $1.75$ Å. In Fig. 3, we show that the side chain dihedral angle distributions predicted using the hard-sphere model for a Val and Ile dipeptide agree with the observed side chain dihedral angle distributions. We have shown similar agreement between the observed and predicted side chain dihedral angle distributions for Cys, Leu, Met, Phe, Thr, Trp, Tyr, and Ser [38]. The atomic radii are similar to values of van der Waals radii reported in other studies, and typically smaller than those used in extended atom models (Table 1) [18, 22, 34, 37, 39, 40, 41, 42, 43, 44, 45, 46, 47].

The packing fraction of each residue in a protein core can be calculated using

[TABLE]

where $\mathrm{V}_{i}$ is the ‘non-overlapping’ volume of atom i, $\mathrm{V}^{v}_{i}$ is the Voronoi volume of atom i, and the summations are over all atoms of a particular residue. We also calculate the packing fraction of a protein core, $\phi_{c}$, where both summations are over all atoms of all residues in a particular protein core. Voronoi cells were obtained for each atom using Laguerre tessellation, where the placement of the Voronoi cell wall is based on the relative radii of neighboring atoms (which is the same as the location of the plane that separates overlapping atoms) [14, 48]. $V_{i}$ was calculated by splitting overlapping atoms by the plane of intersection between the two atoms.

Our analysis focuses on residues in protein cores. We have identified all core residues in a database of high resolution crystal structures (described in [14, 32, 33]) using a method described previously [14, 49]. In brief, non-core atoms are identified as those that are on the surface of the protein or near an interior void with a radius $\geq 1.4$ $\mathrm{\AA}$. In this strict definition, a core residue is defined as any residue containing exclusively core atoms (including hydrogen atoms). This method identifies atoms adjacent to voids in the protein and removes them from the calculation of the packing fraction. According to this definition and using the explicit hydrogen representation, the proteins we considered have an average of $15$ core residues of which $80\%$ are Ala, Cys, Gly, Ile, Leu, Met, Phe, and Val.

As shown in Fig. 4 (a), the average packing fraction of protein cores is $\phi_{c}\approx 0.56$ [14]. This value is much closer to packing fractions obtained for jammed packings of frictional or elongated particles rather than $\phi_{c}=0.71$-$0.74$ for packings with significant FCC crystalline order as proposed in earlier studies [22, 23, 34]. (See Section 4.) The most significant difference between the recent and prior studies is the use of a well-calibrated explicit hydrogen model instead of an extended atom model.

To assess the effect of backbone connectivity on the packing efficiency in protein cores, we performed discrete element simulations to compress amino acid monomers into static (i.e. force and torque-balanced) jammed packings. (See the Appendix for a more detailed description of the packing-generation protocol.) We initialized the system by randomly inserting $N_{R}$ residues into a cubic box (with periodic boundary conditions). We assumed that the residues, which are composed of rigidly connected spherical atoms of different sizes, interact via purely repulsive linear spring forces. We then compress the system by small packing fraction increments $\Delta\phi$, followed by energy minimization. For sufficiently small $\Delta\phi$, the form of the purely repulsive potential does not influence the structure of the final packings. For jammed packings, the total potential energy per residue $U/N_{R}>0$ following energy minimization. In contrast, unjammed packings will possess $U/N_{R}=0$ after energy minimization. In this case, atomic motions can occur in the system without a concomitant increase in the total potential energy. Thus, we can identify the packing fraction at jamming onset $\phi_{J}$ as the one at which the minimized total potential increases above a small threshold [50].

We studied mixtures of $N_{R}$ residues with the fractions of Ala, Ile, Leu, Met, Phe, and Val residues matching the percentages found in protein cores. (We focused on non-polar residues, but because Gly has no side chain and Cys can form disulfide bonds, these were not included in the simulations.) These simulations generate disordered jammed packings with $\phi=0.56$ similar to that found in protein cores (Fig. 4 (b)). These results indicate that the connectivity of the protein backbone does not impose significant constraints on the free volume in protein cores.

To further analyze the packing efficiency in protein cores, we also calculated the distribution of the local packing fractions (i.e. $\phi$ for each residue type) in protein cores for both protein crystal structures and simulations. We find that the distributions of the local packing fractions for each residue type have similar average values, differing by $<5$%. In addition, the average values for the local packing fractions are similar to the global average in the core with standard deviations that are slightly larger, which reflects the fact that the local packing fraction is obtained by averaging over fewer atoms than the global packing fraction. We also find that the average packing fraction of each amino acid type is similar to the average packing fraction in protein cores, except for Ala, which does not have a side chain dihedral angle degree of freedom. The similarity of the average packing fraction for individual residues and the average packing fraction in protein cores suggests that there are only small variations of the packing fraction within each protein core.

We also investigated the packing efficiency of protein-protein interfaces. To do this, we compiled a protein-interface database of $123$ crystal structures containing protein-protein and protein-peptide binding pairs. The structures are composed of both homo- and heterodimers with resolution $\leq 1.5$ Å and less than $50$% sequence identity. A core-interface residue is defined as any residue that is a surface residue in the individual protein monomers, but is completely buried after binding. Several studies have shown that the properties of protein-protein interfaces are similar to those of protein cores [8, 51]. Our analyses of protein cores and interfaces confirm this by showing that they possess a similar distribution of amino acids (i.e. primarily hydrophobic residues with few charged and polar residues). We find that 73% and 68% of the residues in protein cores and interfaces, respectively, are hydrophobic with similar frequencies for each amino acid. In addition, both the distribution of core packing fractions and interface packing fractions are peaked near $0.56$ as shown in Fig. 5. This result demonstrates that protein-protein interfaces are packed similarly to protein cores.

3 Protein core repacking

Computational protein core repacking allows investigation of the uniqueness of the side chain conformation of residues in protein cores. Unique side chain conformations for core residues would imply that protein cores are jammed with very little free volume for rearrangements of side chains. There are two categories of protein core repacking investigations: one starts with all possible sequences and seeks to recover the wild type sequence [52, 53] and the other starts with the wild type sequence and seeks to recover the observed combination of side chain dihedral angles and determine if alternative combinations are possible. Here we focus on the second, where the side chains of core residues are removed simultaneously and all side chain dihedral angle combinations of the starting sequence are sampled. The energy of each conformation is evaluated, the optimal conformation is predicted, and then compared to the observed structure.

To study repacking of protein cores, we again use a hard-sphere plus stereochemistry model. The cores of 221 proteins in the Dunbrack Database [32, 33] were studied. As a way to model the system at non-zero temperature and to improve the statistics, variations in bond lengths and angles are implemented by replacing each side chain with different instances of the side chain taken from high-resolution protein crystal structures [4]. Core residues were identified as described in Section 2. As described in previous work [38, 14], the hard-sphere model treats each atom $i$ as a sphere that interacts pairwise with all other non-bonded atoms $j$ via the purely repulsive Lennard-Jones potential:

[TABLE]

where $r_{ij}$ is the center-to-center separation between atoms $i$ and $j$, $\sigma_{ij}=(\sigma_{i}+\sigma_{j})/2$, $\sigma_{i}/2$ is the radius of atom $i$, $\Theta(\sigma_{ij}-r_{ij})$ is the Heaviside step function, and $\epsilon$ is the strength of the repulsive interactions. Values for the atomic radii are listed in Section 2.

Predictions of the side chain conformations of single amino acids are obtained by rotating each of the side chain dihedral angles $\chi_{1},\chi_{2},...,\chi_{n}$ (with a fixed backbone conformation [54]), and finding the lowest energy conformation of the residue, where the total energy $U(\chi_{1},...,\chi_{n})$ includes both intra- and inter-residue steric repulsive interactions. We then calculate the Boltzmann weight of the lowest energy side chain conformation of the residue, $P_{i}(\chi_{1},....,\chi_{n})\propto e^{-U(\chi_{1},...,\chi_{n})/k_{B}T}$, where the small temperature, $T/\epsilon$=$10^{-2}$, approximates hard-sphere-like interactions. We select $50$ bond length and angle variants, and for each we find the lowest energy dihedral angle conformation and corresponding $P_{i}(\chi_{1},....,\chi_{n})$ values. We average $P_{i}$ over the variants to obtain $P_{m}(\chi_{1},....,\chi_{n})$. We then compare the particular dihedral angle combination $\{\chi_{1}^{HS},...,\chi_{n}^{HS}\}$ associated with the highest value of $P_{m}$ to the side chain of the crystal structure $\{\chi_{1}^{xtal},...,\chi_{n}^{xtal}\}$. To assess the accuracy of the hard-sphere model in predicting the side chain dihedral angles of residues in protein cores, we calculate

[TABLE]

We determine the fraction $F(\Delta\chi)$ of residues of each type with $\Delta\chi$ less than 10o, 20o, and 30o. (See Fig. 6.)

In Fig. 6 (left), we investigate the accuracy of the hard-sphere model in predicting the side chain dihedral angles of single residues in protein cores. For each amino acid (Ile, Leu, Met, Phe, Ser, Thr, Trp, Tyr, and Val), we calculate the fraction of residues, $F(\Delta\chi)$, for which the predicted side chain dihedral angle conformation is within 10o, 20o and 30o of the crystal structure value. Consistent with our prior results, the hard-sphere model accurately predicts the side chain dihedral angle combinations of single residues in the context of the protein for Ile, Leu, Phe, Thr, Trp, Tyr, and Val ($\geq$ 90% within 30o) [49]. This result emphasizes that the purely repulsive hard-sphere model can accurately predict the side chain dihedral angle combinations for nonpolar and uncharged amino acids.

We find that the hard-sphere model is unable to predict with high accuracy the observed side chain conformations for two residues that we studied: Ser and Met. Our results for Met are consistent with those found in Virrueta et al. [55]. In this prior work, we found that local steric interactions were insufficient to predict the shape of the $P(\chi_{3})$ distribution for Met. It was necessary to add attractive atomic interactions to the hard-sphere model to reproduce the observed $P(\chi_{3})$. Here, using only repulsive interactions, we predict $\approx$ 80% of Met residues are within 30o of the crystal structure. Our results for Ser (only 38% within 30o) are also consistent with our prior work in Caballero et al. [49]. We speculate that because the side chain of Ser is small, hydrogen-bonding interactions must be included to correctly place its side chain. In contrast, we suggest that the more bulky Thr and Tyr side chains cause steric interactions to determine the positioning of their side chains, even though they are able to form hydrogen bonds [37].

In addition to single residue rotations, we performed core repacking using combined rotations of interacting core residues in each protein [56]. For the combined rotation method, all residues in an interacting cluster are rotated simultaneously (with fixed backbone conformations), and the global minimum energy conformation is identified. A cluster of interacting residues is defined such that side chain atoms of each residue in the cluster interact with one or more other residues in the cluster, but do not interact with the side chains of other core residues in the protein.

Single and combined rotations have the same prediction accuracy (Figs. 6 and 7), which shows that there are very few arrangements of the residues in a protein core that are sterically allowed and that the side chain conformations of most core residues are dominated by packing constraints. This result implies that there are no alternative sterically allowed conformations of core residues other than those in the crystal structure. If alternative sterically allowed conformations existed, we would have found them using the collective repacking method and thus the prediction accuracy would have dramatically decreased relative to the value for single residue rotations. It does not. Thus, the results for collective repacking reveal that the structures of protein cores are uniquely specified by steric interactions. This conclusion is consistent with those reached by Word et al. [35], where they found that “in a well-packed core region, it is rare that a bond angle can be rotated much in either direction without producing clashes.”

4 Jammed packings of spherical and nonspherical particles

A strict definition of jamming means that a disordered system is solid-like and possesses a static shear modulus [26]. However, jamming also implies that a system is confined to a small region of configuration space, such that little or no motion of the constituent particles can occur. The results presented in Secs. 2 and 3 provide several indications that residues in protein cores are jammed in this latter sense. First, for nearly all protein cores, single and collective repacking give the same side chain dihedral angle combinations found in the protein crystal structures. This result emphasizes that there are no alternative low energy conformations for core residues. Second, the packing fraction of protein cores is $\approx 0.56$, which is similar to those reported for disordered jammed packings of frictional [57] and elongated particles [58, 59, 60].

In this section, we present the results of simulations of jammed packings in three spatial dimensions (3D) for a wide variety of particle shapes including monodisperse spheres, polydisperse spheres, spheres with varying sizes of asperities (or “bumps”), ellipsoids, ellipsoids with varying sizes of asperities, and non-axisymmetric, elongated particles. This range of shapes allows us to study the influence of the particle aspect ratio and surface bumpiness on the packing fraction and determine which particle shapes produce packing fractions that match the packing fraction of residues in protein cores.

We start the discussion with jammed packings of monodisperse spheres. In monodisperse systems, the packing fraction depends on the degree of order that is present in the system. For example, in Fig. 8, we show that the packing fraction varies with the global bond orientational order parameter $Q_{6}$ [61, 62], which measures the degree to which the separation vectors connecting a given particle and its nearest neighbors are consistent with icosohedral symmetry. $Q_{6}\approx 0.57$ for perfect FCC crystalline sphere packings with $\phi\approx 0.74$. The packing fraction for jammed packings of monodisperse spheres decreases as $Q_{6}$ decreases, reaching random close packing $\phi\approx 0.64$ in the limit $Q_{6}\rightarrow 0$ [63]. Jammed packings with different values of $Q_{6}$ can be obtained by varying the rate at which kinetic energy is drained from the system [64]. For the present studies, we consider the limit of fast quenching rates, which gives rise to disordered packings.

Particle size differences can strongly decrease a system’s tendency to order. In previous studies, we focused on jammed packings of bidisperse spheres with half large spheres, half small spheres, and a modest diameter ratio of $d=1.4$ [65, 50]. It is difficult to generate ordered packings of such bidisperse spheres using the packing-generation methods employed here. However, large size ratios ($d\gtrsim 2.4$) can also increase the packing fraction of jammed packings of polydisperse spheres. In this case, small spheres can fill in the gaps between contacting larger spheres. For example, Apollonian sphere packings [66] characterized by a continuous distribution of particle sizes possess packing fractions that approach $1$.

In the hard-sphere model of proteins, we consider seven atom types with differing diameters. The largest diameter ratio is $d=1.8$ between sulfur (which is rare) and hydrogen atoms; the next largest diameter ratio ($d=1.5$) is between sp3 carbon and hydrogen atoms. Thus, we expect that jammed sphere packings composed of mixtures of atoms with the same sizes and number fractions as in protein cores will have packing fraction $\phi\approx 0.64$. This result was shown previously in Ref. [14]. Thus, jammed packings composed of individual spheres with polydispersity that matches atom size differences in protein cores possess packing fractions that are larger than the values we observe in protein cores (Sec. 2).

We now consider jammed packings of symmetric elongated particles, i.e. spherocylinders and ellipsoids, as a function of the aspect ratio $\alpha$. In Fig. 9, we show that the packing fraction $\phi(\alpha)$ is qualitatively the same for jammed packings of spherocylinders and ellipsoids. $\phi\approx 0.64$ for spherical particles with $\alpha=1$, increases for $\alpha>1$, reaches a peak near $\alpha\approx 1.5$ with $\phi>0.7$, and then decreases to a plateau value of $\phi\approx 0.68$ at large $\alpha$.

To compare the results for jammed packings of symmetric, elongated particles to packings of amino acids presented in Sec. 2, we define a generalized aspect ratio and surface bumpiness to characterize the shape of composite particles made from collections of spheres. We define bumpiness by

[TABLE]

where ${\hat{u}}$ is a unit vector with an origin at the geometric center of the composite particle, the integral is over all orientations of ${\hat{u}}$, ${\vec{R}}({\hat{u}})$ gives the location on the surface of the composite particle along ${\hat{u}}$, and ${\vec{\cal R}}({\hat{u}})$ gives the location on the surface of a reference prolate ellipsoid of revolution along ${\hat{u}}$. The bumpiness $B$ for a given composite particle will depend on the orientation of the reference prolate ellipsoid axis ${\hat{e}}$ and the values of the major $a$ and minor $b$ axes.

To define the aspect ratio $\alpha$ for composite particles, we find the reference prolate ellipsoid of revolution that yields the smallest bumpiness. We first fix the reference ellipsoid axis ${\hat{e}}$ to be in the direction that gives the largest distance between the geometric center and the surface of the composite particle. We then minimize $B({\hat{e}},a,b)$ over $a$ and $b$ at fixed ${\hat{e}}$, and define $\alpha=a/b$ for the optimal values of the major and minor axes of the reference ellipsoid.

Fig. 9 shows the packing fraction versus aspect ratio for Ala, Val, Ile, Leu, Met, and Phe residues in protein cores. As discussed in Sec. 2, most core residues have packing fractions near $0.55$-$0.56$. The aspect ratios of amino acids depend on the amino acid type and their backbone and side chain conformations. The average aspect ratios vary from $\alpha\approx 1.4$ for Val to $\approx 2.3$ for Phe. The error bars in both $\phi$ and $\alpha$ are obtained from the root-mean-square fluctuations over different instances of each residue in protein cores.

The packing fraction $\phi\approx 0.55$-$0.56$ observed for amino acids in protein cores with nominal aspect ratios in the range $1.4\lesssim\alpha\lesssim 2.3$ is not consistent with the packing fraction $\phi\approx 0.7$ obtained for jammed packings of ellipsoids and spherocylinders with aspects ratios in the same range. Thus, elongated, smooth, axisymmetric particles are not sufficient to model packings of amino acids in protein cores.

A method for decreasing the packing efficiency of particle packings is to include frictional forces between particles or add asperities (or “bumps”) to the surface of the particles as shown in Fig. 10 (a) and (b). In prior work, we showed in 2D that we could decrease the packing fraction of bidisperse disks from random close to random loose packing (corresponding to more than a $10$% decrease in packing fraction) by increasing the bumpiness or effective friction coefficient between disks [67]. In Fig. 9, we include results from Ref. [57] showing that the packing fraction of frictional spheres (asterisks) in 3D decreases by a similar percentage from $\phi\approx 0.64$ to $\approx 0.55$ as the static friction coefficient $\mu$ increases from $10^{-4}$ to $10$.

We find similar results for bumpy spheres (green squares) in Fig. 9. Here, the bumpy spheres are composite particles made from twelve spheres arranged on the vertices of an icosohedron. We decrease the ratio $r$ of the size of each sphere to the size of the icosohedron to increase the bumpiness $B$. We show in Fig. 11 that for bumpy spheres formed from an icosohedron, we can generate $0\lesssim B\lesssim 0.15$ (corresponding to $5\gtrsim r\gtrsim 0.63$), which accounts for the decrease in packing fraction of the green squares in Fig. 9 from top to bottom.

As discussed above, amino acids cannot be modeled using spherical shapes with $\alpha\approx 1$ or using elongated, smooth particles. Thus, we performed studies of bumpy ellipsoids with $\alpha>1$ to model packings of amino acids in protein cores. For bumpy ellipsoids, we place spheres on the surface of a reference prolate ellipsoid with specified major and minor axes. Two spheres were placed on the ends of the reference ellipsoid and either $3$ or $4$ spheres were placed at equal angular intervals on the ellipsoid surface at distances along the long axis that divide the long axis into $3$ or $4$ equal segments. Thus, the bumpy ellipsoids we studied were made up of either $8$ or $14$ spheres as shown in Fig. 10 (c) and (d). In Fig. 11, we show that we can study bumpiness values $B\lesssim 0.17$ over a wide range of aspect ratios using this method for constructing bumpy axisymmetric elongated particles.

In Fig. 9, we show the packing fraction for jammed packings of bumpy ellipsoids over a range of bumpiness values for two aspect ratios, $\alpha\approx 1.4$ and $2.25$, which spans the range of aspect ratios calculated for amino acids in protein cores. For both aspect ratios, the packing fraction decreases from the values obtained from packings of smooth elongated particles to $\phi\approx 0.55$ as the bumpiness is increased from $B=0.01$ to $0.17$.

An interesting point to note, as shown in Fig. 11, is that amino acids found in protein cores (e.g. Ala and Phe in Fig. 10 (e) and (f)) possess bumpiness values between $B=0.25$ and $0.3$, whereas bumpy axisymmetric shapes have $B\lesssim 0.17$. Thus, we also studied jammed packings composed of the non-axisymmetric composite particles in Fig. 10 (g) and (h). Five spheres make up the composite particle pictured in panel (g). Three are arranged in a straight line, and the other two spheres are placed in a plane perpendicular to the long axis of the composite particle and at an angular separation of $90^{\circ}$. The composite particle in panel (h) contains $7$ spheres with two spheres each placed at the top and bottom of the particle in planes perpendicular to the long axis and in staggered orientations. The bumpiness and aspect ratio of these non-axisymmetric composite particles is varied by changing the size of the bumps compared to the size of the sphere that circumscribes the composite particle. For these two types of non-axisymmetric particles, we were able to increase the maximum bumpiness to $B\approx 0.4$, which is even larger than that of any of the core amino acids (Fig. 11).

As shown in Fig. 9, the packing fractions for jammed packings of the non-axisymmetric particles in Fig. 10 (g) and (h) (with $B=0.33$ and $0.39$) are $\phi\approx 0.56$. These results show that jammed packings of particles with the same $B$ and $\alpha$ as those found for amino acids yield the same packing fraction as amino acids in protein cores.

5 Mutations in protein cores

Additional insight into the packing efficiency in protein cores can be obtained by examining the results from experimental studies of protein core mutations. Several groups have experimentally investigated the potential plasticity of protein cores by performing mutations, i.e. by changing the identities core amino acids. Lim and Sauer simultaneously mutated several hydrophobic residues in the core of a small protein, and used a genetic screen to identify those that were functional and stable. They found that very few combinations of amino acids other than the wildtype set resulted in a stable, folded protein [24]. The functional new cores were dominated by hydrophobic amino acids and the total side chain volumes were within 10% of the original core volume. Combinations of residues outside of these requirements were nonfunctional. Moreover, stereochemical constraints further restricted the allowed sequence space. For example, although many permutations of core residues can maintain the same total volume and hydrophobicity in the core, they do not result in a protein with the same structure and stability [24]. As a result of hydrophobic, volume, and steric constraints, only 0.3% of 60,000 sequences sampled are fully functional [24, 25]. These observations provide experimental support for the dominance of steric interactions in protein cores. Similar experimental results have been found in other proteins [68, 69, 70, 71].

Liu, et al. investigated how mutations from small to large residues in the core affect protein stability [72]. This work illustrates the difficulty in generalizing the effects of a particular type of mutation at different locations and in different proteins. In this work, three Ala residues in the core of a small protein were mutated, individually, to either Cys, Ile, Leu, Met, Phe, Trp, or Val, and the resulting effect on protein stability was determined. They also solved the crystal structures of several of the mutated proteins. They found that in all cases, to varying degrees, to accommodate the larger amino acid side chain, the backbone moved. Interestingly, at two of the three positions, even with backbone movement, the protein with a larger side-chain was destabilized relative to the protein with the original Ala. However, at one position, even large increases in volume (Ala to Phe or Trp) could be accommodated by backbone movement to give a mutated protein with similar stability to that of the parent protein. Liu, et al. hypothesized that this behavior was due to a cavity in the protein near the mutation site, which allowed for more flexibility in this region of the protein [72]. (See also Sec. 6.)

This work shows that the protein core is not able to accommodate mutations to larger residues without significant rearrangement and subsequent destabilization of the original structure. If substantial empty space existed in the protein core, then mutations of this type would likely have small effects because they would fill the existing empty space and not require backbone rearrangements. Instead, backbone rearrangements are necessary to accommodate larger amino acids, supporting the idea that protein cores are tightly packed [72]. This example also illustrates that much is still unknown about protein core packing and how it controls protein stability. The current state of knowledge is such that one can predict neither the backbone movements in response to the incorporation of a larger side chain, nor the changes in stability that result from these structural changes.

6 Conclusions and Future Directions

Our computational studies have established that protein cores are comprised of irregularly shaped objects that are packed into disordered jammed arrangements with $\phi\approx 0.56$ [14]. For a given core, there are no alternative arrangements of the same amino acids that are consistent with a well-packed core with no atomic overlaps [49, 56]. It has also been shown, both experimentally and computationally, that there are a small number of combinations of different core residues that can properly fit in and fill a given core, and thus give rise to a stable folded protein [24, 25, 72, 73, 74]. There are also experimental examples in which amino acids in the core are substituted with ones that are either smaller or larger. Often such substitutions result in changes in the backbone positions. With the current state of understanding in the field, it is not possible to reliably predict such movements. For some mutations, the rearranged protein is as stable as the starting protein, for others it is less stable. Again, the state of the art in computational modeling is such that it is not possible to predict either the structure or the stability of the repacked, rearranged protein.

Even dense packing of amino acids in protein cores results in some void space not occupied by amino acids. There has been some analysis of voids in proteins using a range of probe sizes [75, 23]. Various probe sizes are used to identify void connectivity in the protein and to remove small physically irrelevant voids. Obviously, an exceedingly small probe (e.g. radius $\lesssim 0.05$ Å) will identify a large amount of void space, because even the very smallest voids will be counted. Conversely, a large probe (e.g. radius $\gtrsim 1.4$ Å) will identify few, if any, voids. A ‘reasonable’ probe size to use seems to be around $0.5$ Å. Using such a probe size, Cuff, et al. examined void statistics in a dataset of high-resolution protein structures [75]. They found that the median total void volume was $\approx 15\mathrm{\AA}^{3}$ per residue. To put this into perspective, a CH2 group and a water molecule have a volume of $\approx 25\mathrm{\AA}^{3}$, which indicates that the voids in protein cores are small. In future studies, we will consider the location and size of buried voids to predict the consequences of changes of amino acid size and sequence in protein cores.

7 Acknowledgments

We gratefully acknowledge the support of the Raymond and Beverly Sackler Institute for Biological, Physical, and Engineering Sciences (to L.R., C.S.O., and J.C.G.), National Library of Medicine Training Grant No. T15LM00705628 (to J.C.G.), and National Science Foundation (NSF-PHY-1522467 to L.R., C.S.O. and J.C.G.). This work also benefited from the facilities and staff of the Yale University Faculty of Arts and Sciences High Performance Computing Center and the National Science Foundation (Grant No. CNS-0821132), which in part funded acquisition of the computational facilities. We also thank Mark D. Shattuck for his input on measurements of void volumes in protein cores.

8 Appendix

In this Appendix, we provide additional details that support the results presented in the main text. In Table 2, we provide the volume of the $11$ residues that occur most frequently in protein cores using the explicit hydrogen representation. Gly and Ala have the smallest volumes and Tyr and Trp have the largest. These values differ quantitatively from those obtained using the extended atom model.

We next describe the calculation of the error bars for the fraction $F(\Delta\chi)$ of residues for which the prediction of the hard-sphere model is less than $\Delta\chi$ from the observed side chain conformation that are shown in Figs. 6 and  7. To assess the accuracy of the hard-sphere model in predicting the side chain dihedral angle conformations of residues in protein cores, repacking calculations were performed using $N_{v}=300$ bond length and angle variants for each core residue. For each residue, we randomly select $M$ bond length and angle variants out of the $N_{v}$ variants. For each set of variants, we identified the optimal side chain dihedral angle combination and calculated $\Delta\chi$. We then repeat this process $N$ times, which yields a set of $N\Delta\chi$ values. We then calculated the mean fraction of residues $F(\Delta\chi)$, which satisfy $\Delta\chi<10^{\circ}$, $20^{\circ}$, or $30^{\circ}$, and the standard deviation. We used $N=50$ and $M=50$ for single residue rotations and $N=50$ and $M=30$ for combined rotations.

To understand how particle elongation and surface bumpiness affect packing properties, we generated jammed packings of composite particles formed from spheres. Each composite particle is composed of $n$ spherical asperities placed on the vertices of an icosohedron or locations on the surface of a prolate ellipsoid of revolution. Spherical asperities $i$ and $j$ on composite particles $C$ and $C$’ interact via the pairwise potential $U^{CC^{\prime}}_{ij}=\frac{\epsilon}{2}(1-r_{ij}/\sigma_{ij})^{2}\Theta(\sigma_{ij}-r_{ij})$, where $\epsilon$ is the energy scale of the interaction, $r_{ij}$ is the distance between the centers of asperities $i$ and $j$, $\sigma_{ij}=(\sigma_{i}+\sigma_{j})/2$ is the average diameter of asperities $i$ and $j$, and $\Theta$ is the Heaviside step function. Thus, composite particles $C$ and $C$’ interact via $U^{CC^{\prime}}=\sum_{i,j}U^{CC^{\prime}}_{ij}$. The total potential energy of the system is $U=\sum_{C>C^{\prime}}U^{CC^{\prime}}$.

To find jammed packings, we employ a packing-generation protocol similar to that in Ref. [60]. We first place $N$ composite particles randomly in a cubic periodic cell of unit size. At each step we increase the asperity sizes $\sigma_{i}$ and bond lengths $\delta_{ij}$ between asperities (fixing the ratios between $\sigma_{i}$ and $\delta_{ij}$) corresponding to $\Delta\phi\approx 10^{-3}$, then we relax the system to the nearest potential energy minimum using dissipative dynamics, where the dissipative forces are proportional to the composite particle velocities. If the potential energy is zero after energy minimization (i.e. below a small threshold $U/N<10^{-4}$), we continue compressing; otherwise, we decompress the system, where $\Delta\phi$ is halved each time we switch from compression to decompression. We stop the packing-generation protocol when the potential energy is nonzero and the average particle overlaps are between $0.01$% and $0.1$%. We measure the final packing fraction at jamming onset, which is insensitive to the choice of $\Delta\phi$ and the overlap threshold, provided they are sufficiently small.

9 References