SenseNet, a tool for analysis of protein structure networks obtained from molecular dynamics simulations

Computational methods play a key role for investigating allosteric mechanisms in proteins, with the potential of generating valuable insights for innovative drug design. Here we present the SenseNet (“Structure ENSEmble NETworks”) framework for analysis of protein structure networks, which differs from established network models by focusing on interaction timelines obtained by molecular dynamics simulations. This approach is evaluated by predicting allosteric residues reported by NMR experiments in the PDZ2 domain of hPTP1e, a reference system for which previous computational predictions have shown considerable variance. We applied two models based on the mutual information between interaction timelines to estimate the conformational influence of each residue on its local environment. In terms of accuracy our prediction model is comparable to the top performing model published for this system, but by contrast benefits from its independence from NMR structures. Our results are complementary to experimental data and the consensus of previous predictions, demonstrating the potential of our new analysis tool SenseNet. Biochemical interpretation of our model suggests that allosteric residues in the PDZ2 domain form two distinct clusters of contiguous sidechain surfaces. SenseNet is provided as a plugin for the network analysis software Cytoscape, allowing for ease of future application and contributing to a system of compatible tools bridging the fields of system and structural biology.


Introduction
Protein structure networks map atoms from a protein structure to nodes and define edges to represent atom interactions, e.g. contacts and hydrogen bonds. The resulting networks may be used to predict e.g. allosteric communication pathways [1][2][3] with potential applications in innovative drug design [4][5][6][7][8]. Most commonly, such analyses are based on individual crystal structures and rely on centrality measures such as betweenness centrality (BC) or characteristic path length centrality (CPLC) to identify functionally important residues [1][2][3]9]. However, application of these algorithms to experimental structures of e.g. the PDZ domain did not provide results consistent with experiment [10]. It has been generally recognized that highly dynamic effects such as allostery, which are not always associated with stable conformations, are difficult to study solely on the basis of individual experimentally obtained structures [5,[11][12][13]. Computational methods for analyzing structure ensembles obtained from e.g. molecular dynamics simulations (MD), which capture the dynamic behavior of proteins, are therefore attractive for allosteric prediction [11,[14][15][16][17][18][19][20]. Several tools exist for analysis of structure ensemble networks, among them xPyder [21], PyInteraph [22], MD-TASK [23], gRINN [24], PSN-Ensemble [25], NAPS [26,27], RIP-MD [28], Bio3D [29], MDN [30] and the Cytoscape plugin RINalyzer [31]. A common approach for network analysis of MD data is to define edges by correlation analysis of atomistic motions, which comes at the cost of losing structural and conformational details of the underlying interactions. In addition, many approaches use a rigid mapping of one node per residue, preventing the combination of different levels of resolution, e.g. to separate information flow between backbone and sidechain atoms. Finally, the majority of tools are provided as standalone programs or webservers, making it difficult to combine different algorithms within a single analysis session. To address these limitations, we developed SenseNet, a plugin for the free network analysis software Cytoscape [32]. SenseNet is based on an alternative strategy to scalar correlation coefficients, namely associating edges with MD-based timelines, which allow to track the evolution of interactions during a simulation by checking their existence at predefined timeslots. This representation allows for a larger variety of analyses than correlation-based approaches, like e.g. interaction averages, lifetime analysis, frame clustering, or shared information between timelines.
Ligand binding often modulates protein function by triggering conformational changes distant from the binding site. A major goal of computational allosteric prediction is to identify key residues sensing ligand binding events over long intramolecular distances; in the context of computational predictions, these residues are commonly labeled as "allosteric". For the purpose of evaluating these methods, PDZ domains are a well-established reference system. Members of this abundant domain class commonly bind C-terminal or short internal peptide sequences and participate in allosteric interactions with other domains [33,34], serving as initiators and mediators of protein assembly processes [35][36][37]. Although the domain is allosterically modulated by its peptide ligands, crystal and solution NMR structures of the PDZ2 domain of hPTP1e (human Protein-Tyrosine Phosphatase 1e) show no substantial conformational changes between apo and ligand bound states [38]. Therefore, the relationship between structure, dynamics, and allostery in the PDZ2 domain of hPTP1e was explored by Lee and coworkers, who identified a number of allosteric residues by probing the effects of ligand binding and point mutations on NMR backbone and methyl side chain dynamics [38][39][40]. However, open questions remain concerning the contribution of residues lacking methyl groups and how individual residues act together to form allosteric pathways, motivating structurebased computational prediction as a complementary strategy [41]. Methods previously applied to the PDZ2 system include interaction energy and correlation networks [42,43], elastic network models [44], hydrogen bond heat diffusion pathways [45], relative entropy networks of distance distributions (REDAN) [46], and coordinate fluctuations [47,48]. Furthermore, specialized simulation techniques were employed such as perturbation response scanning [49], rigid residue scan (RRS) [50], and NMR guided simulations [10,51]. However, results reported by computational studies have shown considerable variance, warranting efforts to consolidate and improve prediction models [41].
In this work, we present our network analysis software SenseNet and evaluate two of therein implemented, timeline-focused algorithms to find pathways of allosteric information transfer in the PDZ2 domain. By quantifying how much information the timelines of physical interactions provide about their environment, we obtained accurate models for predicting allosteric residues in PDZ2. Finally, we propose a consolidated allosteric model combining our results with experimental data and the consensus of previous predictions, which suggests that PDZ2 contains two allosteric pathways formed by clusters of contiguous sidechain surfaces.

Algorithms
Protein structure networks based on interaction timelines. In a structure network as implemented in SenseNet, each node (which together form the set of nodes N) represents a single atom or a group of atoms while edges represent interactions between nodes. If several interaction types (e.g. contacts or hydrogen bonds) are present, a node pair may be connected by more than one edge. Every interaction is associated with a timeline, representing the different states of the interaction in the analyzed ensemble of structures, e.g. simulation frames from an MD trajectory. We define an atomistic timeline as the vector " where α, β are nodes representing single atoms, k is an interaction type and t is a simulation time frame (bold type face denotes matrices and vectors). Timelines of edges connecting two atom groups (e. g. residues) are calculated as in which i, j are nodes representing atom groups. The connectivity between nodes is given by the symmetric adjacency matrix for each interaction type k. In combination, the sets of nodes and edges form a network which encodes both the structural topology of the protein system and the fluctuations between different conformational states through its interaction timelines. Those features can then be subjected to further analyses in order to gain insights into the dynamic behavior of the protein system. Note that in cases where the network is based on a single structure instead of an ensemble of structures, the network model reduces to a simple form where each timeline has a length of one and corresponds to the number of interactions between the connected nodes. Allosteric prediction based on correlation between interaction timelines. We propose two novel algorithms, the node correlation factor (NCF) and difference node correlation factor (DNCF), to predict residues associated with allosteric function in proteins. Our model presupposes that in order for a residue to have an observable allosteric function, its conformations must be correlated to conformational changes in its immediate environment. The conformational states of all residues are encoded within the interaction timelines in the network. We define the immediate environment as the interactions represented by neighboring edges, i.e. edges which are separated by at most a single node. Hence, we begin by considering how each interaction is correlated to interactions in its immediate environment. By applying this definition, we obtain the edge neighbor correlation factor (ECF) as with i, j belonging to the node set N, k and l being part of the interaction type set K, and I is the mutual information function in which p(x, y) represents the joint probability of values x and y and p(x) corresponds to the marginal probability of state x in timeline X. The mutual information function is a non-linear measure of correlation quantifying the information shared between timelines, i.e. the increase of predictability of the states in timeline X if the other timeline Y is observed [52]. Furthermore, χ represents an indicator function selecting the neighboring edges of i, j, k and is defined as where d is the Kronecker delta and the δ kl term serves to exclude the self-information of edge i, j, k. The definition ECF score is intuitively illustrated using the network shown in Fig 1. The ECF score of the blue edge is calculated as the sum of mutual information contributions between the blue edge and all its neighboring edges, shown in green. Each contributing mutual information term indicates the strength of correlation between the interaction represented by the blue edge and the respective neighboring interaction. If the interaction states represented in the timeline of the blue edge are strongly correlated to the interaction states of its surrounding edges, it will lead to a high ECF score, suggesting that changes in one interaction may affect its immediate environment; In other words, information about conformational states could then potentially be transmitted via these strongly coupled interactions. Summing up the ECF scores of a node's adjacent edges gives the node correlation factor (NCF) which can be expressed as and highlights residues with strong conformational coupling. These residues, as they participate in interactions that may transfer information to their environment, are thus likely candidates for showing behavior associated with protein allostery.
As an extension to the model, another aspect can be considered for the prediction of allosteric residues, namely the conformational differences between two states of a protein system, e.g. ligand bound and ligand free. The difference node correlation factor (DNCF) quantifies changes in timeline coupling between two networks, each created from a different MD trajectory simulating either the ligand bound or the ligand free state. After selecting one trajectory as the reference and the other as the target, the definition of Eq 5 is adjusted to with X; Ŷ denoting the timelines from the reference simulation matching the locations of X and Y of the target simulation andp representing the probabilities of the reference timelines. Note that edges which exist solely in the reference network do not contribute, therefore the score is not symmetric with respect to interchanging target and reference networks. Substitution of Eq 8 in Eq 4 yields the DNCF score. The DNCF score measures the change in shared information between equivalent interaction timelines in the target and reference systems. This can be illustrated with the following example: Suppose there are two neighboring interactions obtained from MD simulations of the system, and the timelines show that they are strongly correlated. Then the same system is simulated again, but now including a ligand bound to an allosteric site, which are sensed by residues associated with allosteric function. The binding of a ligand to an allosteric binding pocket is likely to change the nature and efficacy of information transfer within the protein, which can manifest stronger or weaker coupling between interaction timelines. The DNCF score is composed of the pointwise mutual information contributions of the allosterically activated system as encoded in timelines X and Y, from which the contributions of the equivalent reference timelines X and Ŷ are subtracted. Thus, high DNCF scores are expected from residues for which the coupling of interactions changes between the target and reference network, i.e. before and after binding of a ligand to an allosteric site. An essential feature of our model emerges from the definitions of the ECF, NCF and DNCF scores, namely the explicit locality of network effects. By limiting our analysis on the shared information between adjacent residues in the network, the influence of spurious correlation is reduced. To illustrate, consider that any pair of residues in a protein, no matter how far apart, would be compared. This would lead to a drastic increase of evaluated correlation terms, and thus more residue pairs showing high correlation by pure chance. At the same time, the probability that two residues influence each other directly in a substantial manner (i.e. without detectable changes in the residues between them) is lower if they are far apart, especially as the physical interactions included in our analysis, i.e. hydrogen bonds and carbon contacts, are of limited range. Adding up contributions of distant residues would thus substantially increase the noise introduced in the analysis. Instead, we propose that in most cases it is more productive to focus on the identification of neighboring residues directly exchanging information, and to analyze how they build chains of signaling residues. However, in instances of allosteric communication lacking this locality of effects, other methods may be more accurate.
Network node centrality methods for allosteric prediction. Measures of node centrality are commonly used to detect functional residues using protein structure networks [1,2,9]. When applying these methods to prediction of allosteric residues, it is postulated that residues important to transferring signals between functional sites are related to the most central nodes in the structure network, i.e. nodes that are essential when walking the shortest path between nodes along network edges. SenseNet implements two centrality functions for this purpose: Betweenness centrality (BC) finds those nodes which are located on the largest number of shortest paths over all possible node pairs [1,53]. It is defined as where i, j, k belong to the set of nodes N, σ jk is the number of shortest paths between j and k, and σ jk|i is the number of shortest paths between j and k passing through i. The second method implemented in SenseNet is characteristic path length centrality (CPLC) [9]. For this method, nodes that are crucial for maintaining the shortest paths are presumed to be key to communication, as measured by the robustness of shortest paths to the removal of individual nodes [9]. In order to determine the robustness of the network, the characteristic path length, i.e. the average length of shortest paths in the network is considered as where N is the set of nodes, N p is the number of node pairs in the network and d(i, j) is the minimum number of edges to be traversed between i and j. The CPLC score corresponds to the effect of removing a node on the characteristic path length of the network, which can be expressed as where L i is the characteristic path length of the network after removal of node i. The BC and CPLC algorithms are commonly applied to individual (crystal or NMR) structures and do not trivially transfer to structure ensembles from MD simulations. This is because the networks obtained from MD simulations contain a large number of additional spurious interactions in the network compared to a crystal structure. Since Eqs 9 and 10 utilize the shortest paths between nodes along a chain of edges without accounting for the stability of the interaction, an interaction present only in a tiny fraction of the simulation could be considered with the same importance as more long-lived, substantial interactions. In contrast, NCF and DNCF methods intrinsically limit the influence of spurious interactions due to the explicit locality of contributing interactions and by definition through the mutual information function. For this reason and the fact that BC and CPLC are most commonly used with individual structures, we applied these methods only to networks obtained from crystal and NMR structures.

Molecular dynamics simulations
MD simulations in this work are based on the crystal structures of hPTP1E-PDZ2 in the apo state (PDB-ID: 3LNX) and bound to the C-terminal peptide of RA-GEF-2 (PDB-ID: 3LNY) as well as the corresponding solution NMR structures 3PDZ and 1D5G, using the first model provided in the files. These NMR structures were chosen to allow for direct comparison with previous studies [10,39]. Protein and ligand residues missing in the crystal structures were added based on their NMR structure analogues using Modeller 9.18 [54], creating 100 candidate structures and selecting the model with the best DOPE score for simulations and network analyses. MD simulations were performed using the Amber16-AmberTools17 software suite [55] with the Amber14SB force field [56] and TIP3P water [57]. The system was solvated in a cubic water box using a minimum solute-face distance of 12 Å and 150 mM NaCl. For the nonbonded interactions a 12 Å direct space cutoff and PME summation for electrostatic interactions were applied. Energy minimization was performed until convergence to 0.01 kcal � mol -1 � Å -1 was reached using the XMIN minimizer. Afterwards, the volume of the solvent box was adjusted to a solvent density of 1.00 kg � m 3 . For all simulations a time step of 1 fs was applied and SHAKE [58] was used for hydrogen-containing bonds. Systems were gradually heated from 0 to 300 K over 1.7 ns using a variant of our published heatup protocol [59], restraining all heavy atoms by 2.39 kcal � mol -1 � Å -2 until 20 K and all backbone atoms until 200 K. For the first 1.2 ns of the heatup a Langevin thermostat was used with a collision frequency of 4 ps -1 and for the last 0.5 ns a Berendsen barostat was employed with a relaxation time of 2 ps. Afterwards the NPT ensemble was used with a slow coupling Berendsen thermostat at 300 K (coupling time: 10 ps) in combination with a Berendsen barostat (relaxation time: 5 ps). For each system, ten independent simulations were performed for 1 μs each (based on separate heatup runs and different randomized Langevin seeds). The initial 100 ns of each replicon were removed before analysis to reduce bias towards initial structures. Trajectory postprocessing was performed with CPPTRAJ [60], using the "nativecontacts" command for contact timelines of carbon atoms (saving both native and nonnative time series), and the "hbond" command for hydrogen bonds (distance cutoff 3.5 Å; angle cutoff 135˚). The data generated by CPPTRAJ provided the interaction timelines for all network analyses based on MD trajectories, i.e. for the NCF and DNCF methods. Interaction data for BC and CPLC analyses were extracted directly from the corresponding PDB files using AIFgen with equivalent settings for interactions and distance/angle cutoffs as detailed for CPPTRAJ (see example script in S2 File).

Protein structure networks
For analyses of protein structure networks and related quantities we used the SenseNet plugin (version 1.0.0) for Cytoscape (version 3.6.1) [32]. In order to create a network, SenseNet requires a list of atom-atom interaction timelines, where each interaction is defined by a minimum of one source atom, one target atom, an interaction type (e.g. hydrogen bond), and a timeline represented as a list of interaction values corresponding to each time frame (e.g. a list where 1 indicates presence of an interaction, while 0 indicates absence in each given frame). As a general input data format for SenseNet, we defined the AIF file format, which provides a list of interaction timelines as a structured text file that can be easily created, inspected and modified using a text editor (see S2 File for an example of the format). SenseNet provides tools for automatic generation of AIF files from multiple sources. Lists of interaction timelines as created by the CPPTRAJ "hbond" and "nativecontacts" analyses can be directly converted into AIF format using the SenseNet GUI or AIFgen, which provides a command line interface to the GUI functions available in SenseNet. Alternatively, SenseNet and AIFgen can extract timelines of pairwise contacts or hydrogen bonds directly from PDB files using the same criteria as implemented in CPPTRAJ. Example scripts demonstrating the workflow for AIFgen for converting CPPTRAJ outputs and extraction of interactions from PDB files are given in S2 File. For this work, we converted CPPTRAJ outputs of contact and hydrogen bond analyses into AIF files using AIFgen (version 1.0.4).
ECF scores were calculated with SenseNet using the therein implemented "Correlation" function set to the "Mutual information" mode. Then, the "Degree" function was used to sum over the ECF scores calculated in the previous step. DNCF scores were calculated after importing first the reference and target systems (see Eq 8) as separate networks. As references in the context of DNCF calculations, we selected the network generated from the corresponding ligand bound simulation for the analysis of the network of the free protein, and vice versa. The DNCF scores were calculated using the "Correlation" function set to "Mutual information difference". The obtained edge scores were then summed up using the "Degree" function. Edges of the two networks were considered equivalent if they connected the same residues and were of the same interaction type (Edge mapping in SenseNet set to "Match Location"). Contact betweenness centralities (BC) [53] and characteristic path length centralities (CPLC) [9] were calculated using the respective modes within the "Centrality" function and normalized using the min-max procedure. For high throughput analyses, we used the CyREST interface of Cytoscape to call the corresponding SenseNet functions. Plots were generated using matplotlib (

Prediction of allosteric residues
Predictions were verified against methyl sidechain dynamics data [39], using classifications as allosterically active and inactive as defined by Cilia et al. ("NMR dataset", n = 25, see S1 Table) [10]. In that study, backbones of NMR structures and Monte Carlo sampling were used to find correlated side chain torsions. As this method was not applicable to alanine residues, the authors evaluated prediction performance using either the complete NMR dataset or a variant excluding alanine residues ("NMR-Ala dataset", n = 21). To be consistent with these former studies, we chose to adopt this scheme in this work. Receiver Operating Characteristic (ROC) curves were generated by plotting, for various prediction score thresholds, the corresponding False Positive Rates (FPR) and True Positive Rates (TPR) with False Positives (FP), True Positives (TP), False Negatives (FN) and True Negatives (TN) according to the NMR datasets. In addition, we generated Precision-Recall (PR) curves based on Precision (PPV) and Recall (equivalent to TPR) scores. The overall prediction performance was evaluated by calculating the area under the curve for both ROC (rocAUC) and PR plots (prAUC) using trapezoidal integration.

Features and Implementation of SenseNet
SenseNet reads interaction data from structure ensemble files in PDB format or MD trajectory analysis outputs generated by CPPTRAJ [60]. By default, each node corresponds to a single amino acid and edges represent interactions on the amino acid level. SenseNet automatically determines the network topology from these timelines (Fig 2A), offering different adjustment options from removing rare interactions to considering only certain interaction types. Different levels of timeline analyses are possible, as users can either scroll through single time frames to investigate e.g. network evolution or time-dependent interactions, or analyze time-averaged networks. At any point during a running session, residue level nodes and associated interactions can be split into individual atoms, allowing for system specific tailoring of different resolution levels. As an example application providing a detailed demonstration of this concept, we refer to our previous study analyzing the recognition of different DNA modifications by the protein UHRF1 [64]. SenseNet's user interface is separated into the main network and three control areas (Fig 2B). The left panel allows access to implemented analysis functions and displays visualization status information, such as the selected edge weighting scheme or a bar to scroll through different time frames of the network. Whenever an analysis is performed, a summary of obtained results appears on the right panel, either as tables or plots. In addition, results are written into the node and edge data tables in the bottom region, from where they can be utilized by other analysis functions, either by SenseNet or other tools. This workflow, in

PLOS ONE
combination with side-by-side network and structure visualization, allows for a rapid explorative cycle of performing quantitative analyses and intuitive exploration of the underlying structural details.
For quantitative analysis of timeline data, SenseNet offers functions for calculating timeline correlation, entropy, autocorrelation, lifetime, clustering, and network comparison. In addition, search algorithms for shortest paths as well as centrality measures are provided. Analysis results are presented as tables or plots and can be exported as raw data or images. For large scale workflows, analyses can be automated via batch script files or the CyREST interface. Network and structure visualization can be carried out in parallel by connecting SenseNet to the PyMOL [63], VMD [62], or UCSF Chimera [65] structure viewers, automatically highlighting selected nodes and edges from the network in the protein structure.

Evaluation of allosteric prediction methods using the PDZ2 domain
First, we reinvestigated the allosteric prediction performance of betweenness centralities (BC) and characteristic path length centralities (CPLC) based on networks generated from NMR and crystal structures, which had previously shown poor prediction performance for the PDZ2 system with CPLC as the best performing centrality model [10]. This allowed us to verify our implementation and to compare different network methods based on the same dataset. In line with the aforementioned work, we determined ROC and PR curves measuring the prediction accuracy of tested models with respect to the NMR dataset, which is composed of allosteric and non-allosteric residues based on methyl sidechain dynamics, and the corresponding NMR-Ala dataset variant excluding alanines [10,39] (S1 Table). In an attempt to replicate the network centrality predictions from Cilia et al. (NMR: 0.54, NMR-Ala: 0.59) [10], we calculated CPLC scores based on the crystal and NMR structures of the PDZ2-RA-GEF-2 complex using a carbon contact distance cutoff of 5 Å. For the NMR structure, resulting rocAUC scores were very close to the previously reported values (NMR: 0.55, NMR-Ala: 0.56) and only modestly higher for the crystal structure (NMR: 0.65, NMR-Ala: 0.69), indicating that the differences are only due to subtly differing details in network implementations.
In contrast to the centrality approach, interaction timelines generated from structure ensembles allow to additionally analyze the correlation between interactions, as quantified by the NCF and DNCF scores (see Materials & Methods). In general, residues with high NCF scores provide information, through linear and nonlinear correlation, about the interaction state of their environment. While the NCF estimates the information of residues within a single simulation, the DNCF score models the corresponding differences between two simulations, e.g. with and without a ligand. In order to obtain the structure ensembles necessary for calculation of these scores, we performed ten 1 μs MD simulations of the free PDZ2 domain and the PDZ2-RA-GEF-2 peptide complex. Timelines of contacts and hydrogen bonds were extracted and converted into protein structure networks using AIFgen and analyzed using SenseNet. First, we systematically evaluated all compared network methods (BC, CPLC, NCF, DNCF) using a grid search of 48 parameter combinations (S2 Table). These combinations were obtained by varying the contact distance cutoff from 4 to 9 Å, the interaction subset settings (all or only inter-sidechain interactions), and networks generated from different sources (apo-or peptide-bound structures; NMR or crystal structures). To understand which parameters are most important for prediction performance, we grouped all data points according to these categories followed by analysis of the obtained rocAUC score distributions. In the following, we focus predominantly on the results obtained for the NMR-Ala dataset, as alanine residues proved to be particularly difficult to predict for all methods tested here as well as those previously published. Fig 3A shows that average rocAUC scores over all combinations were consistently highest for the DNCF method, followed by NCF and finally CPLC and BC, which registered 8-11% lower average AUC scores compared to the former methods. In a more detailed view (Fig 3B), we observed that on average, prediction performances improved if apo PDZ2 was used as starting structure compared to peptide bound systems, with relatively small differences for CPLC, BC, and DNCF (up to 5%), but more substantial improvements for NCF

PLOS ONE
(up to 9%). Interestingly, the NCF prediction performance based on the apo systems was almost as high as the DNCF scores although, in contrast to DNCF, they do not contain any information about the ligand. Regarding the set of included interactions in the network ( Fig  3C), rocAUC scores increased on average by 2-4% if only inter-sidechain interactions were considered. Finally, analysis of contact cutoff distances shows that BC and CPLC method performances appear to peak at 6 Å, whereas a 4 to 5 Å cutoff worked best for the DNCF and NCF methods (Fig 3D). Observing the shape of rocAUC distributions and the lower performance limit for worst-case parameters can give an indication about the sensitivity of a method to choosing inappropriate network parameters. For BC and CPLC methods, several parameter combinations led to essentially random prediction performance (rocAUC~0.5) (Fig 3), indicating a high sensitivity to parameter choices in order to achieve good accuracy. In contrast, NCF and even more so DNCF were consistently more robust, as they showed better performances even for suboptimal parameters over all categories (Fig 3). Many of the observed trends are reflected, to a lesser degree, on the full NMR reference set which includes alanine residues (S1 Fig). In conclusion, we first observe that all parameter categories follow consistent trends, highlighting the importance of parameter choice for prediction quality, which is particularly true for methods based on centrality. Second, this consistency is also observed if the different methods are compared, i.e. the favorable performances of NCF and DNCF models relative to centralities are reflected throughout all parameter settings.
The best performing CPLC model was obtained for the apo PDZ2 crystal structure and a carbon contact cutoff of 6 Å in a sidechain exclusive network, interestingly differing from the original evaluation discussed above (5 Å and including backbone interactions) [10]. Using the optimized parameters, the rocAUC score for the NMR-Ala dataset increased by 5% to 0.74, while performance for the NMR dataset degraded by 1% to 0.64, respectively ( Table 1). The corresponding prAUC scores increased by 2% for the NMR dataset (0.75 to 0.77) and 5% for NMR-Ala (0.78 to 0.83). The BC method performed optimally with the same parameter set as CPLC, but with about 3 to 4% lower rocAUC scores (Table 1). Overall, only modest performance improvements could be achieved for the BC and CPLC methods by variation of network parameters.
For both DNCF and NCF models, the optimal parameter set consisted of a 4 Å contact cutoff in a sidechain exclusive network using simulations of the apo-NMR PDZ2 structure. Of all settings tested in the parameter search, DNCF was found to be the best overall predictor, achieving a rocAUC of 0.71 and prAUC of 0.82 on the full NMR set, which corresponds to a 5 to 7% improvement compared to the CPLC model. Accordingly, the performance on the NMR-Ala set was also higher than for the centrality methods with a rocAUC of 0.81 and a prAUC of 0.88. The best NCF model showed similar overall trends, but individual AUC scores were 1-5% lower (Table 1). In line with most published methods, rocAUC scores were

PLOS ONE
consistently 7-10% lower for the NMR dataset compared to NMR-Ala, which highlights the general difficulty for predicting this residue type (Table 2). In order to obtain sufficient statistical sampling for the determination of optimal model parameters, we performed a total of 10 μs of simulations, which constitutes an increasingly common but still substantial computational effort at this time for a system the size of PDZ2. While such an effort is justified for evaluation studies, for practical and effective application a guideline as to what amounts to a reasonable simulation time should be established. To gain a rough estimate of this and the convergence of our model, we repeated our analysis using the DNCF model with optimal parameters, but with truncated trajectories for each replica. The first analysis was performed on trajectories shortened to contain only the first 100 ns (after removing the initial 100 ns to reduce replica bias towards the initial structure, as detailed above), yielding a cumulative simulation time of 1 μs (10 x 100 ns). Then, subsequent analyses were performed on the first 200 ns yielding a cumulative time of 2 μs, then 300 ns for 3 μs, and so on. This approach was chosen since it shows directly how our results would have changed had we chosen a shorter simulation time for our analysis. The obtained DNCF scores were compared to the NMR-Ala and NMR datasets and rocAUC and prAUC calculated accordingly (Fig 4A and 4C). These data indicate an improvement of prediction performance up until about 3 μs of cumulative simulation time, and remaining approximately constant past that point. Taking those 3 μs as the target time, we proceeded to determine whether it was more beneficial to use fewer replicas with longer individual simulations, or to use more replicas in combination with shorter simulation times. Thus, we compared predictions using between four and ten replicas, taking the appropriate amount of simulation frames from each replica to reach a total simulation time of 3 μs. For example, when using four replicas, each replica trajectory contributed 0.75 μs (total 3 μs from 4 x 0.75 μs), whereas for five replicas each contributed 0.6 μs, and so on. This analysis was performed for each possible combination of replicas, e.g. for four replicas we considered all ways to pick four replicas out of the total of ten replicas. Judging from both the means and standard deviations of rocAUC/prAUC results (Fig 4B and  4D), it is clearly beneficial to use up to 8 replicas, corresponding to 8 replica simulations of 375 ns each, to obtain a cumulative simulation time of 3 μs. With only two data points following after, it is unclear whether this trend would persist further, though we do not expect substantial improvements considering that the values observed at 9 and 10 replicas seem to indicate that a plateau was reached. Based on the totality of the data, we conclude that our DNCF model is adequately converged for the purpose of this study. It should be noted that our analysis constitutes a very rough estimate that is specifically limited to the PDZ2 system, whose allostery does not involve substantial conformational changes.
It has been pointed out that the allosteric residue sets from published computational predictions differ substantially for the PDZ2 system [41], fueling our interest determining how well

PLOS ONE
these models agree with the NMR datasets. However, comparing models based on binary classifications alone can be misleading, since each classification relies on an implicit sensitivity threshold which might differ drastically between models. ROC and PR curves are more suitable for this task since they evaluate prediction performances at all possible thresholds, but require raw prediction scores, which are not always available. Fig 5 shows the ROC and PR curves for the models described above and those for which accompanying literature included the necessary scores. We observed comparably high performances for the DNCF and NMR/ MC [10] models ( Table 2, differences within 1-2%), followed by RRS [50] and REDAN [46]. As the NMR/MC model requires NMR structure data, the DNCF method offers a substantial advantage as the necessary simulations can be based on much more commonly available crystal structures. Thus, although these two methods show comparable accuracy, we expect that the DNCF method can applied to a wider range of systems. We also believe that the method has the potential to show improved results for systems for which induced fit phenomena are important, i.e. for which the conformational ensembles of the apo-and holo-structures differ considerably.

Application of allosteric predictions to the PDZ2 domain
Having established good agreement between DNCF scores and allosteric residues, we investigated the usefulness of these additional features for the biochemical interpretation of our predictions in the PDZ2 structure. Integrating the DNCF scores of the model described above into the structure network (Fig 6A and 6B) reveals two high scoring clusters of residues (clusters I and II). The majority of allosteric residues of the NMR dataset are located in cluster I, which stretches from the top region of the binding pocket towards helix α1 and sheet β1 ( Fig  6B-6D). On the other hand, cluster II encompasses the lower part of the binding pocket surrounding the flexible loop L1 (residues 24-33), including the allosteric residues V26 and V30, furthermore its interaction partners R57, Y36, and finally the C-terminal region.
Comparing these observations to other network scoring methods, the NCF model shows a very similar cluster structure (Fig 7A), whereas for CPLC we observed increased scores for residues located next to the peptide binding groove, e.g. V22, L66, H71, A74, V75 and L78 ( Fig  7B-7D). This can be explained directly by the definition of CPLC (see Algorithms section), which attributes high scores to residues bridging structural modules, e.g. binding grooves. On the other hand, centrality scores for loop L1 (specifically residues 30 to 32) in cluster II are substantially lower than in the timeline-based NCF and DNCF methods, which might be explained by the difficulties of a single structure network to represent the switching contacts of flexible regions. This indicates that centrality methods may fail to account for regions with intrinsic flexibility like the L1 loop, for which methods based on structure ensembles are potentially more appropriate.

Consensus model of allosteric information flow in PDZ2
Finally, we defined a new consensus model of allosteric information flow consolidating our and previous prediction models. For this we first determined a "consensus set" composed of residues predicted as allosteric in � 50% from a selection of published studies (S3 Table) [10, [42][43][44][45][47][48][49][50][51]66]. Next, we obtained a core set of allosteric candidates from our DNCF model, using the score threshold closest to the top left corner in Fig 5B (6.17 bits in S4 Table; TPR: 0.75; FPR: 0.11). This core prediction set (Fig 8 and S5 Table) contains 9 out of 14 residues from the NMR dataset and 11 of the 18 from the consensus set, while 14 residues are complementary predictions. Of these infrequently predicted residues, three form a contiguous surface located on the sheet β1 (F7, V9, L11), connected via L18, V85, and L87 to the peptide binding pocket (Fig 6D). In NMR experiments, V9 was shown to respond to the binding pocket I20F mutation with L11 and L87 as presumed linker residues [40], an interpretation supported by our model. Notably, the clusters surrounding V9 and Y36 agree very well with the DS3 and DS4 regions described previously [10]. Predictions of the C-terminal tail residues  [10,39,66]. Our results suggest the existence of at least two allosteric clusters: Cluster I which encompasses DS1, DS2, and DS3, while cluster II corresponds to DS4.

Discussion
Integration of interaction timelines from molecular dynamics simulations into protein structure networks provides a promising framework for investigating dynamic effects in proteins

PLOS ONE
such as allostery. In this work, we introduce our network analysis tool SenseNet which builds on this theoretical foundation. Using the PDZ2 domain as a reference system, we evaluated four allosteric prediction models implemented in SenseNet, i.e. BC, CPLC, NCF and DNCF, and determined a set of network parameters optimizing their accuracy. Our results are consistent with literature data, as structure networks frequently use carbon contact cutoff distances between 4-6 Å [10, 19,47,67,68], which corresponds approximately to the upper limit of attractive Van-der-Waals interactions. The trend for better prediction results using apo protein states might reflect the observed rigidification of the ligand binding site after binding [39] and is in line with previous suggestions that allosteric mechanisms may be intrinsic properties of apo structures [42,69]. Finally, the improvements observed in sidechain exclusive networks mirror the origins of the NMR dataset, which was obtained from methyl sidechain dynamics [39]. This also highlights an important caveat for comparing prediction models, as some methods might by design match certain types of experimental data more closely than others. Methods based on interaction timelines, i.e. NCF and DNCF, were consistently more accurate than the BC and CPLC methods based on network centrality. This highlights the benefits of using MD simulations to include protein dynamics in protein structure networks, which is achieved by application of methods utilizing interaction timelines. In contrast, centrality-based methods offer the advantage of requiring only a single structure, which makes them uniquely inexpensive in a situation where MD simulations are not feasible. Our data indicate that both BC and CPLC methods could achieve good prediction performances, but were sensitive to the choice of parameters used for network construction. For these in particular, further evaluation studies spanning multiple systems are needed to determine an optimal parameter set that performs well in a wide range of proteins. Of the methods tested, DNCF proved to be the most accurate and robust to changes in network parameters, followed by NCF. This reflects the DNCF method's ability to capture effects from two simulations representing different system states by comparing the changes in shared information. However, the NCF method appears to have potential on its own for predictions based on apo structures alone, for example when there is no known structure of the investigated protein bound to the allosteric ligand. The final allosteric model, based on the DNCF method, was found to be one of the models aligning most closely to experimental data out of those reported in literature, alongside NMR/ MC. However, the DNCF approach offers three distinct advantages to NMR/MC: First, MD simulations for DNCF analyses can be started from only a single, e.g. X-ray, structure, while NMR/MC needs an NMR structure ensemble, which are far rarer and more limited to small proteins. Second, the DNCF method includes all residue types, while NMR/MC by definition cannot predict alanine residues. Third, the DNCF method has the potential to detect induced fit-based conformational changes, which are often not directly detectable in the structural ensembles of the apo-state alone. We determined that 3 μs of total simulation time, spread across 8 replicas and corresponding to 375 ns of simulation for each replica, approximated optimal prediction performance using the DNCF method in the PDZ2 system. These numbers are likely specific to the protein system under investigation and thus can only serve as a guideline for proteins of comparable size and with allosteric effects in the absence of large conformational changes. It should be noted, that fewer replicas and shorter simulation times could still achieve solid performance, which may be relevant when investigating larger proteins for which generating a comparable amount of simulation data may be infeasible. In these cases, additional validation with experimental data is indicated. Our numbers are in agreement with a previous study investigating the reproducibility between replicas in a 10 residue system as well as a 827 residue TCR-p-MHC complex, which recommended using between 5 to 10 replicas for simulations as a rule of thumb [70].
Mapping the results of our DNCF model to the structure of PDZ2 suggests the protein contains two distinct allosteric sites. Most of the experimentally verified allosteric residues from the NMR dataset are located in cluster I, while cluster II has little support from the experimental dataset as the region encompasses only four residues with methyl groups. To fill this gap, alternative experiments may be necessary such as mutational studies connected to changes in PDZ mediated activation. The locations of our observed clusters are matched by several other computational predictions [42,43,45]. Nevertheless, our data contrasts with studies reporting up to four distinct allosteric sites [10,39,66] by suggesting that these four sites are partially overlapping, leaving only two clearly separated allosteric regions. The variance in published allosteric predictions in the PDZ2 domain may be explained by the fact that the experimentally verified data in a single protein are naturally sparse, leading to potentially large error margins for validation. In addition, for many cases quantitative scores are not reported along binary classifications, impeding direct comparison of predictions. To improve prediction models, large scale studies including multiple proteins, computational methods, and experimental data sources will be necessary. With SenseNet we provide a network analysis tool offering considerable advantages over existing implementations: First, by defining edges via interaction timelines, all conformational states of a simulation are readily available for analysis, which is not possible if interactions are reduced to correlation coefficients. Second, adopting a multi-resolution approach via mapping of sub-structures of varying sizes to nodes (from atoms to residues) allows the creation of application-specific network topologies that reduce the underlying structural differences to the most informative level of details. Finally, integration of our tool into Cytoscape allows users to complement their analyses with the community driven ecosystem of biological network analysis plugins, e.g. by connecting structural analysis with system biological or sequence/evolutionary information. Based on these concepts, SenseNet provides an analysis platform implementing a range of well tested analysis algorithms, an easy-to-use UI driven implementation, and interactive side-by-side structure visualization. Together, these features serve as a potential foundation for wide application of timeline-based protein structure networks, paving the way for comparative studies to improve model accuracies and aid experiments in unveiling detailed mechanisms of dynamic processes in biomolecules.  Table. NMR reference set of experimentally verified allosteric and non-allosteric residues. Allosteric residues are represented by a value of 1, non-allosteric residues by a value of 0. (XLSX)