Geometric constraints in protein folding

The intricate three-dimensional geometries of protein tertiary structures underlie protein function and emerge through a folding process from one-dimensional chains of amino acids. The exact spatial sequence and configuration of amino acids, the biochemical environment and the temporal sequence of distinct interactions yield a complex folding process that cannot yet be easily tracked for all proteins. To gain qualitative insights into the fundamental mechanisms behind the folding dynamics and generic features of the folded structure, we propose a simple model of structure formation that takes into account only fundamental geometric constraints and otherwise assumes randomly paired connections. We find that despite its simplicity, the model results in a network ensemble consistent with key overall features of the ensemble of Protein Residue Networks we obtained from more than 1000 biological protein geometries as available through the Protein Data Base. Specifically, the distribution of the number of interaction neighbors a unit (amino acid) has, the scaling of the structure's spatial extent with chain length, the eigenvalue spectrum and the scaling of the smallest relaxation time with chain length are all consistent between model and real proteins. These results indicate that geometric constraints alone may already account for a number of generic features of protein tertiary structures.


Introduction
Proteins consist of sequences of amino acids. The resulting primary structure of a protein, is expected to provide constraints for the folded three-dimensional (3D) structure of a globular protein, its tertiary structure. The problem of predicting the 3D structure of an amino acid sequence in an aqueous solution is known as the protein folding problem consisting of three sub-problems: First, to find the chemically active folded state; second, to uncover the pathway to get to that state; and third, to develop computational tools capable of accurately predicting the folded state [10,18,20,[31][32][33]. Many different avenues have been taken to explore solutions towards this problem, ranging from atomistic models using molecular dynamics approaches [36], to coarse grained models e.g [6], and to machine learning-based and heuristic physical models that disregard the atomistic details of the amino acid sequence [5,9]. While much progress has been made improving molecular dynamics simulations using atomistic detail, the folding process of long chains is computationally highly expensive or even infeasible, and still requires access to purpose build massively parallel computers such as Anton [34], or distributed computing projects such as folding@home in order to generate quantitative data [35]. The other avenue often explored for structure models is tested in community-wide challenges such as the 'Critical Assessment of Protein Structure Prediction' (CASP) [24][25][26]. CASP is run every other year to see if a protein's tertiary structure can be predicted based on its primary sequence of protein structures unresolved at the time of the challenge [Ogorzalek et al.]. Predictions have improved drastically over previous CASP challenges [10], however, often rely on existing structural information in the protein data base (PDB) and homology modelling, comparing new proteins based on existing insights from known template proteins using computational models such as HHPred [4] or I-TASSER [41]. These approaches support accurate prediction of 3D structures, yet by construction limit insights into fundamental physical mechanisms and constraints underlying the folding processes and final structures observed in the many and various proteins observed in nature.
Here, we propose a complementary approach to further understand geometry and formation processes of 3D tertiary structures from chain-like primary protein structures without comparing to specifically chosen protein structures available on the PDB, and without using complex molecular dynamics simulations. First, we analyze 1122 protein structures from the PDB, consider them as an ensemble of protein tertiary structures, and quantify overall properties of this ensemble. In particular we (i) uncover the scaling of the diameter of proteins with their chain length, (ii) reveal the distribution of the number of other amino acids any given amino acid closely interacts with and (iii) find the distribution of second largest eigenvalues of their associated graph Laplacians, characterizing the most persistent time scales on which proteins are dynamically responding to perturbations. Second, we propose and analyze a simple stochastic process modeling the folding of chains of units. The minimal model takes into account geometric constraints only and does not consider any other protein property. The model process keeps connected units connected, forbids geometric overlap of units (volume exclusion) and connects randomly chosen units if geometrically permitted. Based only on such random monomer interactions and geometric constraints, akin to those in Lennard-Jones clusters and sticky hard spheres [37,38], the 3D structures self-organizing through the simple model process are consistent with those of real protein ensembles in all of the above-mentioned features simultaneously.
These results suggest that beyond the details of pairwise interaction of amino acids, from intermediate scales of a few amino acids to the full spatial extent of proteins, geometric constraints play an important role in structure formation and strongly impact the final protein tertiary structure. Our insights may put into perspective the influence of the specific details of sequences of amino acids relative to simpler geometric constraints on structure forming processes of proteins.

Ensemble analysis of Protein Residue Networks
With their modular polymer structure and their complex interaction patterns, proteins lend themselves naturally to a description as ensembles of complex networks. The mathematical object of a graph, simply termed network, represents a structure of nodes (units) and links, each describing an interaction between two units [1,3,27]. Networks and graphs have been used to describe the structure of a wide variety of systems, as different as social networks [15,40] and the global climate system [12,22]. In this article, we analyze an ensemble of 1122 protein tertiary structures of chain lengths ranging from = 8 to = 1500 amino acids. Detailed structures have been experimentally determined to great accuracy and stored in the protein data bank (PDB) [2]. Part of the information stored in the PDB are the coordinates ∈ ℝ 3 of the individual amino acid's central carbon atoms , where indexes the amino acid's position along the chain.
Given such geometric data, the structures resulting from protein folding are commonly expressed as protein residue networks (PRN's) [8,11,13,39], in which the central carbon atom of each amino acid is taken to be a node and a link represents the interaction of two nodes if their spatial distance is small, i.e. less than a distance apart.
Here, the distance between the amino acids indexed and is given by the Euclidean distance metric , = ‖ − ‖. An adjacency matrix encodes the topology of a network, its entries are 1 if , ≤ , i.e. the units are considered connected, and 0 otherwise. The distance matrix resulting from PDB data thus defines the adjacency matrix as The threshold of the PRN is commonly chosen between = 4 Å, the typical length of a covalent bond, and = 15 Å, reflecting an upper bound for a significant interaction to occur between two units. Here, we created the PRNs of 1122 proteins selected from the PDB list in [17], covering a range of chain lengths for comparison to simulations. Their geometric structures have been determined previously via NMR and x-ray studies. We choose a threshold value of = 6.5 Å to calibrate the average degree 1 of nodes in the PRNs to the average degree found in the model simulations in the range of large ∈ [200, 400], Fig. 1a. The average degree ⟨ ⟩ grows with and appears to saturate at a value determined by . The ratio of this cutoff threshold and the unit size in the model, which we take half their mean distance, constitutes the only free structural parameter we employ in the current study. The degree distribution of the resulting network ensemble, displayed in Fig. 1b, is unimodal and covers effective degrees between = 2 and = 11. Interestingly, the degree distribution resulting from simulations of the model ensemble we are about to introduce below is statistically indistinguishable from those of the network of real PRNs (no additional fit parameter), Fig. 1b. Equally, other quantifiers obtained from the simple, geometry-only model ensemble agree surprisingly well with those obtained from our data analysis of the experimentally obtained protein structures.

Simple model focused on geometric constraints
To better understand the impact of geometric constraints on the topology of protein tertiary structures, we introduce a random network formation model that takes into account geometric constraints and leaves out almost all other properties of real proteins, including heterogeneous sequences of amino acids, the amino acids' specifics molecular properties, different forms of electrochemical interactions, conformational details of interactions between nearby amino-acids, and the influence of the fluid environment on protein folding. We find that the simple, geometry-centered model already reproduces a range of overall topological properties of real protein residue networks well.
The model is built on the simple observation that proteins consist of a chain of close-to identical units that interact in complex ways when folding, yet can not intersect, giving rise to geometric constraints. The individual units of the chain interact when they come into contact; typically there is an attraction that is the The average degree of real protein ensemble (red dots) asymptotically saturates to ⟨ ⟩ ≈ 6.8 as the chain length becomes large. The average degree of the nodes resulting from 30 model simulations for each chain length , ranging from = 3 to = 398. b) The degree distribution of the model simulation within the error margin is indistinguishable from that of real proteins (error bars indicate standard deviation of the distribution at each ). stronger the closer they are but repelling once they overlap. Depending on the specific amino acid, size, shape, and electromagnetic properties vary. In our model, however, all amino acids are represented as unit spheres and the interactions between each pair become very simple and identical across all pairs.
The model's initial state consists of a chain of connected spheres, each of diameter and bond length of unity (later rescaled to match the mean distance between neighboring amino acids mean ). A folding proceeds by sequentially picking random pairs of spheres (not connected with each other) and connecting them if possible, given the geometric constraints of volume exclusion. Here, volume exclusion also applies to co-moving other spheres connected either initially along the chain or through a previous steps (see methods section for details). The process repeats until all pairs are either connected or geometrically incapable of connecting. The adjacency matrix sim of the simulated chain keeps track of which spheres are linked to each other. Initially, it contains only zeros except for its secondary diagonal elements which equal 1 since neighboring spheres are connected via the backbone chain. The model is motivated by a two-dimensional model of network-based formation of aggregates where link constraints due to geometry in space have been approximately mapped to purely graph-theoretic constraints during network formation [23].
As described in the method section, the process of moving spheres towards each other is realized in a simple consistent way to satisfy all geometric constraints continuously in time. The forces and potentials employed, however, are not intended to reflect any physical forces or potentials created by amino acids. They plainly help to realize to attempt the joining of two randomly selected spheres.
Snapshots of the folding process are illustrated in Fig. 2, three examples of the final aggregates in Fig. 3. The aggregates are highly compact compared to the straight initial conditions. They are also much more compact than aggregates generated from self-avoiding random walks and close to, yet not quite maximally densely packed (see below), consistent with previous suggestions based on 2D aggregates [23].

Spatial scaling of protein structures
The ensemble of protein tertiary structures exhibits an algebraic scaling law indicating that their radii of gyration depend on their chain length such that: as expected from a number of previous studies [7,17,21,23]. As the overall geometry of a folded protein is often characterized by the locations of the central carbon atoms ( -atoms, one for each amino acid) of its backbone chain, its spatial extension is commonly measured by the radius of gyration quantifying the average distance of units from the center of mass̄ , where is the location of unit ∈ {1, … , }. Our previous study [23] revealed that the scaling law indeed is algebraic and that the exponent is (slightly) larger than for space filling aggregates (where SF = 1 3 = 0.3333 … in 3D) yet (far) smaller than for aggregates created through a self-avoiding random walk (where RW = 3 5 = 0.6 in 3D). That study found = 0.3916 ± 0.0008 for 37162 proteins. For our smaller data set of 1122 proteins, we find = 0.374 ± 0.03, see Fig. 4 for illustration.
To compare the spatial extent of model aggregates, i.e. graph-theoretically defined networks of spheres, to biological proteins on the same footing, we first study how the network diameter compares to the radius of gyration defined through Eq. (3). The graph diameter is defined as the maximum number of links to be taken on the shortest link sequence (also referred to as shortest simple paths) between any pair of units in the PRN. We find that is strongly linearly correlated with the spatial extent of the PRN, Fig. 4. Both the ensemble of biological proteins and the model ensembles studied exhibit a roughly proportional dependence of + 1 on , with the slope obtained from the model data ( = 0.777 Å −1 ) being lower and more precisely determined than that obtained for the PRNs ( = 0.942 Å −1 ). As proportionality factors do not affect the scaling, we thus also find  for both the PDB proteins and geometric-constraint model.
With the cutoff distance for the creation of networks chosen to be = 6.5 Å the resulting average link length in the biological proteins becomes ≈ 5.066 Å, which in Fig. 4 we substituted for the unit length of our model simulations. In the PRNs the network diameters are more dispersed. The lower bound of the experimental data fits well with the simulated structures, suggesting geometric constraints as a major driving mechanism influencing the spatial density during network formation.
Both ensembles show power-law scaling of the diameter. The exponent of = 0.345 ± 0.01 of the simulation is very close to the value of = 0.374 ± 0.03, measured in the PDB data. The plots are shown in Fig. 5. Simulations for heterogeneous systems where the radii of individual units are drawn randomly from the uniform distribution on [1 − , 1 + ] for ∈ {0.0, 0.1, 0.2, 0.3, 0.4, 0.5} increased the variance of the measurements for the radius of gyration, as expected. We did not observe any significant bias in the averages such that the scaling relations stay the same also for heterogeneous systems. The simulated results are found to align very well with the lower bound of folded protein diameters, suggesting that much of the discrepancy (constant factor shifting the measured results up in Fig. 5) can be explained by the fact that the simulation only ceases to make new links when this is no longer geometrically possible. In real proteins on the other hand interactions range from Van-der-Waals interactions to hydrogen bonds and individual monomers vary in size and chemical properties and are subject to thermodynamic fluctuations. All this leads to larger gaps within the folded molecule and hence larger diameters of the PRN's.

Distribution and scaling of Laplacian eigenvalues
Lastly we explore the scaling of the second largest eigenvalue of the graph Laplacian with in Fig. 6 and find that it grows with , approaching a saturation point of ≈ 15 for large .
As two additional features roughly characterizing the dynamic properties of protein residue networks, we consider the distribution and scaling of Laplacian eigenvalues. The Laplacian of a network captures both its  interaction topology and its relaxation and vibration properties [14,29]. If the PRN were made only of the central atoms, the Laplacian would exactly quantify the networks vibrational and relaxational modes. As real PRNs are more complex, the Laplacian spectrum can be taken as a proxy for oscillatory and relaxation dynamics.
Because the eigenvalue spectra intrinsically scale with graph size (here: chain length), we have evaluated the spectra of simulated structures and PRNs of lengths of = 400 ± 30. Fig. 6a shows the histogram of eigenvalues for the 18 PRNs (red) in that length range, accumulating all eigenvalues for each of the 18 PRNs. For comparison, we computed 28 simulated structures (black), that fall in the same length range.
Both eigenvalue spectra exhibit a characteristic unimodal shape. The simulated structures have a more symmetric, slightly broader spectrum with a peak at ≈ 7, while the PRN's have a slightly sharper peak at ≈ 8 and higher probabilities for very small eigenvalues. Similarly, the second largest Laplacian eigenvalue exhibits the same qualitative scaling with chain length for PRNs and geometric-constraint model. The second largest eigenvalue of a network's Laplacian quantifies the time scale of its slowest relaxing mode; as such, its scaling with chain length indicates how intrinsic relaxation time scales change due to the aggregates becoming larger.
The spectra and equally the scaling of the second largest eigenvalues are not indistinguishable between model and biological protein data yet overall exhibit similar properties. Whether or not spectra of model ensemble and PRN ensemble actually agree or disagree cannot be concluded without doubt from the data available, both because at (exactly) fixed chain length there typically is no, one, or only very few proteins available in the real protein data set and because the model realizations at fixed yield very similar spectra due to chain homogeneity. There is no unbiased way we know of to account for uncertainties in and simultaneously inhomogeneities in the chain units such that a unambiguous conclusion can be drawn.

Discussion
In this article we have proposed a simple model of spatial network formation taking into account geometric constraints only. Decoupling the constraints, that drive the folding process (geometry, sequence and solution) and focusing on the geometry allows us insights into the folding mechanisms behind the ensemble features. While this approach does not yield direct predictive power to find the native state of a specific sequence it may narrow down the landscape of possibilities.
We find that geometrically constrained random linking already leads to strong similarities of the resulting structures with protein residue networks in biology. Generalizing a 2D model of purely graph-theoretical network formation presented in [23] to 3D, the model is based upon random link additions with geometric constraints. As the topological shortcut is no longer possible, the geometric constraints are simulated directly. The simulation results were then compared to protein residue networks (PRN's), choosing the threshold such that the mean degrees of simulation results and PRN's matched. As a result, the degree distributions are within the error margins of each other.
The network diameter is linearly related to the radius of gyration in both simulation and data and matches when the simulation results are correctly scaled with the mean connection lengths. The network diameter scales with the chain length as a shifted power law with an exponent of = 0.345 ± 0.01, which is in agreement with value of = 0.374 ± 0.03, measured in PRN's. As in 2D, this is close to space filling. Furthermore, we have studied the Laplacian eigenvalue spectrum and the scaling of the second largest eigenvalue with system size, finding that the two systems are compatible. Using the findings from [14,29] we can infer that the structure of vibrational modes and relaxation properties produced by the model are similar to those found in biological proteins.
These results can be taken as an indication that geometric constraints may be a mechanism behind the scaling behaviour of real protein structures, generating an ensemble also compatible on degree distribution and Laplacian spectrum. Further research, however, is necessary to determine how far the structural similarity reaches. For example by comparing further topological characteristics of PRN's vs. model simu- lations. If the analogy persists, the model could be extended to allow simple sequence features, such as hydrophobicity to attempt to get a simpler predictive model. This may give insights into the folding process, that are otherwise lost in simulation complexity.
Taken together, the above results indicate that coarse ensemble properties of protein tertiary structures are already induced by geometric constraints alone such that only finer scales of the folded structures of individual proteins may be controlled by the details of their amino acid sequences. Such simple models provide a new angle of analyzing protein structures at the coarse scale of ensembles and may help understand core mechanisms underlying the complex folding processes.

Simulation method of the geometric constraint protein model
We have simulated the process modifying the chain geometry in 3D and tested the geometric constraints according to an algorithm consisting of repeated cycles of: 1. A pair ( * , * ) of non-adjacent spheres is randomly chosen from the uniform distribution among the set of untried pairs. 2. The two spheres are attempted to be connected by switching on a force of unit strength pointing towards each other (see Fig. 2), under the geometric constraints: (i) the backbone spheres stay together (ii) no spheres overlap (iii) spheres connected previously stay together. 3. If the selected spheres touch, a new link between them forms and we update the adjacency matrix by setting the elements sim * * = sim * * = 1. Alternatively, if the spheres move less than a velocity threshold Δ ∕Δ , the link is discarded and marked geometrically impossible (see below for details).
and the constraint forces that are gradients of summed potentials quadratic in the distances , = ‖ − ‖. Here the Heaviside step function is defined as Θ( ) = 0 if < 0 and Θ( ) = 1 if ≥ 0. The first term in the final parenthesis in (6) ensures keeping neighboring units along the chain in contact as well as all other pairs of spheres linked so far during the process. The second term causes overlapping spheres to repel each other. is an elastic constant chosen large enough for the constraints to be virtually fulfilled and the final chain statistics being invariant of choosing larger values for , but small enough in order not to limit the allowed numerical time steps unnecessarily. The value = 50 has turned out to meet these conditions. The initial configuration of the chain was drawn from a Boltzmann distribution with probability = −1 exp(− Bend ∕ ) with = 1 and energy where is a normalization constant and is the bending angle at the th unit of the chain, defined as the angle between the adjacent tangential vectors through the scalar (dot) product cos = ( − −1 ) ⋅ ( +1 − ), noting that the sphere diameter equals unity. Initially, generated chains were rejected if any constraint was violated. The prefactor can be interpreted as bending stiffness and determines the persistence length of the initial chains. It was set to = 5 such that initial chains are slightly bent (See Fig. 2 for an example). During a cycle started by selecting the spheres * and * to be pulled together, we monitored their decreasing distance * , * . As soon as * , * ≤ 1, the cycle is considered successful and a new link is formed. We have also periodically checked at intervals Δ whether * , * has shrunk by less than a threshold distance Δ = Δ × × 2∕( ∕2). If this is the case, the cycle is discarded as unsuccessful, because the pair of units cannot make contact due to geometric constraints. The configuration at the beginning of that cycle is then restored. The last factor in Δ is the relative velocity of the spheres * and * in case both -in order to move -have to drag half the other spheres ( ∕2) along. This lower velocity threshold was further decreased by introducing the factor = 0.3 because the final chain statistics weakly varied for larger values but remained the same for smaller values. We have found Δ = 0.15 to be small enough in order not to waste computational time on unsuccessful cycles, but large enough to not abort cycles in which * , * shrinks slowly only temporarily.
The excluded volume forces are nonzero only for pairs of spheres whose distance is less than one. To speed up the simulation, they were only evaluated for spheres that are elements of each others neighbor list listing all spheres within a distance 1 + . We initially generated these lists, then integrated the maximum velocity of all spheres over time and updated the neighbor lists whenever the resulting value exceeded . The value = 0.2 provided the best speed-up. At each integer multiple of 100 cycles, all untried links to a sphere with ∑ =1 sim =∶ ≥ 12 were discarded. This measure was taken to accelerate the simulations as further bonding trials including this sphere are geometrically impossible.

Protein Database protein structure preparation
From reference [17] we obtained the list of PDB files used for their analysis. We split the list into NMR structures and X-ray crystal structures, as the NMR structures would contain multiple protein configurations in their PDB entry. Each PDB was then processed with a custom python script that would count the number of C-atoms found in the structure and order the PDB IDs according to the length of the protein chain. Then from this ordered list every 10th protein was picked, ensuring a good spread of length distributions, a good sample size while also keeping computations easily doable on a workstation. For NMR structures the first structure in the PDB entry was chosen. C-coordinates were then extracted for each protein using MDAnalysis [19,30], from which protein residue contact networks were computed using a cutoff distance of = 6.5 Å. The adjacency matrix PDB was populated according to equation 1. This allows the comparison of the computationally generated adjacency matrix to the PRN generated one. For the network measures and manipulations NetworkX [16] was used.
All simulation details, including the code for reproducing the geometric constraint simulations, as well as the preparation and analysis of PDB files can be found in the following github repository: https: //github.com/ppxasjsm/Geometric-constraints-protein-folding