Markov State Models: To Optimize or Not to Optimize

Markov state models (MSM) are a popular statistical method for analyzing the conformational dynamics of proteins including protein folding. With all statistical and machine learning (ML) models, choices must be made about the modeling pipeline that cannot be directly learned from the data. These choices, or hyperparameters, are often evaluated by expert judgment or, in the case of MSMs, by maximizing variational scores such as the VAMP-2 score. Modern ML and statistical pipelines often use automatic hyperparameter selection techniques ranging from the simple, choosing the best score from a random selection of hyperparameters, to the complex, optimization via, e.g., Bayesian optimization. In this work, we ask whether it is possible to automatically select MSM models this way by estimating and analyzing over 16,000,000 observations from over 280,000 estimated MSMs. We find that differences in hyperparameters can change the physical interpretation of the optimization objective, making automatic selection difficult. In addition, we find that enforcing conditions of equilibrium in the VAMP scores can result in inconsistent model selection. However, other parameters that specify the VAMP-2 score (lag time and number of relaxation processes scored) have only a negligible influence on model selection. We suggest that model observables and variational scores should be only a guide to model selection and that a full investigation of the MSM properties should be undertaken when selecting hyperparameters.


Model summaries
Table 1: Model summaries.A summary of the models discussed in the main text.The logistic transform of contact distances is specified below in terms of the center, c, and steepness, s, as: 'logit(dist.)'with '(c, s)' underneath, in units of A and A −1

Protein
No.

Chignolin model 1
This model has the largest median t 2 after random sampling optimisation.

BBA model 2
This model has the second largest median t 2 after random sampling.

BBA model 3
This model has the largest median t 2 using the distance feature after random sampling.

BBA model 4
This model has the largest median t 2 using the dihedral feature after random sampling.

BBA model 6
This model has the largest median t 2 after Bayesian optimisation of VAMP2 eq (2).

Figure 1 :
Figure1: Logistic transforms in select models.The horizontal axis is the contact distance in A, the vertical axis shows the logistic transform of that distance for BBA model 1 (blue), model 4 (orange) and model 5 (green).

Figure 2 :
Figure 2: Chignolin, model 1 timescales.Text inset shows the MSM hyperparameters; panel (a) shows the implied timescales for the first nine slow relaxation processes (τ = −1 ns); panel (b)shows the gap between successive successive timescales: the gap for process i is defined as t i /t i+1 (τ = −1 ns); panel (c) shows the mean first passage time between the the unfolded and folded state as a function of τ .

Figure 3 :
Figure 3: Chignolin, model 1 free-energy surface.Each panel shows different quantities projected onto the first two TICA components (IC1, IC2).Panel (a) shows the free energy surface; panel (b) shows the PCCA+ clustering into folded, unfolded and intermediate states with the crystal structure (PDB accession code 5AWL) marked with a star; panel (c) shows the 2nd right eigenvector (which corresponds to the slowest relaxation process) with an ensemble of structures corresponding to the extremes values of the eigenvector ('min', 'max' marked with a triangle and cross respectively).

Figure 4 :
Figure4: Chignolin, model 1 validation.The implied timescales plotted as a function of the lag-time.Blue, red, and green lines correspond to the slowest, second slowest, and third slowest process respectively.The solid lines correspond to the timescales from the maximum likelihood MSM.The dashed lines correspond to the mean of Bayesian MSMs.The coloured regions refer to the 0.95 confidence interval.The shaded region shows when the Markov lag time becomes equal to or longer than the implied timescale.

Figure 14 :
Figure 14: BBA, model 1 timescales.Text inset shows the MSM hyperparameters; panel (a) shows the implied timescales for the first nine slow relaxation processes (τ = 41 ns); panel (b) shows the gap between successive successive timescales: the gap for process i is defined as t i /t i+1 (τ = 41 ns); panel (c) shows the mean first passage time between the the unfolded and folded state as a function of τ .

Figure 15 :Figure 16 :
Figure 15: BBA, model 1 free energy surface.Each panel shows different quantities projected onto the first two TICA components (IC1, IC2).Panel (a) shows the free energy surface; panel (b) shows the PCCA+ clustering into folded, unfolded and intermediate states with the crystal structure (PDB accession code 1FME) marked with a star; panel (c) shows the 2nd right eigenvector (which corresponds to the slowest relaxation process) with an ensemble of structures corresponding to the extremes values of the eigenvector ('min', 'max' marked with a triangle and cross respectively).

Figure 32 :
Figure 32: Optimization of MSMs of Chignolin and BBA.The vertical axis is the optimization objective, the horizontal axis is the trial number.The thin line refers to the incumbent over the initialization data ('init.'),and the squares are the incumbent at each trial number.Panels (a) and (c) refer to Chignolin, (b) and (d) to BBA.Panels (a) and (b) refer to Bayesian optimization with either t 2 (red), or dualobjective optimization of both t 2 and the timescale gap, t 2 /t 3 .Panels (c) and (d) refer to optimization of VAMP2 eq (2) (red) and dual-objective optimization of both VAMP2 eq (2) and VAMP2 eq (2)/VAMP2 eq (3).The optimized values of t 2 are shown as labels.

Figure 33 :
Figure 33: Optimisation of the timescale and VAMP2 eq gap.The blue dots show the incumbent gap during optimisation (shown only for the multi-objective optimisation).Panels (a) and (b) refer to the timescale gap t 2 /t 3 while panels (c) and (d) refer to the VAMP2 eq (2)/VAMP2 e q(3).Panels (a) and (c) are for Chignolin optimisation while the panels (b) and (d) refer to BBA.

Figure 34 :Figure 35 :Figure 36 :Figure 37 :
Figure34: Pair plot of VAMP2 eq (k = 2) rank with different lag times.The panel at position (0, 1) plots the rank according toVAMP2 eq (2) against the rank according toVAMP2 eq (3), and similarly for other positions.This This model has the largest median t 2 after random sampling.