To answer our research questions we estimated a large number of Markov state models with different hyperparameters, measured their observables and analysed the results. The workflow may be summarized as follows. We used existing molecular dynamics trajectories of Chignolin and BBA and fit MSMs with 140 randomly sampled hyperparameters (hyperparameter trials) and recorded implied timescales, eigenvalues, and VAMP-2 scores for a range of different lag times (τ ). Each observable was estimated with confidence intervals using bootstrapping. This data constituted our hyperparameter trial data set and was analyzed in the first, second and last results subsections. A 'toy' three-state MSM model was constructed to highlight with the VAMP-2 score, for reversible transition matrix estimation. We then performed Bayesian optimization with a TPE surrogate function, with a variety of different objective functions, and using the hyperparameter trial data set to initialize the surrogate function. These results are discussed in the third part of the results.