Proteins fulfill complex and diverse biological functions through the controlled atomic motions of their structures (Wingert et al., 2021). The protein composition is given by its amino-acid sequence, which was assumed to encode the function. However, the discovery of functional sequence variants proved that the functional encoding does not come down to the sequence, otherwise a change in the sequence would mean a change of function. Likewise, the discovery that function is fulfilled by a set of structures and not by a unique structure showed that the functional encoding does not come down to the structure either (Jaffe, 2005; Jaffe, 2020; Parisi et al., 2015).

The next alternative is that sets of atomic motions, achievable by different sequences and different structures, encode a specific function (Bahar et al., 2010). This is consistent with the multiple dynamic paths that fulfill allostery (Buchenberg et al., 2017). The challenge lies in distinguishing the set of dynamics associated with sequence variants functionally robust or functionally adapted (change of function) from the set of dynamics associated with functional failure. Inferring function from the protein dynamics is also important because pathological variants impacting the protein dynamics but not the protein structure, limit traditional structure-based drug discovery methods (Demir et al., 2021).

Network science is appropriate to study system dynamics from protein structures because it offers multiple avenues to study the complex spatiotemporal relationships between interacting entities (Barrat et al., 2004; Barabási, 2013; Unicomb et al., 2017). Integrative approaches combining experimental data or Molecular Dynamics (MD) simulation with network-based models enable to link protein structures to protein dynamics and function Demir et al. (2011), Leitner and Yamato (2018), Liang et al. (2018), Ponzoni and Bahar (2018), Bourgeat et al. (2019), Gheeraert et al. (2019), Melo et al. (2020), Bourgeat et al. (2021), Di Paola and Leitner (2021) and for a review see (Liang et al., 2020). On one hand, global mode analysis, elastic network models (ENM), dynamics network models (DNM) and protein energy networks (PEN) are used to track multiple scale dynamics in proteins, identify allosteric pathways and residues involved in biological activities. On the other, perturbation response scanning (PRS) and evolutionary network models are used to investigate the impact of mutations on protein features in particular for disease mutations. ELM applied on different protein family members also allows associating scale of motions to various types of activities (Wingert et al., 2021).

The advantages of network based models in probing protein dynamics come from the inference of amino acid and atomic links from the structure. Now, one on-going question is to clarify why network measures pinpoint functional residues (e.g., allostery) or distinguish disease mutations from the rest of the residues in order to better understand what properties amino acids have in a structure that make them functionally tolerant to mutations or not. The comparison of network measures and network models over proteins spanning large dynamics scales from enzyme to pore-forming toxins and over their sequence variants will help validating network measures as hallmarks of functional dynamics and of functional dynamic perturbations related to diseases.

One alternative perspective to network integrative approaches is to find network measures that are relevant to functional dynamics simply from protein structures. This implies a network measure probing collective slow motions and therefore shared across proteins and independent of amino acid features as observed from global modes (Bahar et al., 2010). In addition, a network measure with amino acid specific characteristics is expected if it embed the dynamics of a specific function. An allosteric enzyme and a pore-forming toxin have 3D-structures that share multiple scale collective dynamics but yet they have very different motions to fulfill their functions.

We consider the neighborhoods of each amino acid of a protein as potentially relevant to the problem. This is because on average over its neighbors every amino acid makes moderate and similar number of atomic interactions, a property shared by many different proteins (Dorantes-Gilardi et al., 2018). In addition, each neighborhood is different in terms of number of neighbors and type of neighbors (Dorantes-Gilardi et al., 2018). Thus, neighborhoods satisfy the two conditions to embed protein dynamics. Moreover, the neighborhoods describe the spatial position of the amino acids in the structure, which carves the space occupied by the amino acid atoms and hence uncovers the space left available between amino acids where atomic motions can take place. The relation between the space occupied by entities and the system dynamics is a broad topic from granular material to urban and protein systems (Liang and Dill, 2001; Majmudar and Behringer, 2005; Henzler-Wildman and Kern, 2007; Majmudar et al., 2007; Barthelemy, 2011; Dorantes-Gilardi et al., 2018; Gheeraert et al., 2019; Naganathan, 2019).

To analyze neighborhoods, a protein structure is modeled by an amino acid network (AAN), where the nodes are amino acids and the links are atomic interactions between amino acids, inferred from atomic proximity (Dorantes-Gilardi et al., 2018). The space occupied by neighborhoods is described in terms of amino acids with the node degree and in terms of atoms with the node weight (Methods). Classically, a unique cutoff around the threshold for chemical interactions (5 Å) is used to investigate protein features but here the neighborhoods are computed at different cutoff distances to probe the space occupancy at different spatial scales (Vuillon and Lesieur, 2015; Viloria et al., 2017). This condition is necessary to track the multiple dynamics scales associated with functional dynamics (Henzler-Wildman and Kern, 2007; Munoz and Cerminara, 2016).

Our case study is the third PDZ domain of the synaptic protein PDS-95 (PDB 1BE9) chosen because the functional impact of most of its single amino-acid mutations is known from experiments as well as some double mutations (McLaughlin et al., 2012; Salinas and Ranganathan, 2018). Thus, this case study is appropriate for future validation of the network measures to link dynamics and function.

Starting from the Protein Data Bank (PDB) data, protein structures are modeled using the Amino Acid Network (AAN), an established model in Computational Biology (Dorantes-Gilardi et al., 2018). The AAN is a graph G = (V; E), with V is the set of the N nodes of the network (vertices of the graph) and E the set of links of the network (edges of the graph).

Nodes of the AAN: Each node in the AAN corresponds to one amino acid of the protein’s structure named according to the protein sequence: V = { i | i   is   an   amino   acid } . (1) Links of the AAN: A link is an atomic interaction defined by atomic proximity: two amino acids i and j are connected if there exists at least one couple of atoms, one belonging to i and one belonging to j, at a distance lower or equal to a given threshold (Cutoff distance c), E = { ( i , j ) | i , j ∈ V   with   i ≠ j   and   ∃ ( atom i ∈ i , atom j ∈ j ) with   dist ( atom i , atom j ) ≤ c } . (2)

Link weights of the AAN: Each link is weighted according to the number of atomic couples that respect the cutoff condition: w i j = | { ( atom i ∈ i , atom j ∈ j ) with   dist ( atom i , atom j ) ≤ c   and   i ≠ j } | (3) where the pipe symbol | denotes the cardinality of the set (i.e. the number of elements of the set). When c = 5 Å, that is a threshold for chemical interactions, the link weights measure the number of atomic interactions between two amino acids.

Packing around amino acids: In the AAN, the node degree k _i , defined as the number of amino acid neighbors of a node i, measures the amino-acid packing around the amino acid i, referred to as the amino-acid neighborhood. The node weight w _i is defined as the sum over all the weights of the links that connect the node i to its neighbors ( w i = ∑ j ∈ N ( i ) w i j with N(i) the set of neighbors of node i). The node weight measures the atomic packing around the amino acid i, referred to as the atomic neighborhood.

Cutoff distance: Different cutoff distances are used in this study such that the packing around each amino acid is described at different length-scale via the neighborhoods at variable cutoffs. The cutoffs are integers and range from 3 to 11 Å such that the packing within chemical reach (≤5 Å) and above chemical reach (>5 Å) are monitored. The rational is to distinguish amino acids by their ‘chemical’ neighborhoods and above chemical reach neighborhoods to probe space occupancy involved in multiple spatiotemporal scales.

Plateau versus linear degree dependencies: The degree dependency with the cutoff is plotted. Some amino acids show a linear dependency of the degree with the cutoff and are referred to as linear amino acids. Some amino acids exhibit a plateau over some cutoffs determined with the derivative (here just Δk because the cutoffs are consecutive integers) equals to zero or equals to one if and only if the derivative at cutoffs before the plateau is equal to four or more. We refer to these amino acids as plateau amino acids.

Amino acid side chain length classification: Amino acids are classified by side-chain lengths as follows. Side chain length <3 Å are small amino acids (G, A, P, S, V, I, T, C), side chain length between 3 Å ≤ length <5 Å are medium amino acids (L, E, D, H, N, Q, M) and side chain length ≥5 Å are big amino acids (F, K, R, Y, W).

1D-barcode: The 1D-barcode represents the degree dependency with the cutoff (plateau- or linear-) of the amino acids of the 1BE9 sequence (1D).

2D-barcode: The 2D-barcode represents the degree dependency with the cutoff of the amino acids of the 1BE9 sequence (1D) at variable cutoffs (2D).

The AAN of the PDZ domain of the synaptic protein PDS-95 is generated using the PDB 1BE9 and the weights and degrees of every amino acid nodes are computed at different cutoff distances (Methods). We recall that the weight of a node describes its atomic packing, referred to as its atomic neighborhood while the degree of a node describes its amino-acid packing referred to as its amino-acid neighborhood. The atomic neighborhood takes into account the features of the amino acids. The weight and degree cutoff dependencies are plotted to investigate how atoms and amino acid neighbors occupy the space around each amino acid at different scales in the protein structure. This is a proxy of the dynamics as the more space occupied the less space left available for atomic motions.

The weight cutoff dependencies are quadratic indicating that the weight (the number of atomic interactions) increases with the square of the cutoff distance, i.e., with a surface of contact between atoms rather than a volume of contacts (Figure 1A; Supplementary Figure S1A). The degree cutoff dependencies are not quadratic but exhibit two distinct behaviors: neighborhoods with linear dependency with the cutoff (Figure 1B, top in the protein structure) and neighborhoods with a plateau over some cutoffs (Figure 1B, bottom and Supplementary Figure S1B). We call linear (plateau) amino acids, the amino acids whose degree follows a linear (plateau) dependency with the cutoff. The plateau is due to a lack of amino acid neighbors at higher cutoffs and not to a lack of atoms since the weight cutoff dependencies have no plateau (Supplementary Figure S1A). The plateau can result from having big amino acid neighbors that occupy the space over more than one cutoff versus having small amino acids that lead to a linear increase over the same cutoffs (Figure 2A).

Ten nodes have a plateau between 4 and 5 Å, twenty between 5 and 6 Å, one between 6 and 7 Å and six over 4 to 6 Å (Figure 1B).

All amino-acid types observed in 1BE9 except histidine which appears only twice, adopt plateau and linear neighborhoods (Figure 1C). Methionine, cysteine and tryptophan are not present in 1BE9. Twenty-six surface-exposed amino acids out of fifty-eight are plateau (∼half) and sixteen buried amino acids out of fifty-seven are plateau (∼a third). 27% of the β-strand amino acids are plateau against 35 and 48% for loops and α-helix amino acids, respectively. Thus the plateau and linear neighborhoods are achieved regardless amino-acid type, position in the structure and secondary structure, which makes the degree cutoff dependencies insensitive to amino acid features as global modes. Together with the plateau built at chemical-interaction threshold, it is consistent with the degree cutoff dependencies probing collective motions. Accordingly, a 1D-barcode representing the plateau- or linear-neighborhood of each amino acid of the sequence could be used to characterize the multiple scale dynamic features of a protein (Figure 1C).

To have a plateau or a linear degree cutoff dependency for one amino acid type implies customizing neighborhoods as illustrated by the higher number of big amino acid neighbors for plateau amino acids than linear amino acids (Figure 2A, compare top and bottom). To have it for any amino acid type also implies neighborhood customized to the central amino acids as seen on Figure 2A (compare across amino acid types). Thus, the linear and plateau neighborhoods accessible to all amino acids are nevertheless built from specific amino-acid neighborhoods (Figure 2A; Supplementary Figure S2). This means the protein dynamic-functional specificity could be embedded in the collective motions through specific space occupancy arising from neighborhood diversities. A 2D-barcode representing the amino acid neighborhood size specificities of each amino acid of the protein sequence (Figure 2B, horizontal axis) at different cutoffs (Figure 2B, vertical axis) could be used to characterize the protein specific embedded dynamics (Figure 2B).

We can see from this single case 2D barcode classes of neighborhoods in terms of neighbor sizes which anticipate classes of dynamics and of responses upon mutations supporting the possibility to use the data for mining dynamics and its relation to function (Figure 2B). Some positions are composed of a majority of one-size neighbors over the cutoffs (e.g. small neighbors: R318, T321, G322, L349, N363; medium: F301, G303, R312, G333, G351, A390) while others are a mixture of neighbor sizes (e.g. A308, I341, K355, A383, Q391, E401). In addition, some positions change neighbor sizes over the cutoffs (e.g., E396, E401). However, these sole data do not have the statistics to make hypothesis on which neighborhood classes lead to which dynamic classes or draw conclusion between the neighborhood specificity, dynamics and function.

The study shows that amino acid neighborhoods and not only amino acids and amino acid pairs, contain information relevant to protein dynamics, opening a new perspective to explore the link between dynamics and function. The 1D barcode of, for example, an enzyme and a pore-forming toxin can be compared to determine common and distinct features which can in turn be analyzed with the 2D barcode to survey both the neighborhood diversity of the common 1D barcode features and of the distinct 1D barcode features assuming the former identifies positions functionally insensitive (protein aspecific) and the later positions functionally sensitive (protein specific). The analysis of each protein variants can be used to validate the assumption as well as database analysis.

This will contribute to diagnosing dynamic functional faults and dynamic functional diversity.