{'fraction_non_alphanumeric': 0.034691059716948566, 'fraction_numerical': 0.02099873432286273, 'mean_word_length': 4.754054948692486, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Effective Molecular Dynamics from Neural-Network Based Structure Prediction Models

Recent breakthroughs in neural-network based structure prediction methods, such as AlphaFold2 and RoseTTAFold, have dramatically improved the quality of computational protein structure prediction. These models also provide statistical confidence scores that can estimate uncertainties in the predicted structures, but it remains unclear to what extent these scores are related to the intrinsic conformational dynamics of proteins. Here we compare AlphaFold2 prediction scores with >60 μs of explicit molecular dynamics simulations of 28 one- and two-domain proteins with varying degree of flexibility. We demonstrate a strong correlation between the statistical prediction scores and the explicit motion derived from extensive atomistic molecular dynamics simulations, and further derive an elastic network model based on the statistical scores of AlphFold2 (AF-ENM), which we benchmark in combination with coarse-grained molecular dynamics simulations. We show that our AF-ENM method reproduces the global protein dynamics with improved accuracy, providing a powerful way to derive effective molecular dynamics using neural-network based structure prediction models.


Introduction
AlphaFold2 1 and RoseTTA fold 2 can infer protein structures based on sequence information with high accuracy. As a result, many projects have emerged that aim to improve the accessibility and usability of these structure prediction tools, most notably the AlphaFold Protein Structure Database 3 , which currently provides predicted structures for over 200 million proteins 4 , ColabFold 5 , which combines fast homology search of MMseqs2 6 with AlphaFold2 and RoseTTA fold, and derivations for to predict large protein complexes 7, 8 .
AlphaFold uses the predicted Local Distance Difference Test (pLDDT) to estimate the accuracy study 42 , leading to benchmarking data set of 61.5 µs aMD simulation. The aMD simulations were performed using the a99SB-disp force field, as it accurately captures both folded and disordered protein states 43 . To assess the degree of flexibility of individual residues, we used the average standard deviation of the Cα-distances to the 20 nearest amino acids (Øσd,20). Unlike the widely used root-mean-square-fluctuation (RMSF) metric, the Øσd,20 does not depend on structural alignment, and is therefore well suited for (highly flexible) multidomain proteins. To evaluate the global protein dynamics, we also used the standard deviation of all Cα distances (σd), which can be represented as a matrix similar to the PAE scores. Overall, groups of residues with higher Øσd,20, and σd values can be interpreted as regions with high conformational flexibility.
A comparison between the PAE (in blue) and σd (in red) matrices, as well as the pLDDT and Øσd,20 scores, for the 28 systems are shown in Fig. 2, with a more detailed analysis of the individual proteins in Fig. S1-S28. We note that high pLDDT scores correspond to low Øσd,20 values, with an overall correlation coefficient R of 0.65, with a 1 Å increase in the Øσd,20 corresponding to a 9 unit decrease of the pLDDT score. The overall linear regression across all tested proteins is: pLDDT = 99.1 − 9.0 • Øσ d,20 [Å] [1] We observe that the R values for the regressions of individual systems vary between 0.24 to 0.99, with a median value of 0.74 (Table S2). While the slopes vary from -1. 3  uncertainties, but rather good correlation to the explicit dynamics within the same system (Table   S3). While the overall correlations are not necessarily high enough for a quantitative assessment of protein dynamics, we later show that a threshold-based approach can reliably differentiate between rigid and flexible regions and be used for constructing ENMs (see below).
Overall, the qualitative consistency of the pLDDT and PAE scores across all studied protein systems is striking. We notice that a PAE score of ~8 Å appears to set a good threshold to visibly identify highly flexible regions as well as dynamically separated domains, corresponding to a σd of around 2.5 Å. For example, for protein systems 1, 2, 3, 21, and 27 (arginine repressor, WW domain-binding protein 4, nuclear factor of activated T-cells cytoplasmic 2, calmodulin, and linear diubiquitin), the PAE threshold correctly predicts the high inter-domain flexibility. AlphaFold also appropriately assesses the increased flexibility of N-terminal (system 7 and 16, sorcin and 30S ribosomal) and C-terminal tails (e.g., system 17, 20, 24, and 26 (cell division protein ZapA, riboflavin synthase, ubiquitin, and green fluorescent protein)), as well as dynamic linker regions (e.g., systems 2, 11 and 28 (WW domain-binding protein 4, chloride intracellular channel protein 1, and NTD-MD construct of Hsp90)). Generally, residues with lower pLDDT scores show higher PAE scores for residue pairs, suggesting that the PAE matrix also encodes the local dynamics, e.g., of flexible loop regions. This close link between the scores might explain some of the observed discrepancies between PAE scores and the σd metric in systems 9, 22, and 23 (elongation factor Ts, peptidyl-prolyl cis-trans isomerase FKBP42, and DNA polymerase beta). The differences could be attributable to the over-or underestimation of the conformational dynamics in a small linker section between two domains, with errors that could propagate and lead to significant changes in the PAE scores.
Based on a closer look at our largest and most complicated system within the dataset, the NTD- Thus, our analysis suggests that the pLDDT scores can be closely linked to local protein dynamics, whereas the PAE scores capture large-scale conformational flexibility, with both scores being closely related to each other. We also find that relatively small errors in the pLDDT scores easily propagate large mismatches between the PAE matrix and global protein dynamics.

Using AlphaFold scores to construct elastic network models
To utilize the AlphaFold scores in explicit molecular dynamics simulations, we propose a framework to convert the statistical pLDDT and PAE scores into an elastic network model (ENM) [23][24][25] . An ENM can be constructed by adding harmonic potentials between, e.g., Cα atoms of the protein backbone of residues i and j according to, where is the total potential energy of the elastic network (en), c is the cutoff radius, and en is the force constant for the harmonic bond. The ENM additionally relies on reference distances 0 inferred from an initial reference structure. In the simplest formulation, a generic constant is used for all en s, although this approach becomes inaccurate for proteins with multiple prominent conformational states. ENMs can also be generated based on short-timescale aMD simulations (in the order of tens of ns) in which both the number as well as the strength of bonds are optimized [49][50][51] . However, many processes, such as the dissociation events and large-scale conformational changes, lie outside the short sampling timescales, making this approach less well suited for highly dynamical systems 49 . Here we propose an alternative strategy to adapt bonds and force constant based on AlphaFold scores.
For our AlphaFold optimized elastic network model (AF-ENM), we suggest scaling the force constants based on the PAE scores. To this end, we assume that the distance distribution between the Cα atoms corresponds to a normal distribution, with the predicted aligned error as the standard deviation. The distance distribution P for a given distance between residues i and j can therefore be defined as, Since the PAE matrix is generally not symmetrical, we use the average value between residue PAEij and PAEji. The relative free energy ΔG(r) can be calculated from the distance distribution, By combining equations [3] and [5], the force constant ken for the bond between residues i and j can be defined as, where kb is the Boltzmann constant and T is the temperature.
NMs are typically derived based on a single structure, but this can result in a significant shows no contacts.
To account for this bias, we propose to use the AlphaFold scores to remove network bonds with high expected fluctuation (Fig. 3). To this end, we used the pLDDT scores to identify residues of the rigid protein cores (low Øσd,20) and flexible/disordered regions (high Øσd,20) that we further used to construct an optimized AF-ENM. We identify a threshold that balances the probability of finding a pLDDT score above a specific value, corresponding to a low Øσd,20 (sensitivity or true positive rate) with the likelihood that the pLDDT score is below the threshold associated with a high Øσd,20 value (specificity or true negative rate). The specificity and sensitivity show an inverse relationship (Fig. 3c). We observe an excellent balance between them for a pLDDT threshold of 90, demonstrating a true positive rate of 86.5% and a true negative rate of 80.1% at a Øσd,20 threshold of 2 Å, as compared to 63.7% and 92.2% for a pLDDT threshold of 95, and 95.2% and 66.6% for a pLDDT threshold of 80.
While the pLDDT scores correlate with the local dynamics, the PAE scores are also good indicators of the global dynamics. We therefore use the PAEs to avoid over-constraining the conformational dynamics between two rigid regions in our AF-ENM. Similar to the dependency between pLDDT scores and Øσd,20, we obtain a good balance between the sensitivity and specificity of the PAE/σd relationship (Fig. 3d), as indicated by the likelihood that a low PAE score corresponds to a low σd (sensitivity or true positive test) relative to the probability that a high PAE value corresponds to a high σd (specificity or true negative test). As for the pLDDT Øσd,20 relationship, we observe that the true positive and negative rates have an inverse relationship. Since the addition of a bond that spuriously restricts the dynamics is far more problematic than the absence of bonds in a rigid structure, we suggest that high sensitivity is more important than high specificity within the AF-ENM model. In this regard, we find that a PAE threshold of 2 Å provides a specificity of 98.8% and a sensitivity of 30.8%, compared to 99.8% and 14.0% for a PAE threshold of 1.5 Å, and 86.6% and 70.1% for a PAE threshold of 4 Å.
Finally, we also probed the effect of different rc on our AF-ENMs. We find that increasing the distance cutoffs correlate with increasing PAE values (Fig. S29), with the likelihood of including a given bond in the AF-ENM strongly depending on the Cα-Cα distances (Fig. 3e). By only using a pLDDT threshold of 90, the acceptance rate is around 82% for Cα-Cα distances between 8-10 Å, and 73% for Cα-Cα distances between 18-20 Å (Fig. 3e). In contrast, when a PAE threshold of 2.0 Å is added, these values decrease to 74% and 24%, while the PAE threshold of 1.5 Å leads to 56% and 1.4%. To assess the extent to which the accepted bonds are representative of the structural ensemble, we examined the probability of the discrepancy of Cα-Cα distances from AlphaFold and aMD simulations by more than 20% (Fig. 3e). We observe that the deviation probability for a cutoff radius of 10 Å changes from 8.1% to 3.2% by including only distances with a pLDDT > 90 and PAE < 2 Å.
Based on these results, we used the following framework to construct AF-ENMs for MARTINI3 cgMD simulations: starting from an initial structure, we only considered residues with pLDDT scores > 90 that are separated by at least two residues since the angle and bonded interactions account for their dynamics. Moreover, to balance sensitivity and selectivity, we used a cutoff radius of 9 Å, with force constants determined by the PAE scores based on equation [6]. For the AF-ENM cgMD simulations, we only consider bonds with > 0.75 kJ/(mol • Å 2 ), corresponding to a PAE threshold of 1.85 Å.

Testing the AF-ENM together with Martini3
To test the AF-ENM approach, we compared MARTINI3-cgMD simulations with a standard ENM using a generic ken (5 kJ/(mol • Å 2 ) and the AF-ENM for three single domain (protein 24-  Table S4). A lower KL divergence corresponds to a higher similarity of two distributions and, therefore, a higher agreement between the global dynamics described by aMD and cgMD simulations.
Overall, we find that using AF-ENM instead of standard ENM leads to lower KL divergence (Table S4). The difference is particularly pronounced for systems with higher conformational dynamics. For small rigid ubiquitin (Fig. S30), we observe a KL divergence of 1.1 for the cgMD simulation with ENM compared to 1.8 with AF-ENM. In the case of lysozyme C (Fig. S31), which has more flexible loop regions, we observe a reduction in KL divergence from 4.5 to 1.4. GFP with the flexible N-terminal helix (Fig. S32) shows a KL divergence of 11.0 and 4.8 for ENM and AF-ENM, a similar reduction as for the larger and flexible NTD and MD domains (Fig. S34), where we notice a decrease from 11.9 to 1.2 and from 9.0 to 4.8. This is also consistent if comparing two domains, where we observe a KL divergence of 9.4 and 5.9 for linear diubiquitin ( Fig. S33) and 10.5 and 4.4 for the combined NTD-MD construct (Fig. S34). This trend appears to be consistent for variations in the elastic network model, as we observe similar KL divergences if we double the force constants for our AF-ENM simulations (Table S4).
The advantage of using the AF-ENM with MARTINI3 cgMD simulations is particularly evident when combined with increased protein-water (P-W) interactions 35-39. In Figure 5,

Discussion
While AlphaFold is primarily used to predict protein structures, it is becoming increasingly clear that it can also be a valuable tool for inferring protein dynamics as shown here (cf. also Refs. 11,14,15,18,19 ). Here we performed a comprehensive comparison of 28 different systems to link AlphaFold confidence scores to metrics obtained by atomistic molecular dynamics simulations.
We show that the pLDDT scores correlate with the local residual dynamics, while the PAE scores also provide insight into large-scale conformational dynamics. Although these scores are suited to unambiguously identify highly flexible and rigid protein regions and inter-domain interactions, accurately describing conditionally folded states remains challenging.
Our quantitative analysis allowed us to construct elastic networks based on these statistical scores (Fig. 6), that can be effectively combined with cgMD simulations. For this, we take advantage of AlphaFold's ability to identify flexible regions between rigid domains. We find that a pLDDT threshold of 90 is sufficient to reliably specify well-defined structural elements, which should be conserved by the network. A PAE threshold between 1.5 to 2.0 Å appears to be well suited to avoid overly restrictive bonds between rigid sections. Since the magnitude of the Cα-Cα distance fluctuation correlates with the PAE score, the PAE matrix can be used to scale the force constants of the individual bonds. Our approach allows us to optimize an ENM with easily accessible or available data without relying on, e.g., additional atomistic simulations [49][50][51] .
While we have here tested AF-ENM in combination with the MARTIN3 coarse-grained force field 29, our approach could also increase the accuracy of other ENM-dependent cgMD models.
Furthermore, our AF-ENM method could be used with small modifications to analytically determine effective protein dynamics via normal mode analysis 59 by primarily relying on the PAE scores to scale the force constants of the harmonic potentials. Although this approach would not be applicable for highly dynamic systems, it could nevertheless provide an improved description over conventional ENMs for rigid proteins. Similar to the AF-ENM framework, the AlphaFold scores could possibly be applied to derive Gō-models 40, 41 instead of elastic networks. Hereby, the harmonic potentials would be replaced by van-der-Waals interactions, whose strength and occurrence would be controlled based on pLDDT and PAE scores. This approach could overcome the observed limitations in conditionally folded structural elements, such as the β8-strand of Hsp90 42 . The AF-ENM framework could thus open up new avenues to incorporate structural information into molecular dynamics simulations.

MD simulations
Atomistic molecular dynamics (aMD) simulations were performed for 28 proteins (Figure 1), summarized in Table S1. The aMD simulations were performed with Gromacs 60 using the a99SBdisp force field 43  Coarse-grained MD (cgMD) simulations (Table S4)      c) Sensitivity and specificity for the pLDDT/ Øσd,20 relations. The blue curves show the probability that a Øσd,20 value is below a cutoff value (represented by different transparency levels) when the corresponding pLDDT score is above a threshold value (x-axis). The red curves show the probability that Øσd,20 is above a cutoff value for residues with a pLDDT score below a threshold value. d) Sensitivity and specificity for the σd/PAE relationship. The blue curves show the probability that a PAE score (x-axis) is below a threshold value and is associated with a σd below a threshold value (represented by different transparency levels). The red curves show the probability that a PAE and associated σd are above threshold values (true negative test). e)