Protein folding is an important subject not only for basic research in molecular biology but also for understanding folding diseases and designing new polymeric materials (Dill and MacCallum, 2012). Transient, partially folded states are often encountered on folding pathways, and have been characterized experimentally in solution by methods such as laser temperature jumps, fluorescence labeling, and solution X-ray scattering. Mutagenesis evaluated by Φ-value analysis (Fersht and Daggett, 2002), for instance, has also provided residue-level information on transition states. Recently, state-of-the-art single-molecule (sm) measurements, single-molecule Förster resonance energy transfer (smFRET) (Chung et al., 2012; Oikawa et al., 2013) and force spectroscopy (smFS) (Neupane et al., 2016) have become powerful tools in protein-folding research, providing reliable information on transition-path (barrier-crossing process) times (Chung et al., 2012) and heterogeneity in the unfolded state (Oikawa et al., 2013). A major limitation of smFRET is that the observables are restricted to ‘low-dimensional’ structural data, such as the donor–acceptor distance. Computational modeling should help us to interpret single-molecule time-series measurements and should contribute to solving some remaining puzzles, such as the reduced solvent viscosity dependence of the transition-path times (Chung and Eaton, 2013), the internal viscosity (Chung and Eaton, 2013) and the impact of non-Markovian property (Chung et al., 2015). Theories and computational methods have been developed to extract structural information and dynamics from smFRET data (Haas et al., 2013; Hoefling et al., 2011; Matsunaga et al., 2015; Sun et al., 2016).

Molecular dynamics (MD) simulation is another powerful approach for investigating protein dynamics and folding over relatively long time periods — hundreds of microseconds or longer (Lindorff-Larsen et al., 2011). In theory, ‘low-dimensional’ smFRET measurements are interpreted in terms of the atomic structural models obtained with MD simulations. However, it is still a challenge to achieve quantitative agreement between simulation and experimental data during the entire folding process, owing not only to simulation time limitations but also to inherent force-field biases. In particular, while local interactions are well described by the current force fields, it is still difficult to reproduce energetic balances between unfolded and folded states. Indeed, Piana and coworkers showed, in their protein folding simulation study, that the folding mechanism of the villin headpiece depends substantially on the choice of force field (Piana et al., 2011). It is also known that most force fields produce unfolded states that are more compact and structured than those suggested experimentally (Piana et al., 2014). Methods based on maximum entropy (Beauchamp et al., 2014; Boomsma et al., 2014; Cavalli et al., 2013; Olsson et al., 2013, 2014, 2017; Pitera and Chodera, 2012; Roux and Weare, 2013) or Bayesian statistics (Bonomi et al., 2016) were recently developed for guiding simulations or models to generate ensembles that match experimental ensemble-averaged observations.

Exploiting time-series data from single-molecule experiments is another way to link simulation with experiment and has several advantages over the ensemble-average based approaches: (i) more latent states can be uncovered by inferring states from their historical evolution than from their static ensembles (Li et al., 2008; Matsunaga et al., 2015; Schuetz et al., 2010); (ii) the transition state can be uniquely identified as a dynamic bottleneck by following the actual dynamics. However, the time-scale gap between experiment and simulation previously hampered the direct use of time-series analysis methods in other disciplines, such as data assimilations that are based on the sequential Monte Carlo method (Matsunaga et al., 2015).

Here, we develop a new approach for single-molecule time-series data based on the Markov state model (MSM) (Pande et al., 2010), a statistical model that approximates dynamics by memory-less probabilistic transitions between discrete conformational states. In MSM, the probability of transition from the discrete state i to state j in a lag time of τ is described by a transition-probability matrix $T\left(\tau \right)=\left\{T_{ij}\left(\tau \right)\right\}$. $T\left(\tau \right)$ can be estimated from a set of short simulations instead of a single long simulation. Therefore, MSM is often used to generate long-time dynamics (Silva et al., 2014) whose time-scale is comparable to those of experimental measurements (Feng et al., 2016; Noé et al., 2011; Snow et al., 2002). However, estimation of $T\left(\tau \right)$ is largely dependent on the simulation force-fields, which may have uncertainties or biases for certain conformations (Olsson et al., 2017). To overcome this problem, we propose to "refine" $T\left(\tau \right)$ using high-resolution measurement of single-molecules as time-series data. Specifically, we use a machine-learning approach to estimate $T\left(\tau \right)$ between hidden Markov states from low-dimensional time-series data (Bishop, 2006). To distinguish the original $T\left(\tau \right)\left(=T_{\text{simulation}}\left(\tau \right)\right)$, we refer to the refined $T\left(\tau \right)$ as $T_{\text{experiment}}\left(\tau \right)$ in this paper.

We propose a two-step procedure in our machine-learning approach, which links simulations and single-molecule experiments (Figure 1). (i) Supervised learning step. We first construct an initial MSM from a raw set of simulation data. After defining discrete conformational states by clustering the trajectory snapshots, $T\left(\tau \right)=T_{\text{simulation}}\left(\tau \right)$ is estimated directly by counting transitions between the discrete states in the simulation trajectories. This step is the same as conventional MSM (Pande et al., 2010). (ii) Unsupervised learning step. Using $T_{\text{simulation}}\left(\tau \right)$ as an initial estimate, we perform hidden Markov modeling (Bronson et al., 2009; McKinney et al., 2006; Okamoto and Sako, 2012; Pirchi et al., 2016; Rabiner and Juang, 1986; Schröder and Grubmüller, 2003) to refine the initial MSM using single-molecule measurement time-series data. $T\left(\tau \right)$ is optimized so that the "refined" or data-assimilated MSM with $T_{\text{experiment}}\left(\tau \right)$ can reproduce the time-series data most accurately.

Figure 1 with 1 supplement see all

Download asset Open asset

We applied this procedure to the folding dynamics of the formin-binding protein (FBP) WW domain, a 37-residue three-stranded β-sheet protein. In the construction of the initial MSM, extensive MD simulations of a dye-labeled WW domain were performed for an aggregated time of ~400 μs. Although there are a number of folding simulation studies for the FBP WW domain and its homologs (Ensign and Pande, 2009; Freddolino et al., 2008; Karanicolas and Brooks, 2003; Zanetti-Polzi et al., 2017; Zhou et al., 2014), this may be the first folding simulation to use FRET dyes. High time-resolution smFRET measurements of WW domain folding and unfolding dynamics were used for the unsupervised learning (Chung et al., 2012). The data were previously measured to resolve the durations of the folding and unfolding transitions with microsecond resolution by Chung and coworkers (Chung et al., 2012). The initial and data-assimilated MSM showed different folding pathways and transition-state ensembles of WW domains. An independent mutation experiment with Φ-value analysis (Fersht and Daggett, 2002) validated the data-assimilated MSM. We discuss our time-series analysis method in the context of machine-learning theory and its applicability to conformational transitions in other biomolecules.

We have proposed a two-step procedure for the construction of a data-assimilated MSM with $T\left(\tau \right)=T_{\text{experiment}}\left(\tau \right)$ matching single-molecule time-series data. Using smFRET data for the FBP WW domain, we show that the data-assimilated MSM successfully reproduces the original smFRET data, and yields a transition-state ensemble consistent with an independent mutational experiment (Petrovich et al., 2006). The folding mechanism based on the data-assimilated MSM suggests an interplay between hairpin and hydrophobic formations.

In the context of machine-learning theory, the proposed two-step procedure can be regarded as a semi-supervised learning algorithm, which tries to learn from both labeled and unlabeled data (Rudzinski et al., 2016; Zhu and Goldberg, 2009). In the context of MSM, simulation data correspond to labeled data while experimental data are unlabeled data. In a typical case of semi-supervised learning (e.g., image recognition), the labeled data are correct and usually expensive (e.g., images that are manually labeled by investigator). Therefore, unlabeled data are often ‘de-emphasized’ by scaling their contribution in the likelihood function (Zhu and Goldberg, 2009). By contrast, in our case, labeled data (simulation) may have incorrect transition counts caused by force-field biases whereas unlabeled data (experimental) possess more reliable information on dynamics. Thus, in our two-step procedure, the estimates $T_{\text{simulation}}\left(\tau \right)$ from labeled data (simulation) are replaced with $T_{\text{experiment}}\left(\tau \right)$ "refined" with unlabeled data (experiment). This is regarded as the limit of ‘de-emphasis’ on labeled data (simulation).

When fitting a rather complex model to any experimental data, the model can overreact to noise in the data (the overfitting problem). In particular, our MSM for the FBP WW domain has a rather large number of parameters ($87\times 87$ transition probabilities), which could be easily overfitted to the smFRET data. To assess the overfitting in the unsupervised learning, we divided the smFRET data in half, and unsupervised learning was independently applied to the two subsets. Qualitatively similar network structures appeared in both, and were similar to that obtained with the full data set (Figure 2—figure supplement 4). This is because the effective number of parameters is considerably reduced from $87\times 87$ down to only those involving populated states (ε ~ 0.5–0.8). In order to counteract overfitting, a maximum caliber approach for minimally perturbing the initial MSM could be a promising direction for a future study (Dixit and Dill, 2014, 2015, 2018; Wan et al., 2016; Zhou et al., 2017). Furthermore, to see the dependence on the choice of the Förster radius R₀, we carried out the unsupervised learning using a set of different R₀ values (R₀ = 54, 55, 57, and 58 Å). The overall structure of the MSM network was robust against the choice of R₀ except for R₀ = 58 Å (Figure 2—figure supplement 5).

The initial condition in the unsupervised learning is another issue. Since MSM has a larger number of parameters than in typical hidden Markov modeling, unsupervised learning requires a good initial condition for optimization to avoid being trapped in a local minimum. By using $T_{\text{simulation}}\left(\tau \right)$ as the initial condition, we achieved a likelihood convergence with $\ln L\left(T\left(\tau \right)\right)$ = −584,947 with 10,000-step optimization (taking one week using the parallel implementation of the Baum-Welch algorithm [Rabiner and Juang, 1986]). For comparison, we also performed the optimization using a random matrix as the initial condition (Figure 2—figure supplement 8C). In the figure, the optimization of the likelihood looks stacked after 10,000 iterations and its value is lower than that of $T_{\text{simulation}}\left(\tau \right)$ as the initial condition. This result suggests that the model from a random matrix could be ruled out because of trapping in a local minimum, or at least that global optimization of a random matrix is practically very inefficient. This also emphasizes the importance of transferring knowledge $T_{\text{simulation}}\left(\tau \right)$ learned from simulations for improving the unsupervised learning $T_{\text{experiment}}\left(\tau \right)$ from experimental data (Torrey and Jude, 2009). Although just $T_{\text{simulation}}\left(\tau \right)$ was used as the initial condition here, advanced algorithms in transfer learning (Torrey and Jude, 2009) can be incorporated in a future study.

We here used the conventional constant-time binning for photon counting because the standard MSM is based on constant-time binning. A promising possibility for future studies is to apply continuous-time Markov modeling (McGibbon and Pande, 2015), which may allow us to use information-based binning (Watkins and Yang, 2004) or photon-by-photon analysis (Gopich and Szabo, 2009; Okamoto and Sako, 2012; Pirchi et al., 2016), avoiding photon counting noise.

In conclusion, exploiting the temporal information embedded in experimental time-series data to improve the simulation-based model has provided a rich, dynamic and experimentally consistent picture of the folding mechanism for the FBP WW domain. The data-assimilated MSM pathway could be used to improve the force-field parameters of proteins, nucleic acids, and other biomolecules. The semi-supervised learning combined with MSM method developed here is a quite general framework that can be used to understand conformational transitions in proteins and other biomolecules. It can be extended to interpret other experimental data possibly using more advanced techniques.

All of the in-house program codes for the method proposed in this study, including test datasets, are available at GitHub: https://github.com/ymatsunaga/mdtoolbox.

1. 1. Beauchamp KA
   2. Pande VS
   3. Das R

   (2014) Bayesian energy landscape tilting: towards concordant models of molecular ensembles

   Biophysical Journal 106:1381–1390.

   https://doi.org/10.1016/j.bpj.2014.02.009

   • PubMed
   • Google Scholar
2. 1. Berendsen HJC
   2. Postma JPM
   3. van Gunsteren WF
   4. DiNola A
   5. Haak JR

   (1984) Molecular dynamics with coupling to an external bath

   The Journal of Chemical Physics 81:3684–3690.

   https://doi.org/10.1063/1.448118

   • Google Scholar
3. 1. Best RB
   2. Hofmann H
   3. Nettels D
   4. Schuler B

   (2015) Quantitative interpretation of FRET experiments via molecular simulation: force field and validation

   Biophysical Journal 108:2721–2731.

   https://doi.org/10.1016/j.bpj.2015.04.038

   • PubMed
   • Google Scholar
4. 1. Best RB
   2. Hummer G
   3. Eaton WA

   (2013) Native contacts determine protein folding mechanisms in atomistic simulations

   PNAS 110:17874–17879.

   https://doi.org/10.1073/pnas.1311599110

   • PubMed
   • Google Scholar
5. 1. Best RB
   2. Zheng W
   3. Mittal J

   (2014) Balanced protein-water interactions improve properties of disordered proteins and non-specific protein association

   Journal of Chemical Theory and Computation 10:5113–5124.

   https://doi.org/10.1021/ct500569b

   • PubMed
   • Google Scholar
6. Book

   1. Bishop CM

   (2006)

   Pattern Recognition and Machine Learning

   Springer.

   • Google Scholar
7. 1. Bonomi M
   2. Camilloni C
   3. Cavalli A
   4. Vendruscolo M

   (2016) Metainference: A Bayesian inference method for heterogeneous systems

   Science Advances 2:e1501177.

   https://doi.org/10.1126/sciadv.1501177

   • PubMed
   • Google Scholar
8. 1. Boomsma W
   2. Ferkinghoff-Borg J
   3. Lindorff-Larsen K

   (2014) Combining experiments and simulations using the maximum entropy principle

   PLoS Computational Biology 10:e1003406.

   https://doi.org/10.1371/journal.pcbi.1003406

   • PubMed
   • Google Scholar
9. 1. Bronson JE
   2. Fei J
   3. Hofman JM
   4. Gonzalez RL
   5. Wiggins CH

   (2009) Learning rates and states from biophysical time series: a Bayesian approach to model selection and single-molecule FRET data

   Biophysical Journal 97:3196–3205.

   https://doi.org/10.1016/j.bpj.2009.09.031

   • PubMed
   • Google Scholar
10. 1. Camley BA
    2. Brown FLH
    3. Lipman EA

    (2009) Förster transfer outside the weak-excitation limit

    The Journal of Chemical Physics 131:104509.

    https://doi.org/10.1063/1.3230974

    • Google Scholar
11. 1. Cavalli A
    2. Camilloni C
    3. Vendruscolo M

    (2013) Molecular dynamics simulations with replica-averaged structural restraints generate structural ensembles according to the maximum entropy principle

    The Journal of Chemical Physics 138:094112.

    https://doi.org/10.1063/1.4793625

    • PubMed
    • Google Scholar
12. 1. Chung HS
    2. Eaton WA

    (2013) Single-molecule fluorescence probes dynamics of barrier crossing

    Nature 502:685–688.

    https://doi.org/10.1038/nature12649

    • PubMed
    • Google Scholar
13. 1. Chung HS
    2. McHale K
    3. Louis JM
    4. Eaton WA

    (2012) Single-molecule fluorescence experiments determine protein folding transition path times

    Science 335:981–984.

    https://doi.org/10.1126/science.1215768

    • PubMed
    • Google Scholar
14. 1. Chung HS
    2. Piana-Agostinetti S
    3. Shaw DE
    4. Eaton WA

    (2015) Structural origin of slow diffusion in protein folding

    Science 349:1504–1510.

    https://doi.org/10.1126/science.aab1369

    • PubMed
    • Google Scholar
15. 1. Darden T
    2. York D
    3. Pedersen L

    (1993) Particle mesh Ewald: An N ⋅log( N ) method for Ewald sums in large systems

    The Journal of Chemical Physics 98:10089–10092.

    https://doi.org/10.1063/1.464397

    • Google Scholar
16. 1. Dill KA
    2. MacCallum JL

    (2012) The protein-folding problem, 50 years on

    Science 338:1042–1046.

    https://doi.org/10.1126/science.1219021

    • PubMed
    • Google Scholar
17. 1. Dinner AR
    2. Lazaridis T
    3. Karplus M

    (1999) Understanding beta-hairpin formation

    PNAS 96:9068–9073.

    https://doi.org/10.1073/pnas.96.16.9068

    • PubMed
    • Google Scholar
18. 1. Dixit PD
    2. Dill KA

    (2014) Inferring microscopic kinetic rates from stationary state distributions

    Journal of Chemical Theory and Computation 10:3002–3005.

    https://doi.org/10.1021/ct5001389

    • PubMed
    • Google Scholar
19. 1. Dixit PD
    2. Dill KA

    (2018) Caliber corrected markov modeling (c₂m₂): Correcting equilibrium markov models

    Journal of Chemical Theory and Computation 14:1111–1119.

    https://doi.org/10.1021/acs.jctc.7b01126

    • PubMed
    • Google Scholar
20. 1. Dixit PD
    2. Jain A
    3. Stock G
    4. Dill KA

    (2015) Inferring transition rates of networks from populations in continuous-time markov processes

    Journal of Chemical Theory and Computation 11:5464–5472.

    https://doi.org/10.1021/acs.jctc.5b00537

    • Google Scholar
21. 1. Ensign DL
    2. Pande VS

    (2009) The Fip35 WW domain folds with structural and mechanistic heterogeneity in molecular dynamics simulations

    Biophysical Journal 96:L53–L55.

    https://doi.org/10.1016/j.bpj.2009.01.024

    • PubMed
    • Google Scholar
22. 1. Feng Y
    2. Zhang L
    3. Wu S
    4. Liu Z
    5. Gao X
    6. Zhang X
    7. Liu M
    8. Liu J
    9. Huang X
    10. Wang W

    (2016) Conformational dynamics of apo-glnbp revealed by experimental and computational analysis

    Angewandte Chemie International Edition 55:13990–13994.

    https://doi.org/10.1002/anie.201606613

    • PubMed
    • Google Scholar
23. 1. Fersht AR
    2. Daggett V

    (2002) Protein folding and unfolding at atomic resolution

    Cell 108:573–582.

    https://doi.org/10.1016/S0092-8674(02)00620-7

    • PubMed
    • Google Scholar
24. 1. Freddolino PL
    2. Liu F
    3. Gruebele M
    4. Schulten K

    (2008) Ten-microsecond molecular dynamics simulation of a fast-folding WW domain

    Biophysical Journal 94:L75–L77.

    https://doi.org/10.1529/biophysj.108.131565

    • PubMed
    • Google Scholar
25. 1. Gopich I
    2. Szabo A

    (2012) Theory of single-molecule FRET efficiency histograms

    Advances in Chemical Physics 146:245–297.

    https://doi.org/10.1002/9781118131374.ch10

    • Google Scholar
26. 1. Gopich IV
    2. Szabo A

    (2009) Decoding the pattern of photon colors in single-molecule FRET

    The Journal of Physical Chemistry B 113:10965–10973.

    https://doi.org/10.1021/jp903671p

    • PubMed
    • Google Scholar
27. 1. Graen T
    2. Hoefling M
    3. Grubmüller H

    (2014) AMBER-DYES: Characterization of charge fluctuations and force field parameterization of fluorescent dyes for molecular dynamics simulations

    Journal of Chemical Theory and Computation 10:5505–5512.

    https://doi.org/10.1021/ct500869p

    • PubMed
    • Google Scholar
28. 1. Haas KR
    2. Yang H
    3. Chu JW

    (2013) Expectation-maximization of the potential of mean force and diffusion coefficient in Langevin dynamics from single molecule FRET data photon by photon

    The Journal of Physical Chemistry B 117:15591–15605.

    https://doi.org/10.1021/jp405983d

    • PubMed
    • Google Scholar
29. 1. Harrigan MP
    2. Sultan MM
    3. Hernández CX
    4. Husic BE
    5. Eastman P
    6. Schwantes CR
    7. Beauchamp KA
    8. McGibbon RT
    9. Pande VS

    (2017) MSMBuilder: Statistical models for biomolecular dynamics

    Biophysical Journal 112:10–15.

    https://doi.org/10.1016/j.bpj.2016.10.042

    • PubMed
    • Google Scholar
30. 1. Hoefling M
    2. Lima N
    3. Haenni D
    4. Seidel CA
    5. Schuler B
    6. Grubmüller H

    (2011) Structural heterogeneity and quantitative FRET efficiency distributions of polyprolines through a hybrid atomistic simulation and Monte Carlo approach

    PLoS ONE 6:e19791.

    https://doi.org/10.1371/journal.pone.0019791

    • PubMed
    • Google Scholar
31. 1. Hopkins CW
    2. Le Grand S
    3. Walker RC
    4. Roitberg AE

    (2015) Long-Time-Step Molecular Dynamics through Hydrogen Mass Repartitioning

    Journal of Chemical Theory and Computation 11:1864–1874.

    https://doi.org/10.1021/ct5010406

    • PubMed
    • Google Scholar
32. 1. Hornak V
    2. Abel R
    3. Okur A
    4. Strockbine B
    5. Roitberg A
    6. Simmerling C

    (2006) Comparison of multiple Amber force fields and development of improved protein backbone parameters

    Proteins: Structure, Function, and Bioinformatics 65:712–725.

    https://doi.org/10.1002/prot.21123

    • PubMed
    • Google Scholar
33. 1. Humphrey W
    2. Dalke A
    3. Schulten K

    (1996) VMD: visual molecular dynamics

    Journal of Molecular Graphics 14:33–38.

    https://doi.org/10.1016/0263-7855(96)00018-5

    • PubMed
    • Google Scholar
34. 1. Jung J
    2. Mori T
    3. Kobayashi C
    4. Matsunaga Y
    5. Yoda T
    6. Feig M
    7. Sugita Y

    (2015) GENESIS: a hybrid-parallel and multi-scale molecular dynamics simulator with enhanced sampling algorithms for biomolecular and cellular simulations

    Wiley Interdisciplinary Reviews: Computational Molecular Science 5:310–323.

    https://doi.org/10.1002/wcms.1220

    • PubMed
    • Google Scholar
35. 1. Jäger M
    2. Nguyen H
    3. Crane JC
    4. Kelly JW
    5. Gruebele M

    (2001) The folding mechanism of a beta-sheet: the WW domain

    Journal of Molecular Biology 311:373–393.

    https://doi.org/10.1006/jmbi.2001.4873

    • PubMed
    • Google Scholar
36. 1. Jäger M
    2. Nir E
    3. Weiss S

    (2006) Site-specific labeling of proteins for single-molecule FRET by combining chemical and enzymatic modification

    Protein Science 15:640–646.

    https://doi.org/10.1110/ps.051851506

    • PubMed
    • Google Scholar
37. 1. Karanicolas J
    2. Brooks CL

    (2003) The structural basis for biphasic kinetics in the folding of the WW domain from a formin-binding protein: lessons for protein design?

    PNAS 100:3954–3959.

    https://doi.org/10.1073/pnas.0731771100

    • PubMed
    • Google Scholar
38. 1. Kellogg EH
    2. Lange OF
    3. Baker D

    (2012) Evaluation and optimization of discrete state models of protein folding

    The Journal of Physical Chemistry B 116:11405–11413.

    https://doi.org/10.1021/jp3044303

    • PubMed
    • Google Scholar
39. 1. Kobayashi C
    2. Jung J
    3. Matsunaga Y
    4. Mori T
    5. Ando T
    6. Tamura K
    7. Kamiya M
    8. Sugita Y

    (2017) GENESIS 1.1: A hybrid-parallel molecular dynamics simulator with enhanced sampling algorithms on multiple computational platforms

    Journal of Computational Chemistry 38:2193–2206.

    https://doi.org/10.1002/jcc.24874

    • Google Scholar
40. 1. Le Grand S
    2. Götz AW
    3. Walker RC

    (2013) SPFP: Speed without compromise—A mixed precision model for GPU accelerated molecular dynamics simulations

    Computer Physics Communications 184:374–380.

    https://doi.org/10.1016/j.cpc.2012.09.022

    • Google Scholar
41. 1. Li CB
    2. Yang H
    3. Komatsuzaki T

    (2008) Multiscale complex network of protein conformational fluctuations in single-molecule time series

    PNAS 105:536–541.

    https://doi.org/10.1073/pnas.0707378105

    • PubMed
    • Google Scholar
42. 1. Lindorff-Larsen K
    2. Piana S
    3. Dror RO
    4. Shaw DE

    (2011) How fast-folding proteins fold

    Science 334:517–520.

    https://doi.org/10.1126/science.1208351

    • PubMed
    • Google Scholar
43. 1. Macias MJ
    2. Gervais V
    3. Civera C
    4. Oschkinat H

    (2000) Structural analysis of WW domains and design of a WW prototype

    Nature Structural Biology 7:375–379.

    https://doi.org/10.1038/75144

    • PubMed
    • Google Scholar
44. 1. Mahoney MW
    2. Jorgensen WL

    (2001) Diffusion constant of the TIP5P model of liquid water

    The Journal of Chemical Physics 114:363–366.

    https://doi.org/10.1063/1.1329346

    • Google Scholar
45. 1. Matsunaga Y
    2. Kidera A
    3. Sugita Y

    (2015) Sequential data assimilation for single-molecule FRET photon-counting data

    The Journal of Chemical Physics 142:214115.

    https://doi.org/10.1063/1.4921983

    • Google Scholar
46. Software

    1. Matsunaga Y

    (2018) mdtoolbox, version f8eed18

    GitHub.

    https://github.com/ymatsunaga/mdtoolbox
47. 1. McGibbon RT
    2. Pande VS

    (2015) Efficient maximum likelihood parameterization of continuous-time Markov processes

    The Journal of Chemical Physics 143:034109.

    https://doi.org/10.1063/1.4926516

    • PubMed
    • Google Scholar
48. Conference

    1. McGibbon RT
    2. Ramsundar B
    3. Sultan MM
    4. Kiss G
    5. Pande VS

    (2014a)

    Understanding Protein Dynamics with L1-Regularized Reversible Hidden Markov Models

    Paper Presented at the Proc. 31st Intl. Conf. on Machine Learning (ICML).

    • Google Scholar
49. 1. McGibbon RT
    2. Schwantes CR
    3. Pande VS

    (2014) Statistical model selection for Markov models of biomolecular dynamics

    The Journal of Physical Chemistry B 118:6475–6481.

    https://doi.org/10.1021/jp411822r

    • PubMed
    • Google Scholar
50. 1. McKinney SA
    2. Joo C
    3. Ha T

    (2006) Analysis of single-molecule FRET trajectories using hidden Markov modeling

    Biophysical Journal 91:1941–1951.

    https://doi.org/10.1529/biophysj.106.082487

    • PubMed
    • Google Scholar
51. 1. Metzner P
    2. Schütte C
    3. Vanden-Eijnden E

    (2009) Transition path theory for markov jump processes

    Multiscale Modeling & Simulation 7:1192–1219.

    https://doi.org/10.1137/070699500

    • Google Scholar
52. 1. Miyamoto S
    2. Kollman PA

    (1992) Settle: An analytical version of the SHAKE and RATTLE algorithm for rigid water models

    Journal of Computational Chemistry 13:952–962.

    https://doi.org/10.1002/jcc.540130805

    • Google Scholar
53. 1. Mu Y
    2. Nordenskiöld L
    3. Tam JP

    (2006) Folding, misfolding, and amyloid protofibril formation of WW domain FBP28

    Biophysical Journal 90:3983–3992.

    https://doi.org/10.1529/biophysj.105.076406

    • PubMed
    • Google Scholar
54. 1. Muñoz V
    2. Thompson PA
    3. Hofrichter J
    4. Eaton WA

    (1997) Folding dynamics and mechanism of beta-hairpin formation

    Nature 390:196–199.

    https://doi.org/10.1038/36626

    • PubMed
    • Google Scholar
55. 1. Neupane K
    2. Manuel AP
    3. Woodside MT

    (2016) Protein folding trajectories can be described quantitatively by one-dimensional diffusion over measured energy landscapes

    Nature Physics 12:700–703.

    https://doi.org/10.1038/nphys3677

    • Google Scholar
56. 1. Noé F
    2. Doose S
    3. Daidone I
    4. Löllmann M
    5. Sauer M
    6. Chodera JD
    7. Smith JC

    (2011) Dynamical fingerprints for probing individual relaxation processes in biomolecular dynamics with simulations and kinetic experiments

    PNAS 108:4822–4827.

    https://doi.org/10.1073/pnas.1004646108

    • PubMed
    • Google Scholar
57. 1. Noé F
    2. Schütte C
    3. Vanden-Eijnden E
    4. Reich L
    5. Weikl TR

    (2009) Constructing the equilibrium ensemble of folding pathways from short off-equilibrium simulations

    PNAS 106:19011–19016.

    https://doi.org/10.1073/pnas.0905466106

    • PubMed
    • Google Scholar
58. 1. Noé F
    2. Wu H
    3. Prinz JH
    4. Plattner N

    (2013) Projected and hidden Markov models for calculating kinetics and metastable states of complex molecules

    The Journal of Chemical Physics 139:184114.

    https://doi.org/10.1063/1.4828816

    • PubMed
    • Google Scholar
59. 1. Oikawa H
    2. Suzuki Y
    3. Saito M
    4. Kamagata K
    5. Arai M
    6. Takahashi S

    (2013) Microsecond dynamics of an unfolded protein by a line confocal tracking of single molecule fluorescence

    Scientific Reports 3:2151.

    https://doi.org/10.1038/srep02151

    • PubMed
    • Google Scholar
60. 1. Okamoto K
    2. Sako Y

    (2012) Variational Bayes analysis of a photon-based hidden Markov model for single-molecule FRET trajectories

    Biophysical Journal 103:1315–1324.

    https://doi.org/10.1016/j.bpj.2012.07.047

    • PubMed
    • Google Scholar
61. 1. Olsson S
    2. Frellsen J
    3. Boomsma W
    4. Mardia KV
    5. Hamelryck T

    (2013) Inference of structure ensembles of flexible biomolecules from sparse, averaged data

    PLoS ONE 8:e79439.

    https://doi.org/10.1371/journal.pone.0079439

    • PubMed
    • Google Scholar
62. 1. Olsson S
    2. Vögeli BR
    3. Cavalli A
    4. Boomsma W
    5. Ferkinghoff-Borg J
    6. Lindorff-Larsen K
    7. Hamelryck T

    (2014) Probabilistic determination of native state ensembles of proteins

    Journal of Chemical Theory and Computation 10:3484–3491.

    https://doi.org/10.1021/ct5001236

    • PubMed
    • Google Scholar
63. 1. Olsson S
    2. Wu H
    3. Paul F
    4. Clementi C
    5. Noé F

    (2017) Combining experimental and simulation data of molecular processes via augmented Markov models

    PNAS 114:8265–8270.

    https://doi.org/10.1073/pnas.1704803114

    • PubMed
    • Google Scholar
64. 1. Pande VS
    2. Beauchamp K
    3. Bowman GR

    (2010) Everything you wanted to know about Markov State Models but were afraid to ask

    Methods 52:99–105.

    https://doi.org/10.1016/j.ymeth.2010.06.002

    • PubMed
    • Google Scholar
65. 1. Petrovich M
    2. Jonsson AL
    3. Ferguson N
    4. Daggett V
    5. Fersht AR

    (2006) Phi-analysis at the experimental limits: mechanism of beta-hairpin formation

    Journal of Molecular Biology 360:865–881.

    https://doi.org/10.1016/j.jmb.2006.05.050

    • PubMed
    • Google Scholar
66. 1. Piana S
    2. Donchev AG
    3. Robustelli P
    4. Shaw DE

    (2015) Water dispersion interactions strongly influence simulated structural properties of disordered protein states

    The Journal of Physical Chemistry B 119:5113–5123.

    https://doi.org/10.1021/jp508971m

    • PubMed
    • Google Scholar
67. 1. Piana S
    2. Klepeis JL
    3. Shaw DE

    (2014) Assessing the accuracy of physical models used in protein-folding simulations: quantitative evidence from long molecular dynamics simulations

    Current Opinion in Structural Biology 24:98–105.

    https://doi.org/10.1016/j.sbi.2013.12.006

    • PubMed
    • Google Scholar
68. 1. Piana S
    2. Lindorff-Larsen K
    3. Shaw DE

    (2011) How robust are protein folding simulations with respect to force field parameterization?

    Biophysical Journal 100:L47–L49.

    https://doi.org/10.1016/j.bpj.2011.03.051

    • PubMed
    • Google Scholar
69. 1. Pirchi M
    2. Tsukanov R
    3. Khamis R
    4. Tomov TE
    5. Berger Y
    6. Khara DC
    7. Volkov H
    8. Haran G
    9. Nir E

    (2016) Photon-by-photon hidden markov model analysis for microsecond single-molecule fret kinetics

    The Journal of Physical Chemistry B 120:13065–13075.

    https://doi.org/10.1021/acs.jpcb.6b10726

    • PubMed
    • Google Scholar
70. 1. Pitera JW
    2. Chodera JD

    (2012) On the use of experimental observations to bias simulated ensembles

    Journal of Chemical Theory and Computation 8:3445–3451.

    https://doi.org/10.1021/ct300112v

    • Google Scholar
71. 1. Prinz JH
    2. Wu H
    3. Sarich M
    4. Keller B
    5. Senne M
    6. Held M
    7. Chodera JD
    8. Schütte C
    9. Noé F

    (2011) Markov models of molecular kinetics: generation and validation

    The Journal of Chemical Physics 134:174105.

    https://doi.org/10.1063/1.3565032

    • PubMed
    • Google Scholar
72. 1. Rabiner L
    2. Juang B

    (1986) An introduction to hidden Markov models

    IEEE ASSP Magazine 3:4–16.

    https://doi.org/10.1109/MASSP.1986.1165342

    • Google Scholar
73. 1. Roux B
    2. Weare J

    (2013) On the statistical equivalence of restrained-ensemble simulations with the maximum entropy method

    The Journal of Chemical Physics 138:084107.

    https://doi.org/10.1063/1.4792208

    • PubMed
    • Google Scholar
74. 1. Rudzinski JF
    2. Kremer K
    3. Bereau T

    (2016) Communication: Consistent interpretation of molecular simulation kinetics using Markov state models biased with external information

    The Journal of Chemical Physics 144:051102.

    https://doi.org/10.1063/1.4941455

    • PubMed
    • Google Scholar
75. 1. Ryckaert J-P
    2. Ciccotti G
    3. Berendsen HJC

    (1977) Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes

    Journal of Computational Physics 23:327–341.

    https://doi.org/10.1016/0021-9991(77)90098-5

    • Google Scholar
76. 1. Salomon-Ferrer R
    2. Götz AW
    3. Poole D
    4. Le Grand S
    5. Walker RC

    (2013) Routine microsecond molecular dynamics simulations with AMBER on GPUs. 2. Explicit Solvent Particle Mesh Ewald

    Journal of Chemical Theory and Computation 9:3878–3888.

    https://doi.org/10.1021/ct400314y

    • PubMed
    • Google Scholar
77. 1. Schröder GF
    2. Grubmüller H

    (2003) Maximum likelihood trajectories from single molecule fluorescence resonance energy transfer experiments

    The Journal of Chemical Physics 119:9920–9924.

    https://doi.org/10.1063/1.1616511

    • Google Scholar
78. 1. Schuetz P
    2. Wuttke R
    3. Schuler B
    4. Caflisch A

    (2010) Free energy surfaces from single-distance information

    The Journal of Physical Chemistry B 114:15227–15235.

    https://doi.org/10.1021/jp1053698

    • PubMed
    • Google Scholar
79. 1. Schwantes CR
    2. Pande VS

    (2013) Improvements in markov state model construction reveal many non-native interactions in the folding of NTL9

    Journal of Chemical Theory and Computation 9:2000–2009.

    https://doi.org/10.1021/ct300878a

    • PubMed
    • Google Scholar
80. 1. Senne M
    2. Trendelkamp-Schroer B
    3. Mey AS
    4. Schütte C
    5. Noé F

    (2012) EMMA: A Software Package for Markov Model Building and Analysis

    Journal of Chemical Theory and Computation 8:2223–2238.

    https://doi.org/10.1021/ct300274u

    • PubMed
    • Google Scholar
81. 1. Silva DA
    2. Weiss DR
    3. Pardo Avila F
    4. Da LT
    5. Levitt M
    6. Wang D
    7. Huang X

    (2014) Millisecond dynamics of RNA polymerase II translocation at atomic resolution

    PNAS 111:7665–7670.

    https://doi.org/10.1073/pnas.1315751111

    • PubMed
    • Google Scholar
82. 1. Snow CD
    2. Nguyen H
    3. Pande VS
    4. Gruebele M

    (2002) Absolute comparison of simulated and experimental protein-folding dynamics

    Nature 420:102–106.

    https://doi.org/10.1038/nature01160

    • PubMed
    • Google Scholar
83. 1. Sun X
    2. Morrell TE
    3. Yang H

    (2016) Extraction of protein conformational modes from distance distributions using structurally imputed bayesian data augmentation

    The Journal of Physical Chemistry B 120:10469–10482.

    https://doi.org/10.1021/acs.jpcb.6b07767

    • PubMed
    • Google Scholar
84. Book

    1. Torrey L
    2. Jude S

    (2009)

    Transfer Learning

    In: Soria E, Martin J, Magdalena R, Martinez M, Serrano A, editors. Handbook of Research on Machine Learning Applications and Trends: Algorithms, Methods, and Techniques 1. IGI Global. pp. 242–264.

    • Google Scholar
85. 1. Wako H
    2. Saitô N

    (1978) Statistical mechanical theory of the protein conformation. Ii. Folding pathway for protein

    Journal of the Physical Society of Japan 44:1939–1945.

    https://doi.org/10.1143/JPSJ.44.1939

    • Google Scholar
86. 1. Wan H
    2. Zhou G
    3. Voelz VA

    (2016) A maximum-caliber approach to predicting perturbed folding kinetics due to mutations

    Journal of Chemical Theory and Computation 12:5768–5776.

    https://doi.org/10.1021/acs.jctc.6b00938

    • PubMed
    • Google Scholar
87. 1. Watkins LP
    2. Yang H

    (2004) Information bounds and optimal analysis of dynamic single molecule measurements

    Biophysical Journal 86:4015–4029.

    https://doi.org/10.1529/biophysj.103.037739

    • PubMed
    • Google Scholar
88. 1. Zanetti-Polzi L
    2. Davis CM
    3. Gruebele M
    4. Dyer RB
    5. Amadei A
    6. Daidone I

    (2017) Parallel folding pathways of Fip35 WW domain explained by infrared spectra and their computer simulation

    FEBS Letters 591:3265–3275.

    https://doi.org/10.1002/1873-3468.12836

    • PubMed
    • Google Scholar
89. 1. Zheng W
    2. Borgia A
    3. Buholzer K
    4. Grishaev A
    5. Schuler B
    6. Best RB

    (2016) Probing the Action of Chemical Denaturant on an Intrinsically Disordered Protein by Simulation and Experiment

    Journal of the American Chemical Society 138:11702–11713.

    https://doi.org/10.1021/jacs.6b05443

    • PubMed
    • Google Scholar
90. 1. Zhou G
    2. Pantelopulos GA
    3. Mukherjee S
    4. Voelz VA

    (2017) Bridging microscopic and macroscopic mechanisms of p53-MDM2 binding with kinetic network models

    Biophysical Journal 113:785–793.

    https://doi.org/10.1016/j.bpj.2017.07.009

    • PubMed
    • Google Scholar
91. 1. Zhou R
    2. Maisuradze GG
    3. Suñol D
    4. Todorovski T
    5. Macias MJ
    6. Xiao Y
    7. Scheraga HA
    8. Czaplewski C
    9. Liwo A

    (2014) Folding kinetics of WW domains with the united residue force field for bridging microscopic motions and experimental measurements

    PNAS 111:18243–18248.

    https://doi.org/10.1073/pnas.1420914111

    • PubMed
    • Google Scholar
92. 1. Zhu X
    2. Goldberg AB

    (2009) Introduction to semi-supervised learning

    Synthesis Lectures on Artificial Intelligence and Machine Learning 3:1–130.

    https://doi.org/10.2200/S00196ED1V01Y200906AIM006

    • Google Scholar
93. 1. Ziv G
    2. Thirumalai D
    3. Haran G

    (2009) Collapse transition in proteins

    Physical Chemistry Chemical Physics 11:83–93.

    https://doi.org/10.1039/B813961J

    • PubMed
    • Google Scholar

Thank you for submitting your article "Linking time series of single-molecule experiments with molecular dynamics simulations by machine learning" for consideration by eLife. Your article has been reviewed by three peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Arup Chakraborty as the Senior Editor. The following individual involved in review of your submission has agreed to reveal his identity: John Straub (Reviewer #2).

The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.

Summary:

The manuscript describes a novel approach based on hidden Markov models to connect single-molecular FRET measurements with MD simulations, and application to the folding of formin-binding protein WW domain. The manuscript address two key issues in improving the validity of Markov state models used in protein folding studies, namely, 1) the construction of the discrete states, and 2) the construction of the continuous time Markov transition matrix T(\tau). Initially, discrete conformational states were obtained by clustering trajectory snapshots, and transition rates were then estimated from trajectories. To correct potential bias due to problems in the force field, the authors then used a hidden Markov model to model smFRET time series measurement and a) modified the discrete states and b) estimated the T(\tau) values between the discrete states so the smFRET is better explained. Indeed, the refined Markov state model can better account for the single-molecule time-series data: The improvement in reproducing the smFRET histogram of 50 uS as shown in Figure 3 is dramatic. It further revealed a different folding pathway of the WW-domain, as well as a different transition state ensemble.

The approach of using MD-derived MSM as basis and initial input for likelihood optimization with respect to smFRET data is very interesting and novel, and the calculations and uncertainty analysis were conducted carefully. The insight into the folding mechanism is very interesting. However, quite a few outstanding issues, including some serious concerns, remain to be addressed to assess whether publication is warranted.

Essential revisions:

1) An overall question about the current practice of Markov state models in MD-based protein folding studies. Assuming both simulation and experimental studies are under the same conditions (salt, dye, temperature, etc.), they are probing the same underlying physical process. Therefore, they should reflect the same ground-truth reality. As MD can provide atomistic details within the time scale of simulation, MD studies should provide far more details than smFRET. However, this is not the case. Rather, MD simulation cannot stand on its own, and needs to be corrected by smFRET. As described by the authors, the MSM is shifted before and after incorporating smFRET data. The authors showed that the discrete states are very different when smFRET data are brought into consideration. While it is not explicitly stated (please clarify), it seems that the refined Markov states are derived from the original Markov states by clustering of MD trajectories, but updated with new equilibrium probability (hence different bubble areas in Figure 2) and new transition rates. How can the ground truth, for example, equilibrium probabilities of the Markov states, change if it is examined one way (by MD) or another way (MD+smFRET)? This would suggest there is no objective ground truth that is accessible by current MD simulations.

Lending support to this conclusion, is the discussion you provided, that in other studies prior to this manuscript, the number of discrete states and the nature of these states are all malleable, depending on which metric one examines (RMSD or contact map). This would again suggest that there is no objective truth in determining what constitutes a set of discrete states that are Markovian and what their transition rates are. Rather, all depends on which metric and/or which additional experimental data one chooses to examine.

One cannot help but speculate: would one expect some other alternative discrete Markovian states different from the currently reported ones to emerge, when other types of experimental time series data other than \epsilon is incorporated. Is it likely that yet more different Markovian states, different dynamics of transition rates, different folding dynamics, and different transition state ensembles will be identified, all for the same WW-domain under the same condition (dye, pH, salt, temperature, if well controlled)?

This aspect of shifting ground truth in MD simulation and Markov state modeling is rather unsettling. A pessimistic, but not unreasonable, view is that without a principled approach in defining true reaction coordinates for protein folding one has to make do with rather ad hoc and heuristic approaches in defining the discrete states with unexamined consequences, and hence this unsettling aspect may be with us for a long while and the truth may be elusive.

Here is a suggestion in this respect: The authors employ Equation 5 for a model of FRET efficiency, however, more detailed models exist such as http://aip.scitation.org/doi/pdf/10.1063/1.3230974

The authors could evaluate their approach and determine if a more detailed model might be justified or needed in making contact with experimental measurements.

2) While concern 1 is more general, there is another question more specific to connecting MD and the smFRET measurements. It will be necessary for the authors to show how well the HMM learned T(\tau) conform with MD simulation generated T(\tau) after the refined discrete states are obtained, in a statistically significant way. For example, the authors should compare the HMM rates and the MD rates for the subsets of state-state transitions with the highest and relevant lowest rates among the 87x87 parameters. In addition, the flux and pathway analysis should also be carried out using MD-derived rates after the refined states are obtained by the HMM model, for key paths carrying most of the flux, and compared with the analysis using HMM rates. The authors may have already done some of these in the transition state ensemble analysis, and may already have all necessary data.

This is important for validating the connection between the sm-FRET (HMM modeling) time series and the MD trajectories the authors are proposing. It would be fantastic if all works out, as this would naturally suggest that in the future one no longer needs to do MD simulation. Rather, one could use smFRET to fit dynamic rates, as long as the discrete states can be obtained, for example, by some other means. However, one should be cautious about being overly optimistic here, as there are significant problems in the enormous search space of very high dimension. It will be very interesting if the authors can show that the HMM rates and MD rates after refinement are the same.

3) Potential over-fitting remains a concern, as there is a large number (87x87 = 7,569) parameters that need to be fitted. While the authors qualitatively compared results using each half of the data through plotting, it would be necessary to do some cross-validation tests, which are standard in machine learning. Specifically, the authors could use (n-1)/n balanced smFRET and MD data (e.g. similar number of folding trajectories) to identify refined states, and estimating \tau values, then test whether the remaining 1/n data falls within these defined states, and whether the unseen epsilon histograms can be reproduced accurately.

4) In Equation 5, the average values are used in calculating \epsilon. It seems that kappa has little difference at different regions of Q. Is the difference along the dimension of \epsilon that is used for clustering due to difference in distances r(s_i)? It will be helpful to also show both mean and variance of r(s_i) for each of the major discrete states, and justify whether the degree of in-cluster and between-cluster heterogeneity/homogeneity conforms to the Markovian assumption of each discrete state. This relates to the issue of validating the defined discrete states.

5) The calculations of FRET efficiencies from MD trajectories need to be stated explicitly. Furthermore, the estimation of FRET efficiencies from experimentally measured photon sequences also needs to be presently clearly. The latter in fact is a long-standing issue in the single-molecule community, for example, the time binning convention and the information-base binning by Haw Yang.

6) It is not obvious how the folding pathways were decomposed into individual pathways for both MSM models. How this procedure was conducted would critically impact the interpretation of results. This part thus needs careful and clear explanation.

7) Certain key references of interring dynamics from smFRET trajectories were missed in the context of this work, such as Haas et al., 2013.

8) Obtaining a solution of similar likelihood value starting from a random matrix is a bit worrisome. Can the authors elaborate more on how this model was ruled out? Furthermore, the convergence of likelihood optimization should also be discussed to reveal the robustness of this approach.

[Editors' note: further revisions were requested prior to acceptance, as described below.]

Thank you for resubmitting your work entitled "Linking time series of single-molecule experiments with molecular dynamics simulations by machine learning" for further consideration at eLife. Your revised article has been favorably evaluated by Arup Chakraborty (Senior Editor), a Reviewing Editor, and three reviewers.

The manuscript has been improved but there are some remaining issues that need to be addressed before acceptance, as outlined below:

Overview:

The unsettling aspects of shifting ground truth, or rather, the lack of ground truth remains. This is, however, an issue the community of protein folding simulations has to deal with.

Using Q as an adequate reaction coordinate is similarly problematic and inadequate, as Q is no substitute to reaction coordinates in the strict sense of earlier classic studies of Chandler and Dinner. But this again is a community problem, and we do not wish to specifically penalize the authors of the current work. Therefore, we accept the authors' arguments, although at least one of the reviewers had some fundamental concerns about the entire set of approaches in this area.

Specific points to address:

1) There is one issue that the authors did not provide answers to but they can easily do, namely, regarding the question of providing details of the 87x87 rates from the Hidden Markov model and the corresponding rates from MD simulation. It will be necessary to make first a general statement on how similar/different these two sets of matched/paired rates are, and second, provide details, e.g., in the Appendix, with side-by-side or overlaid histograms. If these two distributions are similar, it would be a comforting result. However, even if these two distributions do not agree, it will be important for readers to know this fact so they can directly draw their own inferences. This can be easily done with existing data.

2) A grammar and expression overhaul is highly recommended to enhance the readability of the work.

Essential revisions:

1) An overall question about the current practice of Markov state models in MD-based protein folding studies. Assuming both simulation and experimental studies are under the same conditions (salt, dye, temperature, etc.), they are probing the same underlying physical process. Therefore, they should reflect the same ground-truth reality. As MD can provide atomistic details within the time scale of simulation, MD studies should provide far more details than smFRET. However, this is not the case. Rather, MD simulation cannot stand on its own, and needs to be corrected by smFRET. As described by the authors, the MSM is shifted before and after incorporating smFRET data. The authors showed that the discrete states are very different when smFRET data are brought into consideration. While it is not explicitly stated (please clarify), it seems that the refined Markov states are derived from the original Markov states by clustering of MD trajectories, but updated with new equilibrium probability (hence different bubble areas in Figure 2) and new transition rates. How can the ground truth, for example, equilibrium probabilities of the Markov states, change if it is examined one way (by MD) or another way (MD+smFRET)? This would suggest there is no objective ground truth that is accessible by current MD simulations.

We agree with the reviewers that both simulation and experimental data should reflect the same ground-truth reality. However, despite improvements over the decades, force field parameters used in MD simulation are not perfectly correct as reaching the ground-truth reality. In particular, while the local interactions are well described by the current force fields it is still difficult to reproduce the energetic balance between global states (such as unfolded and folded states). Indeed, Piana and coworkers (Piana et al., 2011) showed, in their protein folding simulation study, the folding mechanism of the villin headpiece depends substantially on the choice of force field. It is also known that most force fields produce unfolded states that are more compact and structured than those suggested experimentally (Piana et al., 2014). As suggested by these studies, a subtle free energy balance between the unfolded and folded states is crucial for determining the folding mechanism. Considering this background, the motivation of this work is to make MD data toward the ground-truth reality in a data-driven way. We reweighted the transition probabilities by using smFRET data to achieve the subtle balance in free energy. In the future study, using the free energy landscape and folding pathways obtained by our scheme, it is possible to improve force field parameters closer to the ground truth reality. According to this discussion, we have revised the Introduction section (second paragraph).

Lending support to this conclusion, is the discussion you provided, that in other studies prior to this manuscript, the number of discrete states and the nature of these states are all malleable, depending on which metric one examines (RMSD or contact map). This would again suggest that there is no objective truth in determining what constitutes a set of discrete states that are Markovian and what their transition rates are. Rather, all depends on which metric and/or which additional experimental data one chooses to examine.

One cannot help but speculate: would one expect some other alternative discrete Markovian states different from the currently reported ones to emerge, when other types of experimental time series data other than \epsilon is incorporated. Is it likely that yet more different Markovian states, different dynamics of transition rates, different folding dynamics, and different transition state ensembles will be identified, all for the same WW-domain under the same condition (dye, pH, salt, temperature, if well controlled)?

This aspect of shifting ground truth in MD simulation and Markov state modeling is rather unsettling. A pessimistic, but not unreasonable, view is that without a principled approach in defining true reaction coordinates for protein folding one has to make do with rather ad hoc and heuristic approaches in defining the discrete states with unexamined consequences, and hence this unsettling aspect may be with us for a long while and the truth may be elusive.

We agree with the reviewers that Markov state modeling depends on the choice of coordinate or metric one examines. It generally quite difficult to find the true reaction coordinate for defining states in complex systems. On the other hand, in the case of protein folding problem, theoretical and simulation studies over the decades have shown that the fraction of native contacts, Q, is the best reaction coordinate for describing folding process. For example, Best and coworkers, by analyzing folding MD simulation data of various proteins, showed that Q successfully captures the transition states for all the proteins (Best, Hummer and Eaton, 2013). The FRET efficiency \epsilon was chosen for matching with smFRET data and for differentiating the compact and extended states. These two coordinates (Q and \epsilon) are necessary and sufficient for describing the progress of folding as well as refining the bias of the current force field parameters toward compact structures. Also, since the structure of WW domain is rather simple (three-strands, betasheet-only), these coordinates are enough for uniquely identifying tertiary structures. Considering these reasons and the validation through the comparison with the Phi-value analysis, we think that the picture obtained by the current MSM may be the robust against the choice of coordinates. We have added this discussion in the Results section (subsection “Molecular dynamics simulations”, second paragraph).

Here is a suggestion in this respect: The authors employ Equation 5 for a model of FRET efficiency, however, more detailed models exist such as http://aip.scitation.org/doi/pdf/10.1063/1.3230974

The authors could evaluate their approach and determine if a more detailed model might be justified or needed in making contact with experimental measurements.

We thank the reviewers for letting us know the important paper. Following the reviewers’ suggestion, we applied our scheme with the new definition for FRE efficiency 𝜀 = 1⁄[Λ + (𝑟⁄𝑅_*)^-] where Λ > 1. Here, we optimized the likelihood function with Λ = 1.065, a value used in the paper. The result has been added as Figure 2—figure supplement 6. In the figure, although the transition states look more stabilized, the locations of stabilized states are qualitatively same as the original one (Λ = 1.0). This suggests that the story of the manuscript would be robust against the definition of ε (subsection “Markov state model and semi-supervised learning”, seventh paragraph).

2) While concern 1 is more general, there is another question more specific to connecting MD and the smFRET measurements. It will be necessary for the authors to show how well the HMM learned T(\tau) conform with MD simulation generated T(\tau) after the refined discrete states are obtained, in a statistically significant way. For example, the authors should compare the HMM rates and the MD rates for the subsets of state-state transitions with the highest and relevant lowest rates among the 87x87 parameters. In addition, the flux and pathway analysis should also be carried out using MD-derived rates after the refined states are obtained by the HMM model, for key paths carrying most of the flux, and compared with the analysis using HMM rates. The authors may have already done some of these in the transition state ensemble analysis, and may already have all necessary data.

This is important for validating the connection between the sm-FRET (HMM modeling) time series and the MD trajectories the authors are proposing. It would be fantastic if all works out, as this would naturally suggest that in the future one no longer needs to do MD simulation. Rather, one could use smFRET to fit dynamic rates, as long as the discrete states can be obtained, for example, by some other means. However, one should be cautious about being overly optimistic here, as there are significant problems in the enormous search space of very high dimension. It will be very interesting if the authors can show that the HMM rates and MD rates after refinement are the same.

In this study, the states were defined through the structural clustering of 400 microseconds MD simulation data. As 400 microseconds should be long enough for sampling all possible structures of WW domain, we think that the defined states are transferable during the HMM refinement. We would like to point that we already characterized the difference in T(\tau) from the MD and HMM refinement: Figure 4 shows the result of flux analysis of folding trajectories, which reveals the bottleneck of folding (the transition state ensemble). Figure 5 shows the results of the pathway analysis, which reveals folding pathways with largest fluxes.

As the reviewers pointed out, it would be nice if all transition probabilities are accurately estimated by the HMM from smFRET data. However, as shown by the HMM refinement starting from a random matrix (in response to comment 8), it is still difficult to estimate accurate kinetics along the Q axis from the limited number of smFRET observations. This suggests an importance of complementary use of both data (MD and smFRET), as proposed by our scheme (the MD information as the initial condition for HMM).

3) Potential over-fitting remains a concern, as there is a large number (87x87 = 7,569) parameters that need to be fitted. While the authors qualitatively compared results using each half of the data through plotting, it would be necessary to do some cross-validation tests, which are standard in machine learning. Specifically, the authors could use (n-1)/n balanced smFRET and MD data (e.g. similar number of folding trajectories) to identify refined states, and estimating \tau values, then test whether the remaining 1/n data falls within these defined states, and whether the unseen epsilon histograms can be reproduced accurately.

We thank the reviewers for the important suggestion. Following the reviewers’ comment, we have performed a k-fold cross validation test (with k=4) by partitioning the smFRET data set into k subsets (k-1 subsets as the training data, a single subset as the test data). We have evaluated the mean squared errors between two normalized FRET efficiency histograms calculated from the model and the subsets of smFRET data. The prediction error (the so-called cross validation error) for the test data was found to be 11.8 × 10⁷⁸, which is quite small. This suggests that the overfitting may not be a major issue in our modeling. We have added this result (Figure 3—figure supplement 1) in the Materials and methods section (subsection “Analysis of Markov state models”, second paragraph).

4) In Equation 5, the average values are used in calculating \epsilon. It seems that kappa has little difference at different regions of Q. Is the difference along the dimension of \epsilon that is used for clustering due to difference in distances r(s_i)? It will be helpful to also show both mean and variance of r(s_i) for each of the major discrete states, and justify whether the degree of in-cluster and between-cluster heterogeneity/homogeneity conforms to the Markovian assumption of each discrete state. This relates to the issue of validating the defined discrete states.

Following the reviewers’ comment, we have calculated the orientation factor k² along the e (FRET efficiency) axis. The average (standard deviation) in each region was k² = 0.66 (0.69) for e = 0.00-0.25, 0.64 (0.68) for e = 0.25-0.50, 0.61 (0.66) for e = 0.50-0.75, and 0.63 (0.64) for e = 0.75-1.00. Since the orientation factor k² show little dependence on both Q and e axes, it would be reasonable to use the isotropic average approximation as employed in the manuscript. We have added this result in the Materials and methods section (subsection “Markov state model and semi-supervised learning”, third paragraph).

Also, we have calculated both mean and variance of donor-acceptor distance r(s_i) for each state. It has been plotted in the 2D space spanned by Q and r (added as Figure 2—figure supplement 2). The figure shows that there are distance gaps between states, which is supportive for the Markovian assumption for the states. Large standard deviations are observed within states of small donor-acceptor distances. This is due to the worsening of spatial resolution in FRET efficiency when r(s_i) is small compared to R₀ as discussed in the information-binning study by Haw-Yang and coworkers (Watkins and Yang, 2004). We have added this result (Figure 2—figure supplement 2) in the Materials and methods section (subsection “Molecular dynamics simulations”, last paragraph).

5) The calculations of FRET efficiencies from MD trajectories need to be stated explicitly. Furthermore, the estimation of FRET efficiencies from experimentally measured photon sequences also needs to be presently clearly. The latter in fact is a long-standing issue in the single-molecule community, for example, the time binning convention and the information-base binning by Haw Yang.

We have explicitly described the calculations of FRET efficiencies from MD trajectories and photon sequences in the Results and Materials and methods sections. In particular, we have now used different terms for two FRET efficiencies (‘FRET efficiency \epsilon’ estimated from MD, ‘Measured FRET efficiency E’ from photon counts). We here used the conventional constant-time binning because the Markov state modeling is based on the constant-time binning. As a future direction, it is promising to apply a continuous-time Markov modeling (McGibbon and Pande, 2015), which may allow us to use the information-based binning (Watkins and Yang, 2004) or photon-by-photon analysis (Gopich and A. Szabo, 2009). We have described this direction in the Discussion section, adding citations for these works (fifth paragraph).

6) It is not obvious how the folding pathways were decomposed into individual pathways for both MSM models. How this procedure was conducted would critically impact the interpretation of results. This part thus needs careful and clear explanation.

We have described how the folding pathways were decomposed and its impact in the Results and Materials and methods sections. In the Results section we have added the following:

“The flux of folding trajectories was decomposed into individual pathways for both models (Figures 5A and 5B). […] The figures show dominant folding pathways with each of largest fluxes contributing to 50% of the total flux.”

7) Certain key references of interring dynamics from smFRET trajectories were missed in the context of this work, such as Haas et al., 2013.

We thank the reviewers for letting us know the important paper. We have added following references for recent advances in extracting more dynamical and structural information from smFRET trajectory (Introduction, first paragraph):

Haas, Yang and Chu, 2013; Matsunaga, Kidera, and Sugita, 2015; Sun, Morrell and Yang, 2016; Hoefling et al., 2011.

8) Obtaining a solution of similar likelihood value starting from a random matrix is a bit worrisome. Can the authors elaborate more on how this model was ruled out? Furthermore, the convergence of likelihood optimization should also be discussed to reveal the robustness of this approach.

We have extended the optimization calculation of a random matrix and examined the convergence of the log likelihood function. The result has been added as Figure 2—figure supplement 8. In the figure, the optimization of the log likelihood from random matrix looks stacked after 10000 iterations and its value is lower than that of MD data. This suggests that the optimization gets trapped in a local minimum in the parameter space. From these observations, we conclude that the model from a random matrix could be ruled out, or at least, global optimization of a random matrix is practically very inefficient. We have added this result (Figure 2—figure supplement 8) in the Materials and methods section (subsection “Markov state model and semi-supervised learning”, last paragraph).

[Editors' note: further revisions were requested prior to acceptance, as described below.]

The manuscript has been improved but there are some remaining issues that need to be addressed before acceptance, as outlined below:

Overview:

The unsettling aspects of shifting ground truth, or rather, the lack of ground truth remains. This is, however, an issue the community of protein folding simulations has to deal with.

Using Q as an adequate reaction coordinate is similarly problematic and inadequate, as Q is no substitute to reaction coordinates in the strict sense of earlier classic studies of Chandler and Dinner. But this again is a community problem, and we do not wish to specifically penalize the authors of the current work. Therefore, we accept the authors' arguments, although at least one of the reviewers had some fundamental concerns about the entire set of approaches in this area.

Specific points to address:

1) There is one issue that the authors did not provide answers to but they can easily do, namely, regarding the question of providing details of the 87x87 rates from the Hidden Markov model and the corresponding rates from MD simulation. It will be necessary to make first a general statement on how similar/different these two sets of matched/paired rates are, and second, provide details, e.g., in the Appendix, with side-by-side or overlaid histograms. If these two distributions are similar, it would be a comforting result. However, even if these two distributions do not agree, it will be important for readers to know this fact so they can directly draw their own inferences. This can be easily done with existing data.

Following the reviewers’ comment, we have added new plots comparing the rates of the two MSMs as Figure 2—figure supplement 4.

First (in panel A), we decomposed each transition probability matrix to obtain characteristic time scales before and after the hidden Markov modeling. We can see that the slowest time scale increases from 2.6 microseconds to 100 microseconds after the hidden Markov modeling. These time scales are related to folding/unfolding transitions of WW domain. Generally, simulation time scale is faster than experiments because of lower viscosity of TIP3P water model. This gap was improved by the information of the single-molecule experimental data.

In panel B, 87x87 T_ij rates before and after the hidden Markov modeling are directly compared in a scatter plot where each rate is colored using the FRET efficiency e_i. T_ij rates from states with high FRET efficiencies are well correlated with each other. This means that T_ij rates between compact states derived from MD simulations are consistent with the experimental data. In contrast, T_ij rates from states with middle or low FRET efficiencies are less correlated. This means that hidden Markov modeling mainly “refined” the rates between elongated states.

In this revision, we have added these explanations in the main text.

2) A grammar and expression overhaul is highly recommended to enhance the readability of the work.

In this revision, a native speaker has carefully edited the language. We have also checked a number of expressions.

Author details

1. Yasuhiro Matsunaga

   1. Computational Biophysics Research Team, RIKEN Center for Computational Science, Kobe, Japan
   2. JST PRESTO, Kawaguchi, Japan

   Contribution

   Conceptualization, Data curation, Software, Formal analysis, Funding acquisition, Validation, Investigation, Visualization, Methodology, Writing—original draft, Writing—review and editing

   Competing interests

   No competing interests declared

   "This ORCID iD identifies the author of this article:" 0000-0003-2872-3908

2. Yuji Sugita

   1. Computational Biophysics Research Team, RIKEN Center for Computational Science, Kobe, Japan
   2. Theoretical Molecular Science Laboratory, RIKEN Cluster for Pioneering Research, Wako, Japan
   3. Laboratory for Biomolecular Function Simulation, RIKEN Center for Biosystems Dynamics Research, Kobe, Japan

   Contribution

   Conceptualization, Data curation, Supervision, Funding acquisition, Writing—original draft, Project administration, Writing—review and editing

   For correspondence

   sugita@riken.jp

   Competing interests

   No competing interests declared

   "This ORCID iD identifies the author of this article:" 0000-0001-9738-9216

• 3,344

  Page views

• 538

  Downloads

• 30

  Citations

Article citation count generated by polling the highest count across the following sources: Crossref, Scopus, PubMed Central.