Research Article

Green function of correlated genes in a minimal mechanical model of protein evolution

See allHide authors and affiliations

PNAS May 15, 2018 115 (20) E4559-E4568; first published April 30, 2018; https://doi.org/10.1073/pnas.1716215115

  1. Contributed by Albert Libchaber, November 20, 2017 (sent for review September 25, 2017; reviewed by Mukund Thattai and Massimo Vergassola)

Significance

Many protein functions involve large-scale motion of their amino acids, while alignment of their sequences shows long-range correlations. This has motivated search for physical links between genetic and phenotypic collective behaviors. The major challenge is the complex nature of protein: nonrandom heteropolymers made of 20 species of amino acids that fold into a strongly coupled network. In light of this complexity, simplified models are useful. Our model describes protein in terms of the Green function, which directly links the gene to force propagation and collective dynamics in the protein. This allows for derivation of basic determinants of evolution, such as fitness landscape and epistasis, which are often hard to calculate.

Abstract

The function of proteins arises from cooperative interactions and rearrangements of their amino acids, which exhibit large-scale dynamical modes. Long-range correlations have also been revealed in protein sequences, and this has motivated the search for physical links between the observed genetic and dynamic cooperativity. We outline here a simplified theory of protein, which relates sequence correlations to physical interactions and to the emergence of mechanical function. Our protein is modeled as a strongly coupled amino acid network with interactions and motions that are captured by the mechanical propagator, the Green function. The propagator describes how the gene determines the connectivity of the amino acids and thereby, the transmission of forces. Mutations introduce localized perturbations to the propagator that scatter the force field. The emergence of function is manifested by a topological transition when a band of such perturbations divides the protein into subdomains. We find that epistasis—the interaction among mutations in the gene—is related to the nonlinearity of the Green function, which can be interpreted as a sum over multiple scattering paths. We apply this mechanical framework to simulations of protein evolution and observe long-range epistasis, which facilitates collective functional modes.

A common physical basis for the diverse biological functions of proteins is the emergence of collective patterns of forces and coordinated displacements of their amino acids (1⇓⇓⇓⇓⇓⇓⇓⇓⇓⇓⇓–13). In particular, the mechanisms of allostery (14⇓⇓⇓–18) and induced fit (19) often involve global conformational changes by hinge-like rotations, twists, or shear-like sliding of protein subdomains (20⇓–22). An approach to examine the link between function and motion is to model proteins as elastic networks (23⇓⇓–26). Decomposing the dynamics of the network into normal modes revealed that low-frequency “soft” modes capture functionally relevant large-scale motion (27⇓⇓–30), especially in allosteric proteins (31⇓–33). Recent works associate these soft modes with the emergence of weakly connected regions in the protein (Fig. 1 A and B)—“cracks,” “shear bands,” or “channels” (21, 22, 34⇓–36)—that enable viscoelastic motion (37, 38). Such patterns of “floppy” modes (39⇓⇓–42) emerge in models of allosteric proteins (36, 43⇓–45) and networks (46⇓–48).

Fig. 1.
Fig. 1.

Protein as an evolving machine and propagation of mechanical forces. (A) Formation of a softer shear band (red) separating the protein into two rigid subdomains (light blue). When a ligand binds, the biochemical function involves a low-energy hinge-like or shear motion (arrows). (B) Shear band and large-scale motion in a real protein: the arrows show the displacement of all amino acids in human glucokinase when it binds glucose (Protein Data Bank ID codes 1v4s and 1v4t). The coloring shows a high-shear region (red) separating two low-shear domains that move as rigid bodies (shear calculated as in refs. 21 and 36). (C) The mechanical model. The protein is made of two species of amino acids, polar (P; red) and hydrophobic (H; blue), with a sequence that is encoded in a gene. Each amino acid forms weak or strong bonds with its 12 near neighbors (Right) according to the interaction rule in the table (Left). (D) The protein is made of 10×20=200 amino acids with positions that are randomized from a regular triangular lattice. Strong bonds are shown as gray lines. Evolution begins from a random configuration (Left) and evolves by mutating one amino acid at each step, switching between H and P. The fitness is the mechanical response to a localized force probe (pinch) (2). After ∼103 mutations (Center; intermediate stage), the evolution reaches a solution (Right). The green arrows show the mechanical response: a hinge-like, low-energy motion with a shear band starting at the probe and traversing the protein, qualitatively similar to B. L, left; R, right.

Like their dynamic phenotypes, proteins’ genotypes are remarkably collective. When aligned, sequences of protein families show long-range correlations among the amino acids (49⇓⇓⇓⇓⇓⇓⇓⇓⇓⇓⇓–61). The correlations indicate epistasis, the interaction among mutations that takes place among residues linked by physical forces or common function. By inducing nonlinear effects, epistasis shapes the protein’s fitness landscape (62⇓⇓⇓⇓⇓–68). Provided with sufficiently large data, analysis of sequence variation can predict the 3D structure of proteins (50⇓–52), allosteric pathways (53⇓–55), epistatic interactions (56, 57), and coevolving subsets of amino acids (58⇓–60, 69).

Still, the mapping between sequence correlation and collective dynamics—and in particular, the underlying epistasis—is not fully understood. Experiments and simulations provide valuable information on protein dynamics, and extensive sequencing accumulates databases required for reliable analysis; however, there remain inherent challenges: the complexity of the physical interactions and the sparsity of the data. The genotype-to-phenotype map of proteins connects spaces of huge dimension, which are hard to sample, even by high-throughput experiments or natural evolution (70⇓–72). A complementary approach is the application of simplified coarse-grained models, such as lattice proteins (73⇓–75) or elastic networks (24), which allow one to extensively survey the map and examine basic questions of protein evolution. Such models have been recently used to study allosteric proteins (35, 36, 43⇓–45) and in networks (46⇓–48). Our aim here is different: to construct a simplified model of how the collective dynamics of functional proteins directs their evolution and in particular, to give a mechanical interpretation of epistasis.

This paper introduces a coarse-grained theory that treats protein as an evolving amino acid network with topology that is encoded in the gene. Mutations that substitute one amino acid with another tweak the interactions, allowing the network to evolve toward a specific mechanical function: in response to a localized force, the protein will undergo a large-scale conformational change (Fig. 1 C and D). We show that the application of a Green function (76, 77) is a natural way to understand the protein’s collective dynamics. The Green function measures how the protein responds to a localized impulse via propagation of forces and motion. The propagation of mechanical response across the protein defines its fitness and directs the evolutionary search. Thus, the Green function explicitly defines the map: gene → amino acid network → protein dynamics → function. We use this map to examine the effects of mutations and epistasis. A mutation perturbs the Green function and scatters the propagation of force through the protein (Fig. 2). We quantify epistasis in terms of “multiple scattering” pathways. These indirect physical interactions appear as long-range correlations in the coevolving genes.

Fig. 2.
Fig. 2.

Force propagation, mutations, and epistasis. (A) The Green function G measures the propagation of the mechanical signal, depicted as a “diffraction wave,” across the protein (blue) from the force source f (pinch) to the response site v. (B) A mutation δHi deflects the propagation of force. The effect of the mutation on the propagator δG can be described as a series of multiple scattering paths (6). (C) The epistasis between two mutations, δHi and δHj, is equivalent to a series of multiple scattering paths (9).

Using a Metropolis-type evolution algorithm, solutions are quickly found, typically after ∼103 steps. Mutations add localized perturbations to the amino acid network, which are eventually arranged by evolution into a continuous shear band. Protein function is signaled by a topological transition, which occurs when a shearable band of weakly connected amino acids separates the protein into rigid subdomains. The set of solutions is sparse: there is a huge reduction of dimension between the space of genes to the spaces of force and displacement fields. We find a tight correspondence between correlations in the genotype and phenotype. Owing to its mechanical origin, epistasis becomes long ranged along the high-shear region of the channel.

Model: Protein as an Evolving Machine

The Amino Acid Network and Its Green Function.

We use a coarse-grained description in terms of an elastic network (23⇓⇓⇓–27, 39) with connectivity and interactions that are encoded in a gene (Fig. 1 C and D). Similar vector elasticity models were considered in refs. 35 and 36 (app. B3 therein). The protein is a chain of na=200 amino acids: ai (i=1,…,na) folded into a 10×20 2D hexagonal lattice (d=2). We follow the HP model (73, 74) with its two species of amino acids, hydrophobic (ai=H) and polar (ai=P). The amino acid chain is encoded in a gene c, a sequence of 200 binary codons, where ci=1 encodes an H amino acid and ci=0 encodes a P amino acid.

We consider a constant fold, and therefore, any particular codon ci in the gene encodes an amino acid ai at a certain constant position ri in the protein. The positions ri are randomized to make the network amorphous. These nd=d⋅na=400 dfs are stored in a vector r. Except the ones at the boundaries, every amino acid is connected by harmonic springs to z=12 nearest and next nearest neighbors. There are two flavors of bonds according to the chemical interaction, which is defined as an AND gate: a strong H−H bond and weak H−P and P−P bonds. The strength of the bonds determines the mechanical response of the network to a displacement field u, when the amino acids are displaced as ri→ri+ui. The response is captured by Hooke’s law that gives the force field f induced by a displacement field, f=H(c) u . The analogue of the spring constant is the Hamiltonian H(c), a nd×nd matrix, which records the connectivity of the network and the strength of the bonds. H(c) is a nonlinear function of the gene c, reflecting the amino acid interaction rules of Fig. 1C (Eq. 11, Materials and Methods).

Evolution searches for a protein that will respond by a prescribed large-scale motion to a given localized force f (“pinch”). In induced fit, for example, specific binding of a substrate should induce global deformation of an enzyme. The response u is determined by the Green function G (76):u=G(c) f .[1]G is the mechanical propagator that measures the transmission of signals from the force source f across the protein (Fig. 2A). Eq. 1 constitutes an explicit genotype-to-phenotype map from the genotype c to the mechanical phenotype u: c→u(c)=G(c)f. This reflects the dual nature of the Green function G: in the phenotype space, it is the linear mechanical propagator that turns a force into motion, u=G f, whereas it is also the nonlinear function that maps the gene into a propagator, c→G(c).

When the protein is moved as a rigid body, the lengths of the bonds do not change, and the elastic energy cost vanishes. A 2D protein has n0=3 such zero modes (Galilean symmetries), two translations, and one rotation, and H is, therefore, always singular. Hence, Hooke’s law and [1] imply that G is the pseudoinverse of the Hamiltonian, G(c)=H(c)+ (78, 79), which amounts to inversion of H in the nonsingular subspace of the nd−n0=397 nonzero modes (Materials and Methods). A related quantity is the resolvent, G(ω)=(ω−H)−1, with poles at the energy levels of H, ω=λk.

The fitness function rewards strong mechanical response to a localized probe (pinch in Fig. 1D): a force dipole at two neighboring amino acids p′ and q′ on the left side of the protein (L in Fig. 1D), fq′=−fp′. The prescribed motion is specified by a displacement vector v, with a dipolar response, vq=−vp, on the right side of the protein (R in Fig. 1D). The protein is fitter if the pinch f produces a large deformation in the direction specified by v. To this end, we evolve the amino acid network to increase a fitness function F, which is the projection of the displacement u=Gf on the prescribed response v:F(c)=v⊺u=v⊺G(c) f .[2]Eq. 2 defines the fitness landscape F(c). Here, we examine particular examples for a localized pinch f and prescribed response v, which drive the emergence of a hinge-like mode. This approach is general and can as well treat more complex patterns of force and motion.

Evolution Searches in the Mechanical Fitness Landscape.

Our simulations search for a prescribed response v induced by a force f applied at a specific site on the left side (pinch). The prescribed dipolar response may occur at any of the sites on the right side. This gives rise to a wider shear band that allows the protein to perform general mechanical tasks (unlike specific allostery tasks of communicating between specified sites on L and R). We define the fitness as the maximum of F [2] over all potential locations of the channel’s output (typically 8–10 sites) (Materials and Methods). The protein is evolved via a point mutation process where, at each step, we flip a randomly selected codon between zero and one. This corresponds to exchanging H and P at a random position in the protein, thereby changing the bond pattern and the elastic response by softening or stiffening the amino acid network.

Evolution starts from a random protein configuration encoded in a random gene. Typically, we take a small fraction of amino acid of type P (about 5%) randomly positioned within a majority of H (Fig. 1D, Left). The high fraction of strong bonds renders the protein stiff and therefore, of low initial fitness F≃0. At each step, we calculate the change in the Green function δG (by a method explained below) and use it to evaluate from [2] the fitness change δF:δF=v⊺δG f.[3]The fitness change δF determines the fate of the mutation: we accept the mutation if δF≥0; otherwise, the mutation is rejected. Since fitness is measured by the criterion of strong mechanical response, it induces softening of the amino acid network.

The typical evolution trajectory lasts about 103 steps. Most are neutral mutations (δF≃0) and deleterious ones (δF<0); the latter are rejected. About a dozen or so beneficial mutations (δF>0) drive the protein toward the solution (Fig. 3A). The increase in the fitness reflects the gradual formation of the channel, while the jump in the shear signals the emergence of the soft mode. The first few beneficial mutations tend to form weakly bonded P-enriched regions near the pinch site on the left side and close to the right boundary of the protein. The following ones join these regions into a floppy channel (a shear band), which traverses the protein from left to right. We stop the simulation when the fitness reaches a large positive value Fm∼5. The corresponding gene c* encodes the functional protein. The ad hoc value Fm∼5 signals slowing down of the fitness curve toward saturation at F>Fm, as the channel has formed and now only continues to slightly widen. In this regime, even a tiny pinch will easily excite a large-scale motion with a distinct high-shear band (Fig. 1D, Right).

Fig. 3.
Fig. 3.

The mechanical Green function and the emergence of protein function. (A) Progression of the fitness F during the evolution run shown in Fig. 1D (black) together with the fitness trajectory averaged over ∼106 runs ⟨F⟩ (red). Shown are the last 16 beneficial mutations toward the formation of the channel. The contribution of the emergent low-energy mode ⟨F*⟩ (blue) dominates the fitness [4]. (B) Landscape of the fitness change δF [3] averaged over ∼106 solutions for all 200 possible positions of point mutations at a solution. Underneath, the average amino acid configuration of the protein is shown in shades of red (P) and blue (H). In most sites, mutations are neutral, while mutations in the channel are deleterious on average. L, left. (C) The average magnitude of the two-codon correlation |Qij| [5] in the shear band (amino acids in rows 7–13; red) and in the whole protein (black) as a function of the number of beneficial mutations, t. (Inset) Profile of the spatial correlation g(r) within the shear band (after t=1,11,16 beneficial mutations). (D) The mean shear in the protein in a single run (black) and averaged over ∼106 solutions (red) as a function of the fraction of P amino acids, p. The values of p are shifted by the position of the jump, pc. (Inset) Distribution of pc.

Results

Mechanical Function Emerges at a Topological Transition.

The hallmark of evolution reaching a solution gene c* is the emergence of a new zero-energy mode, u*, in addition to the three Galilean symmetry modes. Near the solution, the energy of this mode λ* almost closes the spectral gap, λ*→0, and G(ω) has a pole at ω≈0. As a result, the emergent mode dominates the Green function, G≃u*u*⊺/λ*. The response to a pinch will be mostly through this soft mode, and the fitness [2] will diverge asF(c*)≃F*=v⊺u*u*⊺fλ*.[4]On average, we find that the fitness increases exponentially with the number of beneficial mutations (Fig. 3A). However, beneficial mutations are rare and are separated by long stretches of neutral mutations. This is evident from the fitness landscape (Fig. 3B), which shows that, in most sites, the effect of mutations is practically neutral. The vanishing of the spectral gap, λ*→0, manifests as a topological change in the system: the amino acid network is now divided into two domains that can move independently of each other at low energetic cost. The soft mode appears at a dynamical phase transition, where the average shear in the protein jumps abruptly as the channel is formed and the protein can easily deform in response to the force probe (Fig. 3D).

As the shear band is taking shape, the correlation among codons builds up. To see this, we align genes from the ∼106 simulations in analogy to sequence alignment of real protein families (49⇓⇓⇓⇓⇓⇓⇓⇓⇓⇓⇓–61), albeit without the phylogenetic correlation that hampers the analysis of real sequences. At each time step, we calculate the two-codon correlation Qij between all pairs of codons ci and cj:Qij≡⟨cicj⟩−⟨ci⟩⟨cj⟩,[5]where brackets denote ensemble averages. We find that most of the correlation is concentrated in the forming channel (Fig. 3C), where it is 10-fold larger than in the whole protein. In the channel, there is significant long-range correlation shown in the spatial profile of the correlation g(r) (Fig. 3C, Inset). Analogous regions of coevolving residues appear in real protein families (53⇓–55, 58⇓–60) as well as in coarse-grained models of protein allostery (35, 36, 43, 44) and allosteric networks (46, 47).

Point Mutations Are Localized Mechanical Perturbations.

A mutation may vary the strength of no more than z=12 bonds around the mutated amino acid (Fig. 2B). The corresponding perturbation of the Hamiltonian δH is, therefore, localized, akin to a defect in a crystal (80, 81). The mechanics of mutations can be further explored by examining the perturbed Green function, G′=G+δG, which obeys the Dyson equation (77, 82) (Materials and Methods):G′=G−G δH G′.[6]The latter can be iterated into an infinite seriesδG=G′−G=−G δH G+G δH G δH G−⋯.This series has a straightforward physical interpretation as a sum over multiple scatterings: as a result of the mutation, the elastic force field is no longer balanced by the imposed force f, leaving a residual force field δf=δH u=δH G f. The first scattering term in the series balances δf by the deformation δu=G δf=G δH Gf. Similarly, the second scattering term accounts for further deformation induced by δu and so forth. In practice, we calculate the mutated Green function using the Woodbury formula [12], which exploits the localized nature of the perturbation to accelerate the computation by a factor of ∼104 (Materials and Methods).

Epistasis Links Protein Mechanics to Genetic Correlations.

Our model provides a calculable definition of epistasis, the nonlinearity of the fitness effect of interacting mutations (Fig. 2C). We take a functional protein obtained from the evolution algorithm and mutate an amino acid at a site i. This mutation induces a change in the Green function δGi (calculated by [12]) and hence, in the fitness function δFi [3]. One can similarly perform another independent mutation at a site j, producing a second deviation, δGj and δFj. Finally, starting again from the original solution, one mutates both i and j simultaneously, with a combined effect δGi,j and δFi,j. The epistasis eij measures the departure of the double mutation from additivity of two single mutations:eij≡δFi,j−δFi−δFj.[7]To evaluate the average epistatic interaction among amino acids, we perform the double-mutation calculation for all 106 solutions and take the ensemble average Eij=⟨eij⟩. Landscapes of Eij show significant epistasis in the channel (Fig. 4). Amino acids outside the high-shear region show only small epistasis, since mutations in the rigid domains hardly change the elastic response. The epistasis landscapes (Fig. 4 A–C) are mostly positive, since the mutations in the channel interact antagonistically (83): after a strongly deleterious mutation, a second mutation has a smaller effect.

Fig. 4.
Fig. 4.

Mechanical epistasis. The epistasis [7], averaged over ∼106 solutions Eij=⟨eij⟩, between a fixed amino acid at position i (black arrow) and all other positions j. Here, i is located at (A) the binding site, (B) the center of the channel, and (C) slightly off the channel. Underneath, the average amino acid configuration of the protein is drawn in shades of red (P) and blue (H). Significant epistasis mostly occurs along the P-rich channel, where mechanical interactions are long ranged. Although epistasis is predominantly positive, negative values also occur, mostly at the boundary of the channel (C). Landscapes are plotted for specific output site at right. L, left. (D) The two-codon correlation function Qij [5] measures the coupling between mutations at positions i and j [5]. The epistasis Eij and the gene correlation Qij show similar patterns. Axes are the positions of i and j loci. Significant correlations and epistasis occur mostly in and around the channel region (positions ∼70–130, rows 7–13).

Definition [7] is a direct link between epistasis and protein mechanics: the nonlinearity (“curvature”) of the Green function measures the deviation of the mechanical response from additivity of the combined effect of isolated mutations at i and j, ΔGi,j≡δGi,j−δGi−δGj. The epistasis eij is simply the inner product value of this nonlinearity with the pinch and the response:eij=v⊺ ΔGi,j f.[8]Relation [8] shows how epistasis originates from mechanical forces among mutated amino acids.

In the gene, epistatic interactions are manifested in codon correlations (56, 57) shown in Fig. 4D, which depicts two-codon correlations Qij from the alignment of ∼106 functional genes c* [5]. We find a tight correspondence between the mean epistasis Eij=⟨eij⟩ and the codon correlations Qij. Both patterns exhibit strong correlations in the channel region with a period equal to channel’s length: 10 amino acids. The similarity in the patterns of Qij and Eij indicates that a major contribution to the long-range correlations observed among aligned protein sequences stems from the mechanical interactions propagating through the amino acid network.

Epistasis as a Sum over Scattering Paths.

One can classify epistasis according to the interaction range. Neighboring amino acids exhibit contact epistasis (49⇓–51), because two adjacent perturbations, δHi and δHj, interact nonlinearly via the AND gate of the interaction table (Fig. 1C), ΔHi,j≡δHi,j−δHi−δHj≠0 (where δHi,j is the perturbation by both mutations). The leading term in the Dyson series [6] of ΔGi,j is a single scattering from an effective perturbation with an energy ΔHi,j, which yields the epistasiseij=−v⊺G ΔHi,jGf+⋯.Long-range epistasis among nonadjacent, noninteracting perturbations (ΔHi,j=0) is observed along the channel (Fig. 4). In this case, [6] expresses the nonlinearity ΔGi,j as a sum over multiple scattering paths, which include both i and j (Fig. 2C):eij=v⊺G δHiG δHjG+G δHjG δHiGf−⋯.[9]The perturbation expansion directly links long-range epistasis to shear deformation: near the transition, the Green function is dominated by the soft mode, G≃u*u*⊺/λ*, with fitness F given by [4]. From [6] and [8], we find a simple expression for the mechanical epistasis as a function of the shear:eij≃F⋅hi1+hi+hj1+hj−hi+hj1+hi+hj.[10]The factor hi≡u*⊺δHiu*/λ* in [10] is the ratio of the change in the shear energy due to mutation at i (the expectation value of δHi) and the energy λ* of the soft mode, and it is similar for hj. Thus, hi and hj are significant only in and around the shear band, where the bonds varied by the perturbations are deformed by the soft mode. When both sites are outside the channel, hi,hj≪1, the epistasis [10] is small, eij≃2hihjF. It remains negligible even if one of the mutations, i, is in the channel, hj≪1≪hi, and eij≃hjF. Epistasis can only be long ranged along the channel when both mutations are significant, hi,hj≫1, and eij≃F1−hi−1−hj−1+(hi+hj)−1≃F[1−1/min(hi,hj)]. We conclude that epistasis is maximal when both sites are at the start or end of the channel as illustrated in Fig. 4. The nonlinearity of the fitness function gives rise to antagonistic epistasis.

Geometry of Fitness Landscape and Gene-to-Function Map.

With our mechanical evolution algorithm, we can swiftly explore the fitness landscape to examine its geometry. The genotype space is a 200D hypercube with vertices that are all possible genes c. The phenotypes reside in a 400D space of all possible mechanical responses u. The Green function provides the genotype-to-phenotype map [1]. A functional protein is encoded by a gene c* with fitness that exceeds a large threshold, F(c*)≥Fm≃5, and the functional phenotype is dominated by the emergent zero-energy mode, u(c*)≃u* (Fig. 3A). We also characterize the phenotype by the shear field s* (Materials and Methods).

The singular value decomposition (SVD) of the 106 solutions returns a set of eigenvectors with ordered eigenvalues that show their significance in capturing the data (Materials and Methods). The SVD spectra reveal strong correspondence between the genotype c* and the phenotype, u and s* (Fig. 5). In all three datasets, the largest eigenvalues are discrete and stand out from the bulk continuous spectrum. These are the collective dfs, which show loci in the gene and positions in the “flow” (i.e., displacement) and shear fields that tend to vary together.

Fig. 5.
Fig. 5.

From gene to mechanical function: spectra and dimensions. (A) The first four SVD eigenvectors (in the text) of the gene Ck (Top), the displacement flow field Uk (Middle), and the shear Sk (Bottom). (B) Cross-sections through the set of solutions in the genotype space (Upper) and the phenotype space (Lower). Density of solutions is color coded. The genotype cross-section is the plane defined by the eigenvectors C3−C100, and in the phenotype space, it is defined by the eigenvectors U3−U100 (in the text). The dimensional reduction is manifested by the discoid geometry of the phenotype cloud compared with the spheroid shape of the genotype cloud. (C) Genetic correlations Qij show similarity to correlations in the shear field, s* (color coded in log scale). Corr, correlation.

We examine the correspondence among three sets of eigenvectors: {Uk} of the flow, {Ck} of the gene, and {Sk} of the shear. The first eigenvector of the flow, U1, is the hinge motion caused by the pinch, with two eddies rotating in opposite directions (Fig. 5A). The next two modes, shear (U2) and breathing (U3), also occur in real proteins, such as glucokinase (Fig. 1B). The first eigenvectors of the shear S1 and of the gene sequence C1 show that the high-shear region is mirrored as a P-rich region, where a mechanical signal may cause local rearrangement of the amino acids by deforming weak bonds. In the rest of the protein, the H-rich regions move as rigid bodies with insignificant shear. The higher gene eigenvectors, Ck (k>1), capture patterns of correlated genetic variations. The striking similarity between the sequence correlation patterns Ck and the shear eigenvectors Sk shows a tight genotype-to-phenotype map, as is further shown in the likeness of the correlation matrices of the amino acid and shear flow (Fig. 5C).

In the phenotype space, we represent the displacement field u in the SVD basis, {Uk} (Fig. 5B). Since ∼90% of the data are explained by the first ∼15 Uk, we can compress the displacement field without much loss into the first 15 coordinates. This implies that the set of solutions is a 15D discoid, which is flat in most directions. In contrast, representation of the genes c* in the SVD frame of reference (with the {Ck} basis) reveals that, in genotype space, the solution set is an incompressible 200D spheroid (Fig. 5B). The dramatic dimensional reduction in mapping genotypes to phenotypes stems from the different constraints that shape them (36, 84⇓⇓⇓⇓–89).

Discussion

Theories of protein need to combine the many-body physics of the amino acid matter with the evolution of genetic information, which together, give rise to protein function. We introduced a simplified theory of protein, where the mapping between genotype and phenotype takes an explicit form in terms of mechanical propagators (Green functions), which can be efficiently calculated. As a functional phenotype, we take cooperative motion and force transmission through the protein [2]. This allows us to map genetic mutations to mechanical perturbations, which scatter the force field and deflect its propagation [3 and 6] (Fig. 2). The evolutionary process amounts to solving the inverse scattering problem: given prescribed functional modes, one looks for network configurations that yield this low end of the dynamical spectrum. Epistasis, the interaction among loci in the gene, corresponds to a sum over all multiple scattering trajectories or equivalently, the nonlinearity of the Green function [7 and 8]. We find that long-range epistasis signals the emergence of a collective functional mode in the protein. The results of this theory (in particular, the expressions for epistasis) follow from the basic geometry of the amino acid network and the localized mutations and are, therefore, applicable to general tasks and fitness functions with multiple inputs and responses.

Materials and Methods

The Mechanical Model of Protein.

We model the protein as an elastic network made of harmonic springs (23, 24, 39, 90). The connectivity of the network is described by a hexagonal lattice with vertices that are amino acids and edges that correspond to bonds. There are na=10×20=200 amino acids indexed by Roman letters and nb bonds indexed by Greek letters. We use the HP model (73) with two amino acid species, hydrophobic (ai=H) and polar (ai=P). The amino acid chain is encoded in a gene c, where ci=1 encodes H and ci=0 encodes P, i=1,…,na. The degree zi of each amino acid is the number of amino acids to which it is connected by bonds. In our model, most amino acids have the maximal degree, which is z=12, while amino acids at the boundary have fewer neighbors, z<12 (Fig. 1C). The connectivity of the graph is recorded by the adjacency matrix A, where Aij=1 if there is a bond from j to i and Aij=0 otherwise. The gradient operator ∇ relates the spaces of bonds and vertices (and is, therefore, of size nb×na): if vertices i and j are connected by a bond α, then ∇αi=+1 and ∇αj=−1. As in the continuum case, the Laplace operator Δ is the product Δ=∇⊺∇. The nondiagonal elements Δij are −1 if i and j are connected and 0 otherwise. The diagonal part of Δ is the degree Δii=zi. Hence, we can write the Laplacian as Δ=Z−A, where Z is the diagonal matrix of the degrees zi.

We embed the graph in Euclidean space Ed (d=2) by assigning positions ri∈Ed to each amino acid. We concatenate all positions in a vector r of length na⋅d≡nd. Finally, to each bond, we assign a spring with constant kα, which we keep in a diagonal nb×nb matrix K. The strength of the spring is determined by the AND rule of the HP model’s interaction table (Fig. 1C), kα=kw+(ks−kw)cicj, where ci and cj are the codons of the amino acid connected by bond α. This implies that a strong H−H bond has kα=ks, whereas the weak bonds H−P and P−H have kα=kw. We usually take ks=1 and kw=0.01. This determines a spring network. We also assume that the initial configuration is such that all springs are at their equilibrium length, disregarding the possibility of “internal stresses” (39), so that the initial elastic energy is E0=0.

We define the “embedded” gradient operator D (of size nb×nd), which is obtained by taking the graph gradient ∇ and multiplying each nonzero element (±1) by the corresponding direction vector nij=ri−rj/ri−rj. Thus, D is a tensor (D=∇αinij), which we store as a matrix (α is the bond connecting vertices i and j). In each row vector of D, which we denote as mα≡Dα,:, there are only 2d nonzero elements. To calculate the elastic response of the network, we deform it by applying a force field f, which leads to the displacement of each vertex by ui to a new position ri+ui (39). For small displacements, the linear response of the network is given by Hooke’s law, f=Hu. The elastic energy is E=u⊺Hu/2, and the Hamiltonian, H=D⊺KD, is the Hessian of the elastic energy E, Hij=δ2E/(δuiδuj). By rescaling, D→K1/2D, which amounts to scaling all distances by 1/kα, we obtain H=D⊺D. It follows that the Hamiltonian is a function of the gene H(c), which has the structure of the Laplacian Δ multiplied by the tensor product of the direction vectors. Each d×d block Hij (i ̸ = j) is a function of the codons ci and cj:Hij(ci,cj)=Δijnijnij⊺=−Aijkw+ks−kwcicjnijnij⊺.[11]The diagonal blocks complete the row and column sums to zero, Hii=−∑j ≠ iHij.

The Inverse Problem: Green Function and Its Spectrum.

The Green function G is defined by the inverse relation to Hooke’s law, u=Gf [1]. If H were invertible (nonsingular),G would have been just G=H−1. However, H is always singular owing to the zero-energy (Galilean) modes of translation and rotation. Therefore, one needs to define G as the Moore–Penrose pseudoinverse (78, 79), G=H+, on the complement of the space of Galilean transformations. The pseudoinverse can be understood in terms of the spectrum of H. There are at least n0=d(d+1)/2 zero modes: d translation modes and d(d−1)/2 rotation modes. These modes are irrelevant and will be projected out of the calculation (note that these modes do not come from missing connectivity of the graph Δ itself but from its embedding in Ed). H is singular but is still diagonalizable (since it has a basis of dimension nd), and it can be written as the spectral decomposition, H=∑k=1ndλkukuk⊺, where {λk} is the set of eigenvalues and {uk} are the corresponding eigenvectors (note that k denotes the index of the eigenvalue, while i and j denote amino acid positions). For a nonsingular matrix, one may calculate the inverse simply as H−1=∑k=1ndλk−1ukuk⊺. Since H is singular, we leave out the zero modes and get the pseudoinverse H+, G=H+=∑k=n0+1ndλk−1ukuk⊺. It is easy to verify that, if u is orthogonal to the zero modes, then u=GHu. The pseudoinverse obeys the four requirements (78): (i) HGH=H, (ii) GHG=G, (iii) (HG)⊺=HG, and (iv) (GH)⊺=GH. In practice, as the projection commutes with the mutations, the pseudoinverse has most virtues of a proper inverse. The reader might prefer to link G and H through the heat kernel, K(t)=∑keλktukuk⊺. Then, G=∫0∞dt K(t) and H=ddt K|t=0.

Pinching the Network.

A pinch is given as a localized force applied at the boundary of the “protein.” We usually apply the force on a pair of neighboring boundary vertices, p′ and q′. It seems reasonable to apply a force dipole (i.e., two opposing forces fq′=−fp′), since a net force will move the center of mass. This pinch is, therefore, specified by the force vector f (of size nd), with the only 2d nonzero entries being fq′=−fp′. Hence, it has the same structure as a bond vector mα of a “pseudobond” connecting p′ and q′ A normal pinch f has a force dipole directed along the rp′−rq′ line (the np′q′ direction). Such a pinch is expected to induce a hinge motion. A shear pinch will be in a perpendicular direction ⊥np′q′ and is expected to induce a shear motion.

Evolution tunes the spring network to exhibit a low-energy mode, in which the protein is divided into two subdomains moving like rigid bodies. This large-scale mode can be detected by examining the relative motion of two neighboring vertices, p and q, at another location at the boundary (usually at the opposite side). Such a desired response at the other side of the protein is specified by a response vector v, and the only nonzero entries correspond to the directions of the response at p and q. Again, we usually consider a “dipole” response vq=−vp.

Evolution and Mutation.

The quality of the response (i.e., the biological fitness) is specified by how well the response follows the prescribed one v. In the context of our model, we chose the (scalar) observable F as F=v⊺u=vp⋅up+vq⋅uq=v⊺Gf [2]. In an evolution simulation, one would exchange amino acids between H and P, while demanding that the fitness change δF is positive or nonnegative. By this, we mean δF>0 is thanks to a beneficial mutation, whereas δF=0 corresponds to a neutral one. Deleterious mutations δF<0 are generally rejected. A version that accepts mildly deleterious mutations (a finite temperature Metropolis algorithm) gave similar results. We may impose a stricter minimum condition δF≥ε F with a small positive ε, say 1%. An alternative, stricter criterion would be the demand that each of the terms in F, vp⋅up and vq⋅uq, increases separately. The evolution is stopped when F≥Fm∼5, which signals the formation of a shear band. When simulations ensue beyond Fm∼5, the band slightly widens, and the fitness slows down and converges at a maximal value, typically Fmax∼8.

Evolving the Green Function Using the Dyson and Woodbury Formulas.

The Dyson formula follows from the identity δH≡H′−H=G′+−G+, which is multiplied by G on the left and G′ on the right to yield [6]. The formula remains valid for the pseudoinverses in the nonsingular subspace. One can calculate the change in fitness by evaluating the effect of a mutation on the Green function, G′=G+δG, and then examining the change, δF=v⊺δGf [3]. Using [6] to calculate the mutated Green function G′ is an impractical method, as it amounts to inverting at each step a large nd×nd matrix. However, the mutation of an amino acid at i has a localized effect. It may change only up to z=12 bonds among the bonds α(i) with the neighboring amino acids. Thanks to the localized nature of the mutation, the corresponding defect Hamiltonian δHi is, therefore, of a small rank, r≤z=12, equal to the number of switched bonds (the average r is about 9.3). δHi can be decomposed into a product δHi=MBM⊺. The diagonal r×r matrix B records whether a bond α(i) is switched from weak to strong (Bαα=ks−kw=+0.99) or vice versa (Bαα=−0.99), and M is a nd×r matrix with r columns that are the bond vectors mα for the switched bonds α(i). This allows one to calculate changes in the Green function more efficiently using the Woodbury formula (91, 92):δG=−GMB−1+M⊺GM+M⊺G.[12]

The two expressions for the mutation impact δG, [6 and 12], are equivalent, and one may get the scattering series of [6] by a series expansion of the pseudoinverse in [12]. The practical advantage of [12] is that the only (pseudo-)inversion that it requires is of a low-rank tensor (the term in parentheses). This accelerates our simulations by a factor of (na/r)3≃104.

Pathologies and Broken Networks.

A network broken into disjoint components exhibits floppy modes owing to the low energies of the relative motion of the components with respect to each other. The evolutionary search might end up in such nonfunctional unintended modes. The common pathologies that we observed are (i) isolated nodes at the boundary that become weakly connected via H→P mutations, (ii) “sideways” channels that terminate outside the target region (which typically includes around 8–10 sites), and (iii) channels that start and end at the target region without connecting to the binding site. All of these are some easy to understand floppy modes, which can vibrate independent of the location of the pinch and cause the response to diverge (>Fm) without producing a functional mode. We avoid such pathologies by applying the pinch force to the protein network symmetrically: pinch the binding site on face left, look at responses on face right, and vice versa. Thereby, we not only look for the transmission of the pinch from the left to right but also, from right to left. The basic algorithm is modified to accept a mutation only if it does not weaken the two-way responses and enables hinge motion of the protein. This prevents the vibrations from being localized at isolated sites or unwanted channels.

Dimension and SVD.

To examine the geometry of the fitness landscape and the genotype-to-phenotype map, we looked at the correlation among numerous solutions, typically Nsol∼106. Each solution is characterized by three vectors: (i) the gene of the functional protein, c* (a vector of length na=200 codons); (ii) the flow field (displacement), u(c*)=G(c*)f (a vector of length nd=400 of x and y velocity components); and (iii) the shear field s* (a vector of length na=200). We compute the shear as the symmetrized derivative of the displacement field using the method in ref. 21. The value of the s* field is the sum of squares of the traceless part of the strain tensor (Frobenius norm). These three types of vectors are stored along the rows of three matrices WC, WU, and WS. We calculate the eigenvectors of these matrices, Ck, Uk, and Sk, via SVD (as in ref. 36). The corresponding SVD eigenvalues are the square roots of the eigenvalues of the covariance matrix W⊺W, while the eigenvectors are the same. In typical spectra, most eigenvalues reside in a continuum bulk region that resembles the spectra of random matrices. A few larger outliers, typically around a dozen or so, carry the nonrandom correlation information.

The Protein Backbone.

A question may arise as to what extent the protein’s backbone might affect the results described so far. Proteins are polypeptides, linear heteropolymers of amino acids, linked by covalent peptide bonds, which form the protein backbone. The peptide bonds are much stronger than the noncovalent interactions among the amino acids and do not change when the protein mutates. We, therefore, augmented our model with a “backbone”: a linear path of conserved strong bonds that passes once through all amino acids. We focused on two extreme cases: a serpentine backbone either parallel to the shear band or perpendicular to it (Fig. 6).

Fig. 6.
Fig. 6.

The effect of the backbone on evolution of mechanical function. The backbone induces long-range mechanical correlations, which influence protein evolution. We examine two configurations: parallel (A and B) and perpendicular (C and D) to the channel. (A and B) Parallel. (A) The backbone directs the formation of a narrow channel along the fold (compared with Fig. 5A). (B) The first four SVD eigenvectors of the gene Ck (Top), the flow Uk (Middle), and the shear Sk (Bottom). (C and D) Perpendicular. (C) The formation of the channel is “dispersed” by the backbone. (D) The first four SVD eigenvectors of Ck (Top), Uk (Middle), and shear Sk (Bottom).

The presence of the backbone does not interfere with the emergence of a low-energy mode of the protein with a flow pattern (i.e., displacement field) that is similar to the backboneless case with two eddies moving in a hinge-like fashion. In the parallel configuration, the backbone constrains the channel formation to progress along the fold (Fig. 6A). The resulting channel is narrower than in the model without backbone (Figs. 1D and 5). In the perpendicular configuration, the evolutionary progression of the channel is much less oriented (Fig. 6C). While the flow patterns are similar, closer inspection shows noticeable differences, as can be seen in the flow eigenvectors Uk (Fig. 6 B and D). The shear eigenvectors Sk represent the derivative of the flow and therefore, highlight more distinctly these differences.

As for the correspondence between gene eigenvectors Ck and shear eigenvectors Sk, the backbone affects the shape of the channel in concert with the sequence correlations around it. Transmission of mechanical signals seems to be easier along the orientation of the fold (parallel configuration) (Fig. 6A). Transmission across the fold (perpendicular configuration) necessitates significant deformation of the backbone and leads to “dispersion” of the signal at the output (Fig. 6C). We propose that the shear band will be roughly oriented with the direction of the fold, but this requires further analysis of structural data. Overall, we conclude from our examination that the backbone adds certain features to patterns of the field and sequence correlation without changing the basic results of our model. The presence of the backbone might constrain the evolutionary search, but this has no significant effect on the fast convergence of the search and on the long-range correlations among solutions.

Acknowledgments

We thank Jacques Rougemont for calculations of shear in glucokinase (Fig. 1B) and for helpful discussions. We thank Stanislas Leibler, Michael R. Mitchell, Elisha Moses, Giovanni Zocchi, and Olivier Rivoire for helpful discussions and encouragement. We also thank Alex Petroff, Steve Granick, Le Yan, and Matthieu Wyart for valuable comments on the manuscript. J.-P.E. is supported by European Research Council Advanced Grant Bridges, and T.T. is supported by Institute for Basic Science Grant IBS-R020 and the Simons Center for Systems Biology of the Institute for Advanced Study, Princeton.

Footnotes

  • ↵1To whom correspondence may be addressed. Email: libchbr{at}rockefeller.edu or tsvitlusty{at}gmail.com.

  • Author contributions: S.D., J.-P.E., A.L., and T.T. designed and performed research and wrote the paper.

  • Reviewers: M.T., Tata Institute of Fundamental Research, National Center for Biological Sciences; and M.V., University of California, San Diego.

  • The authors declare no conflict of interest.

Back to top
Article Alerts
Email Article
Citation Tools
Request Permissions
Green function of correlated genes in a minimal mechanical model of protein evolution
Sandipan Dutta, Jean-Pierre Eckmann, Albert Libchaber, Tsvi Tlusty
Proceedings of the National Academy of Sciences May 2018, 115 (20) E4559-E4568; DOI: 10.1073/pnas.1716215115
del.icio.us logo Digg logo Reddit logo Twitter logo CiteULike logo Facebook logo Google logo Mendeley logo

Article Classifications

Proceedings of the National Academy of Sciences: 115 (20)
Table of Contents

Submit

Jump to section

You May Also be Interested in

Microplastic particles in atmospheric dust.
Microplastics in the atmosphere
Even after atmospheric microplastics settle on land or in water, they may reenter the atmosphere, a study suggests.
Image credit: Janice Brahney.
Human-modified landscape.
Anthropogenic transformation of terrestrial nature
A study examines how humans have reshaped terrestrial nature and suggests that restoring Indigenous peoples and local communities to positions of environmental stewardship may help conserve biodiversity.
Image credit: Erle C. Ellis.

Similar Articles