As Simple As Possible, but Not Simpler: Exploring the Fidelity of Coarse-Grained Protein Models for Simulated Force Spectroscopy
[['No other scientific discipline has been so challenged to match the standard of physics-based simplicity as molecular and cell biology , perhaps in parts due to the inherent complexity of the systems under study and to our incomplete knowledge of the structure and function of the living cell .', 'In narrowing this gap , minimal models of proteins have been developed as a step towards the goal of finding an “irreducible element” that still captures at least some of the essential physics and can thus reproduce and predict experimental measurements [1 , 2] .', 'In this regard , minimal models have enjoyed success in testing , refining , and validating the conceptual foundations of the energy landscape theory of protein folding [3–7] as well as forced unfolding mechanisms [8] .', 'A minimal model attempts to capture the essential dynamical behavior of a protein , while upholding the notion of simplicity along with its concommitant computational efficiency .', 'In practice this involves coarse-grained ( CG ) representations of a protein with fewer degrees of freedom than the atomic level of description , simpler , phenomenological interaction potentials , and classical rather than quantum dynamics .', 'Various semi-quantitative comparisons between CG models and experiments have been made [9–11] .', 'At present however , systematic tests comparing the accuracy of coarse-grained models with fully atomistic models are still in need .', 'Fully-atomistic models of proteins have their own shortcomings , including the inability of current atomistic force-fields to fold some proteins such as ubiquitin , a problem which has however been addressed recently and at least partially resolved [12] .', 'However , all-atom models have now been successful in folding small proteins [13 , 14] , elucidating the binding properties of small-molecule drugs [15] , and characterizing complex molecular processes such as ribosomal translation [16] .', 'Steered molecular dynamics ( SMD ) simulations can provide an in silico realization of experimental force microscopy studies [17–19] , where a force can be applied to a single protein– by optical tweezers for example– to unfold it [9 , 10 , 20] .', 'Such computational studies can reveal details of the conformations of proteins during forced unfolding at atomic resolution .', 'Force-extension curves obtained from atomic force microscopy ( AFM ) or optical trap assays generally display a saw-tooth pattern , where each partial unfolding event corresponds to a sudden drop in resistive force [9 , 10 , 20–22] .', 'Our objective in this paper is to evaluate several CG models in SMD simulations by comparing the unfolding mechanisms predicted by each model to those predicted by a reference all-atom simulation under the same conditions .', 'To this end , we construct scaling procedures such that the time , energy , and force scales can be meaningfully compared , and we develop several different metrics that each provide a different viewpoint of the unfolding dynamics .', 'It has been shown that the dynamics of small , globular proteins is well-depicted by all-atom force fields with CHARMM22* with explicit TIP3P water molecules as solvent [14] .', 'Atomistic simulations with explicit solvent , however , are limited in length and time scales of order 100nm and a few μs , unless specialized hardware is used [23] .', 'Simulating the complete unfolding process of a full protein in explicit solvent is currently unfeasible if one wishes to simulate the unfolding mechanism with the same pulling rates as in experiments , and obtain comparable statistics .', 'Thus , to simulate and sample large systems , coarse-grained models are required , because the energy function can be evaluated rapidly and the resulting molecular dynamics does not require a short time step .', 'Various aspects of the protein dynamics and folded structures are successfully captured by structure-based Gō-like models [1 , 24–27] , in which the protein is biased towards its native folded state by native interactions .', 'An interesting question is whether structure-based models can accurately capture the dynamics and the intermediate conformations of partially unfolded proteins during the mechanical unfolding process [9 , 10] .', 'Here we consider three Gō-like models at different levels of resolution: the Associative memory , Water mediated , Structure and Energy Model ( AWSEM-Gō ) [27]; a heavy-atom Gō model [25] that considers all atoms except hydrogen; and a one bead per residue Cα-based Gō model [24] .', 'Several previous studies have compared CG models to all-atom simulations and experiments [9 , 10] .', 'Nevertheless , none of these studies have taken into account that effective time and energy scales must be normalized for meaningful comparison .', 'There is some disagreement whether or not the unfolding pathways predicted by structure based models agree with all-atom simulations or experimental observations [9 , 10] .', 'The authors of ref .', '[10] propose that the unfolding pathway from both CG models of titin I27 domain protein and all-atom implicit solvent simulations are not consistent with the experimental results even at low pulling speeds .', 'On the other hand , CG pulling simulations of T4 lysozyme in Ref .', '[9] qualitatively agree with the experimental findings [28–30] .', 'Sun et al . [31] have compared structure-based Gō models and experiments using force-clamp simulations; these comparisons show general agreement but often fail when sequence details are important in determining the weights of folding intermediates .', 'In this paper , we study the forced unfolding process of a monomer of a loop-truncated variant of superoxide dismutase ( SOD1 ) .', 'SOD1 was the first protein discovered in which mutations had an autosomal-dominant causal relationship to amyotrophic lateral sclerosis ( ALS ) [32 , 33] , an invariably fatal motor neuron degenerative disease characterized by progressive loss of motor neurons [34] , with a lifetime risk by age 70 of about 1/1000 [35] .', 'The loop-truncated variant of SOD1 has loops IV ( residues 49–81 ) and VII ( residues 124–139 ) replaced with short Gly-Ala-Gly tripeptide linkers; here we denote this variant simply as tSOD1 [36 , 37] .', 'tSOD1 consists of a β-barrel tertiary fold containing 8 β-strands and 110 residues .', 'While full-length SOD1 readily forms a homodimer , tSOD1 is obligately monomeric .', 'Moreover , the disulfide bond between C57 and C146 is no longer formed due to the truncation of loop IV and removal of the putative C57 .', 'In experimental protein constructs , the remaining cysteines are mutated ( C6A/C111S/C146S ) to avoid intermolecular crosslinking; we employ the same construct here .', 'In what follows , we first present the details of each model and the simulation set-up .', 'We next describe the normalization of time and energy across models , by calibrating the pulling-rate , temperature , and force in the CG models with respect to the all-atom model .', 'Then we discuss the force-extension curves we obtained , the evolution of structure as the protein is unfolded , and the predictions of the unfolding pathways provided by each model .', 'We finally conclude and briefly discuss the implications of our results .'], ['The aim of this study is to simulate the pulling process of the loop truncated SOD1 protein [37] , and compare the results of an all-atom model with several coarse-grained ( CG ) models .', 'The experimental structure of the tSOD1 monomer can be found as chain A of PDB ID 4BCZ .', 'Force spectroscopy simulations were carried out by tethering both termini with a harmonic potential .', 'The last residue ( C-terminus ) is then moved along the vector from C- to N-terminus with constant velocity of 1 m/s .', 'The stiffness of the spring that imparts the pulling force on the protein was set to 1000 kJ/ ( mol · nm2 ) .', 'Experimental pulling speeds in atomic force microscopy ( AFM ) vary widely between 10−8-10−2 m/s [38–40] , while typical speeds in atomistic simulations are significantly faster , also varying widely between 1–1000 m/s [9 , 10 , 22 , 41] .', 'Simulating and sampling the unfolding mechanism of a full protein in explicit solvent with the same pulling rate as in experiments is currently not feasible .', 'The faster pulling rates in simulations may preempt slow dynamical transitions on the unfolding pathway that would otherwise occur at slower rates .', 'A systematic study of the dependence of the unfolding mechanism on pulling rate for the present system is an interesting topic for future research .', 'Four different types of force fields and protein models were considered: an all-atom ( AA ) simulation in explicit solvent , a heavy atom Gō model ( HA-Gō ) [25] , the Associative memory , Water mediated , Structure and Energy Model ( AWSEM-Gō ) [27] , and a Cα-Gō model [24] in order of decreasing resolution .', 'In the HA-Gō model [25] , all heavy atoms are present .', 'The AWSEM-Gō [27] model is an associative memory Hamiltonian model with a three-bead representation per amino acid .', 'In the Cα-Gō [24] model , each amino acid is represented by only one bead [11 , 24] .', 'Note that in Gō models , only native interactions are attractive , while non-native interactions are purely repulsive .', 'Further description of the Gō model including interaction potentials is given in the specific models sections below .', 'Fig 1 shows a representation of four amino acids in each of the models .', 'Pulling simulations were repeated 20 times for each model , with the same initial structure but different random seeds .', 'To compare the mechanical unfolding pathway of the protein in the all-atom and coarse-grained models , we computed the number of native contacts of all configurations during the pulling simulations .', 'The definition of a native contact is the same throughout this paper .', 'We calculated the native contacts for pairwise distances of all the moieties i and j in each model for any protein structure ( these may be heavy atoms , or coarse-grained residues ) .', 'The fraction of the native contacts Q for conformation X , Q ( X ) , is defined as', 'Q ( X ) = 1 | S | ∑ ( i , j ) ∈ S 1 1 + exp [ β 0 ( r i j ( X ) - λ r i j 0 ) ] , ( 3 )', 'where rij ( X ) is the distance between moieties i and j in conformation X , r i j 0 is the distance between the corresponding moieties i to j in the native state conformation , S is the set of all pairs of native contacts ( i , j ) belonging to the native structure .', 'Amino acids having a native contact must be separated by four or more residues in the primary sequence and r i j 0 < r c u t ( rcut is a model-dependent cutoff distance given in Table 1 ) in the native state [58] , β0 is a smoothing parameter and the factor λ takes into account the fluctuations of the contacts .', 'As a result of adjusting rcut , different models exhibit approximately the same native contact map , and a scatter plot of the number of native contacts present during the pulling trajectory for the AA model vs the CG model exhibits a slope of unity ( y ≈ x ) , ( see Supporting Information S1 Fig ) .', 'Table 1 summarizes the values of β0 , λ , and rcut for each model .', 'For the all-atom and HA-Gō models , the same set of parameters were used as the models share the same structure .', 'The number of contacts for pairs of residues in the Cα-Gō model were weighted with respect to the number of contacts between the same pair in the native state of the protein in the all-atom model , i . e . a given pair of residues could have more than one contact between them , in proportion to how many of their heavy atoms were in contact .', 'All new contacts that are formed during the simulations between moieties i and j are considered non-native contacts if the distance between i , j in the PDB structure is larger than rcut , see Table 1 for values of rcut in each model .', 'To count the total number of non-native contacts in configuration X , we introduce a smooth function that interpolates between 1 and 0 as distance between i and j is increased , with a characteristic length scale R0 given by the mean of the distances between native pairs in the PDB structure: R 0 = 〈 r i j 0 〉 .', 'The smoothing parameter β0 and the factor λ are the same as for native contacts ( see Table 1 ) .', 'The number of non-native contacts in configuration X is then:', 'N n n ( X ) = ∑ ( i , j ) ′ 1 1 + exp [ β 0 ( r i j ( X ) - λ R 0 ) ] , ( 4 ) R0 = 0 . 24 , 0 . 46 , 0 . 91 nm for the AA & HA-Gō , AWSEM , and Cα-Gō model , respectively .', 'We use Q in our analysis of all models as a convenient order parameter on which to project the unfolding mechanism , independent of its accuracy as a kinetic reaction coordinate .', 'In what follows , we will also look at other quantities describing unfolding , such as β sheet dissociation , and structural alignment of remaining parts of the native fold .', 'The interpretation of “time” and “energy” in a CG model must be carefully considered .', 'The energy landscape of CG models is generally smoother , due to softer interaction potentials , reduced degrees of freedom , and lack of explicit solvent molecules .', 'A smoother potential energy surface leads to faster dynamics in comparison to all-atom forcefields .', 'Therefore , the meaning of time in CG models is not the same as in all-atom explicit simulations .', 'When comparing time , velocity , energy , and forces in CG models and all-atom force fields , we should interpret the results with respect to an “effective” energy and time in the system .'], ['Fig 3 depicts representative snapshots of a pulling simulation in the all-atom model .', 'The N- and C-termini are shown in red and blue spheres , respectively , and the structure of the protein is color-coded based on the residue index in the primary sequence .', 'The reported values are the change in separation distance δx = xi − x0 , where x0 and xi are the initial and instantaneous separation distance between tether points respectively ( see Fig 3b ) .', 'As we strain the protein , destabilized contacts between residues break , and regions of secondary structure in the protein are disrupted and dissociate .', 'β-strands lose their native contacts , and locally unfold .', 'The residues in the dissociated regions are then free to form turns or coil structures .', 'In the all-atom simulations , the dissociation of the C-terminus at δx = 4 . 2 nm is the first unfolding event , see Fig 3b .', 'In the unfolding trajectory shown in Fig 3c , we observe the dissociation of part of the N-terminus ( β1-strand ) at δx = 9 . 2 nm .', 'In the snapshot shown in Fig 3f , the β5 and β6 sheets unravel last .', 'At δx = 30 nm , the protein loses all its native contacts and forms a coiled chain .', 'In force spectroscopy simulations , the force ramps up until multiple contacts break , releasing the applied load .', 'We observe multiple force drops ( corresponding to multiple unfolding events ) in the force extension curves .', 'Fig 4 shows a force extension curve for one run of the AA ( black line ) , AWSEM ( red line ) , HA-Gō ( blue line ) , and Cα model ( cyan line ) .', 'Contact maps of the protein averaged over the four runs in Fig 5 are depicted in Fig 6 .', 'In this work , native contacts are defined from the initial PDB structure .', 'The upper triangle shows all native contacts at Q = 0 . 8 , Q = 0 . 5 , and Q = 0 . 1 , respectively from left to right .', 'The bottom triangle shows all non-native contacts , i . e . all new contacts that are formed during the course of the simulation .', 'Since some of these residue pairs may also posess native contacts , they will appear in both maps .', 'It is clear from the figure that native contacts induce the formation of many nearby contacts in the contact map when thermal fluctuations are taken into account .', 'Native contacts between residues k and l are color coded by the thermal average number of contacts divided by the total number of contacts in the PDB structure , 〈Qkl ( Q ) 〉 .', 'Non-native contacts do not have a particular reference structure to normalize with respect to .', 'We thus color code the non-native contact between residues k and l by the frequency of occurrence of any non-native contacts between those residues in the ensemble of structures at Q , i . e . the fraction of conformations at Q that have at least one non-native contact between residues k and l .', 'Here “at Q” means within the bin Q − δQ , Q + δQ , where δQ = 0 . 01 .', 'At Q = 0 . 8 , the native contact maps are approximately the same for all the models ( see first column , Q = 0 . 8 ) .', 'As the protein unfolds from Q = 0 . 8 to Q = 0 . 5 ( second column ) , the C-terminal domain unfolds completely in all the models .', 'The contact maps predict the same general unfolding events until Q ≈ 0 . 5 .', 'As the protein unfolds further to Q = 0 . 1 ( third column ) , the unfolding processes begin to take different pathways across models .', 'For the HA-Gō and Cα-models , the largest folded domain is located at residues 50–70 , while for the AWSEM model the folded domain lies in residues 10–30 .', 'The remaining structured domain in the AA model is larger but only partially folded , consisting of residues 10–60 .', 'We wish to emphasize that Fig 6 is not intended to illustrate the dominant unfolding mechanism for each model , but is simply an analysis of a subset of the unfolding trajectories , chosen only because they were distinct .', 'A further analysis of the dominant unfolding mechanism will be discussed in the subsequent text and corresponding figures .', 'Fig 6 shows that , for all models having more than one interaction site per amino acid , non-native contacts consist largely of what one might call “near-native” contacts .', 'For example in the AA and HA-Gō models , pairs of amino acids have several native contacts between their constituent atoms , however some atom pairs exist between these same amino acids that are not in contact in the native PDB structure .', '“Near-native” contacts would involve these particular atom pairs , and the non-native contact map , which does not have any native interactions by construction , appears quite similar to the native contact map as a result .', 'The presence of native interactions increases the likelihood of proximal non-native interactions .', 'On the other hand , the Cα-Gō model has only one interaction site per amino acid and so cannot exhibit near-native contacts .', 'The non-native contact map is thus sparser than the other models , and involves distinct amino acid pairs .', 'The short-range contacts reminiscent of α-helical structure that are observed at Q = 0 . 1 in the Cα-Gō model are a consequence of the lenient cutoff used for contacts between Cα residues– the other models would show these non-native contacts as well , but because they have more degrees of freedom their cutoff distance for non-native interactions are shorter .', 'Non-native interactions between amino acid pairs wherein one amino acid has been shifted in primary sequence by one , i . e . from amino acids ( m , n ) to ( m ± 1 , n ) or ( m , n ± 1 ) , can be induced by the shear forces between β-strands in the present assay , so that strands may slide over each other or reptate .', 'Similar reptation has been observed in unbiased folding simulations of a β-hairpin [65] .', 'Here , such “off-native” contacts are relatively common for all models that have more than one interaction site per amino acid; the relative numbers of amino acid pairs that partake in off-native contacts compared to the number of amino acid pairs partaking in native contacts , at the values of Q in Fig 6 , are given in Table 4 .', 'To determine the sequence of the unfolding residues , we monitored the number of native contacts for each residue during the pulling simulations .', 'Fig 7A plots the average fraction Qk ( Q ) of a given residue k as a function of total Q for all models .', 'To calculate Qk ( Q ) , we normalize the number of contacts at Q by the number of contacts that residue k possesses in the native structure where Q = 1 .', 'Red color corresponds to Qk ( Q ) = 1 and white indicates Qk ( Q ) = 0 , i . e . the residue has lost all its native contacts .', 'The color scheme in Fig 7B represents the sequence ( in terms of the global order parameter Q ) by which residues lose more than 50% of their contacts during unfolding .', 'The most persistent residues are colored dark blue , and the residues that are broken first in sequence are colored white .', 'From Fig 7 , it is clear that all models predict as first event the dissociation of the C-terminus , residues 100–110 ( β8 ) .', 'Then , in the AA , HA-Gō , and Cα-Gō model , the N-terminus detaches .', 'The average unfolding pathways predicted by the AA model are very similar to the HA-Gō model , where residues in the N- and C-terminus dissociate first , and the contacts of residues 50–74 are broken last .', 'In contrast , the sequence of unfolding in the AWSEM model starts from β8 and β7 , and the last domains to rip off are β3 and β2 .', 'In summary , the similarity between unfolding events depicted in Fig 7B may be quantified by computing the correlation coefficient between the degree of remaining structure for individual β strands ( the similarity of the darkness of the bands for each model in Fig 7B ) .', 'This gives the following correlation coefficients: between AA and HA-Gō: 0 . 94 , between AA and Cα-Gō: 0 . 86 , and between AA and AWSEM: 0 . 62 .', 'In order to determine whether there exists a well defined unfolding pathway of the tSOD1 protein , and if so , to compare it across models , we used the template modeling score ( TM-score ) [66] to compare the similarity between the protein structures of different pulling trajectories at the same Q . The TM-score for the alignment of two structures is defined as [66]:', 'TM=1L∑N11+ ( did ) 2 , d=1 . 24 ( N−15 ) 3−1 . 8 , ( 8 )', 'where N is the number of residue pairs , di is the distance between identical residues i in two structures , and L is the number of residues in the reference structure .', 'The TM-score lies between 0 and 1; a TM-score of one indicates that the two protein structures are perfectly matched .', 'Usually , two structures with TM-score higher than 0 . 5 are considered to have the same folded conformations , while uncorrelated protein structures have a TM-score of less than 0 . 2 [66] .', 'Measuring the TM-alignment , as well as clustering of structures by TM-score , was performed by using Maxcluster ( http://www . sbg . bio . ic . ac . uk/maxcluster ) [67] .', 'TM-scores of an all-against-all structure comparison of folded segment of protein structures obtained from each run for Q = 0 . 8 , Q = 0 . 4 , and Q = 0 . 2 are shown in Fig 8 . The color code quantifies the TM-score of pairs of structures at the same value of Q , obtained from all pairs of trajectories: red color indicates perfectly matched structures , and white represents a TM-score of zero .', 'For comparing the conformations , we only considered Cα-atoms in the backbone for the folded region of the protein .', 'This folded region at each Q-value was defined as a contiguous sequence of n residues with residue index i ≤ j ≤ i + n , where 〈Qi ( Q ) 〉 > 0 . 5 and 〈Qi+n ( Q ) 〉 > 0 . 5 .', 'The average here corresponds to the ensemble of states of all trajectories .', 'If there is an unfolded region with more than 10 residues in between i and i + n , then the largest contiguous sequence of residues with 〈Qi ( Q ) 〉 > 0 . 5 was considered .', 'In Fig 8 , TM-scores for Q = 0 . 8 ( see left column in Fig 8 ) are high for all four models , which indicates that at the beginning of the unfolding process , the backbone of the protein is very similar in the unfolding trajectories .', 'The Cα-model and the HA-Gō model exhibit slightly larger deviations between trajectories at this value of Q . As the protein unfolds further , at Q = 0 . 4 ( second column in Fig 8 ) , the TM-scores drop to lower values .', 'In the AA model , the average TM-score of one trajectory ( run 20 ) is 0 . 33 , while other runs have higher TM-scores .', 'For the HA-Gō model , values of the TM-score range between 0 . 3–0 . 6 .', 'In the Cα-Gō model , the TM-scores range between 0 . 5–0 . 76 .', 'At the same Q = 0 . 4 , the TM-scores in the AWSEM model are still much higher and vary between 0 . 6–0 . 94 , which indicates the presence of one dominant pathway .', 'It is clear from the large number of trajectories with high TM-scores that the AWSEM model exhibits a much stronger pathway behavior than the other models , which begin to balkanize into clusters of residual structure .', 'This can also be clearly seen by plotting the mean TM-score between all M ( M − 1 ) /2 trajectories ( M = 20 here ) as a function of Q , for all four models , see Fig 9 . At Q = 0 . 2 , the TM-scores for AA , HA-Gō , and Cα-Gō models have reached about 0 . 2 , which is comparable to the TM-score of a random coil ensemble .', 'This indicates a highly diverse residual structure between trajectories .', 'The length of the residual folded structures at Q = 0 . 2 is only about 24 , 38 , 21 , and 27 for AA , HA-Gō , AWSEM , and Cα-Gō models .', 'Thus , the AA , HA-Gō , and Cα-Gō models predict multiple unfolding pathways for lower values of Q . On the other hand , the AWSEM model still has a fairly high TM-score; indicating that it predicts only one main unfolding pathway .', 'Two structures that are nearly folded at Q ≈ 1 are obliged to have a high TM score , while two structures at low Q are not so obliged .', 'We thus also plot in Fig 9 a reference curve to compare the structural overlap .', 'We construct this curve by taking a window containing a given number of residues ( e . g . 50 ) , and slide this window along all possible locations of the primary sequence ( 1–50 , 2–51 , etc . ) , to obtain a set of partial native structures , one structure for each window position .', 'This process is repeated for all window sequence lengths .', 'The native contacts Q are calculated for all of the structures , binned , and TM-aligned .', 'This gives a randomized set of partially unfolded structures , which nevertheless lack thermal fluctuations and strain distortions , and so would tend to have larger TM-alignments when they overlap .', 'Interestingly , this curve lies roughly between the AWSEM model and all other models , consistent with the strong pathway-like unfolding mechanism of the AWSEM model .', 'In order to more clearly render the unfolding pathways predicted by each model , we clustered the protein conformations based on the TM-scores during the unfolding at several different Q-values , see Figs 10 and 11 .', 'The structures shown are centroids of the corresponding clusters that emerge from the clustering analysis .', 'A TM-score cut-off of 0 . 6 is used to define when configurations no longer belong to a given cluster .', 'The coloring is based on the residue index , where the C-terminus of the structured protein is in red and the N-terminus is colored blue .', 'The thickness of the lines is proportional to the fraction of total trajectories in each cluster .', 'As can be seen in Figs 10 and 11 , each model predicts a dominant unfolding route , which is shown with a thick black line .', 'All models predict one unique unfolding pathway until Q ≈ 0 . 44 .', 'Along this pathway , β strand 8 at the C-terminus loses structure first , however subsequent events differ between models .', 'As the structure continues to unfold from Q ≈ 0 . 44 to Q ≈ 0 . 2 , we observe multiple unfolding pathways in all models but the AWSEM model; see Fig 10 panel', 'a ) AA ,', 'b ) HA-Gō , and Fig 11b ) Cα-Gō models .', 'The protein structures from different pulling simulations in the above 3 models are distributed in multiple diverse conformations .', 'For the AA , HA , and Cα-Gō models , β strand 1 on the N-terminus generally dissociates after β strand 8 at the C-terminus .', 'In 3 out of 20 trajectories of the AA model however , β strands 1 and 2 were the last to unfold .', 'This mechanism with β strands 1 and 2 unfolding last is the pathway observed in the AWSEM model .', 'Generally , the last unfolding events involve breakage of contacts in β strands 5 and 6 in the AA model .', 'The sequence of unfolding events along the main forced unfolding pathway in the AA model is β strand 8 , then β1 and 7 , β2 , then β3 and 4 , β6 , and then finally β5 .', 'In the HA-Gō model , the sequence of unfolding of events is β8 and β1 , then β2 , β7 , then β3 and 4 , then β6 and finally β5 is the last domain to unfold , which is similar to the AA model .', 'In the Cα-Gō , the first unfolding event is also dissociation of C-terminal β strand 8 , then β1 , β strands 2 and 7 , then β3 , β4 , β6 , and finally β5 .', 'In contrast to the above three models , the AWSEM model ( Fig 11c ) predicts only one unfolding pathway .', 'In this pathway , the unfolding of the protein starts from the C-terminal β strand 8 , then β7 , β4 , the C-terminal portion constituting roughly half of β strand 3 , the N-terminal portion constituting roughly half of β strand 1 , β strands 5 and 6 and the remainder of β strand 3 , the remainder of β1 , and β2 .', 'Strands 1 and 2 were the last to dissociate in all the 20 trajectories .', 'We also compare the main pathway of unfolding of the AA-model with other models by calculating the TM-score between the AA model and the three CG models .', 'For comparison across different models , TM-score was calculated using the program TM-align [68] .', 'The conformations of the most populated cluster at Q in the AA model was compared to the corresponding conformations in the other models at the same value of Q . In order to compare CG with AA models , the TM-alignment only includes the Cα atoms in the backbone of the folded segment of the protein as described above .', 'The TM-score versus Q , for pairs of two models , AA with HA-Gō ( black line ) , AA with AWSEM ( red line ) , and AA with Cα-Gō ( blue line ) , is depicted in Fig 12a .', 'A TM-score with a value of > 0 . 5 for a pair of proteins means that the structures are similar [68] .', 'The observed high TM-scores between AA and all CG models for Q > 0 . 45 indicate that all CG models predict unfolding pathways similar to the AA model by this metric .', 'Interestingly , in the range of Q between 0 . 45 and 1 , the AWSEM model shows the best agreement with the AA model , and the HA-Gō model shows the least agreement .', 'As the protein is unfolded below Q ≤ 0 . 44 , the TM-score shows a more sensitive dependence upon models .', 'At Q less than about 0 . 25 , the TM-scores have reached small values that would be expected for the alignment of random dissimilar structures .', 'We conclude therefore that all models predict similar unfolding pathways until the protein is about half unfolded , at which point the mechanisms begin to diverge from the AA model .', 'The AWSEM model does not predict multiple pathways as the other models do , but the dominant pathway observed for the AWSEM model is structurally as similar to the AA model as any of the other CG models .', 'None of the CG models can completely capture the unfolding mechanism at the lower values of Q for the AA model .', 'The above conclusion is recapitulated by analyzing the corresponding alignment between models using the more conventional metric of RMSD .', 'Comparing the folded core of the AA model in the most populated cluster , as defined in Section “Residue Contacts” , to the same region in the CG models ( most populated cluster , same sequence length as in the AA model ) yields a plot of RMSD vs . Q , as shown in Fig 12b .', 'By this metric , the AWSEM model again shows the best structural alignment ( lowest RMSD ) until Q ≈ 0 . 3 , while the HA-Gō model shows the worst structural alignment .', 'In this paper we explored the limits of validity of several structural-based coarse-grained ( CG ) models by comparing the unfolding mechanisms of a truncated variant of superoxide dismutase , when the protein is subjected to force-induced unfolding .', 'An all-atom ( AA ) , explicit-solvent model is used as the benchmark standard to which the other models are compared .', 'A more desirable comparison would be with experimental data , however no experimental data exists for this particular system , and moreover the data that does exist for other systems does not have the atomic resolution that we have measured and compared with here .', 'Unfortunately then such a comparison is not possible at present .', 'One may entertain the 
['Introduction', 'Methods', 'Results/Discussion']
['Mechanical unfolding of a single domain of loop-truncated superoxide dismutase protein has been simulated via force spectroscopy techniques with both all-atom ( AA ) models and several coarse-grained models having different levels of resolution: A Gō model containing all heavy atoms in the protein ( HA-Gō ) , the associative memory , water mediated , structure and energy model ( AWSEM ) which has 3 interaction sites per amino acid , and a Gō model containing only one interaction site per amino acid at the Cα position ( Cα-Gō ) .', 'To systematically compare results across models , the scales of time , energy , and force had to be suitably renormalized in each model .', 'Surprisingly , the HA-Gō model gives the softest protein , exhibiting much smaller force peaks than all other models after the above renormalization .', 'Clustering to render a structural taxonomy as the protein unfolds showed that the AA , HA-Gō , and Cα-Gō models exhibit a single pathway for early unfolding , which eventually bifurcates repeatedly to multiple branches only after the protein is about half-unfolded .', 'The AWSEM model shows a single dominant unfolding pathway over the whole range of unfolding , in contrast to all other models .', 'TM alignment , clustering analysis , and native contact maps show that the AWSEM pathway has however the most structural similarity to the AA model at high nativeness , but the least structural similarity to the AA model at low nativeness .', 'In comparison to the AA model , the sequence of native contact breakage is best predicted by the HA-Gō model .', 'All models consistently predict a similar unfolding mechanism for early force-induced unfolding events , but diverge in their predictions for late stage unfolding events when the protein is more significantly disordered .']
['Although experimentalists can now unfold single proteins in the lab by pulling them apart and measuring the force and extension , a clear idea of how the protein changes shape and loses structure during this process is currently missing .', 'Molecular dynamics simulations can offer insight as to what is actually happening structurally when you pull a protein apart .', 'However , typical simulations of processes that happen nearly instantaneously in the lab take weeks to perform , when every atom must be accounted for .', 'Researchers have thus resorted to much faster “coarse-grained models” , where the system is simplified by removing select atoms and the remaining interactions rescaled , but the accuracy of such simulations are known to suffer as a result .', 'How accurate or inaccurate are the current coarse-grained models in capturing the unfolding mechanisms of proteins ?', 'Our findings upon investigating this question suggest that , while coarse-grained models successfully capture early unfolding events of nearly-folded proteins , they suffer when trying to describe the late stages of unfolding in mostly-disordered proteins .', 'By showing how coarse-grained models may fail to capture the accuracy of their more sophisticated but cumbersome counterparts , we can shed light on how to improve their reliability , increase their speed , and enhance their relevance in capturing biologically-relevant phenomena .']
['enzymes', 'condensed matter physics', 'enzymology', 'dismutases', 'simulation and modeling', 'protein structure prediction', 'protein structure', 'research and analysis methods', 'relaxation (physics)', 'proteins', 'melting', 'phase transitions', 'molecular biology', 'physics', 'protein structure comparison', 'biochemistry', 'biochemical simulations', 'superoxide dismutase', 'biology and life sciences', 'physical sciences', 'computational biology', 'macromolecular structure analysis']