Tram convergence check

[bonimba87]

Hi all, I have been running Tram on a set of biased and unbiased simulation data, and I argue the results are not converging but I cannot figure out why. I have a set of 24 umbrella sampling simulations (generated with a 1d umbrella coordinate, each with 9000 frames) and a set of 20 regular MD trajectories (each with 2000 frames). I ran TICA on the us data to extract independent components to use (together with the biasing coordinate) to generate the clusters. I then use the following syntax:

lags = [1, 2, 3, 5, 7, 10]

tram = pyemma.thermo.estimate_umbrella_sampling(
    us_trajs, us_dtrajs, us_centers = umbrella_centers, us_force_constants = spring_constants, 
md_trajs = w_md_trajs, md_dtrajs = w_md_dtrajs,  maxiter=10000000, maxerr=max_err, 
save_convergence_info=10000, estimator='tram',  lag=lags, init='mbar', init_maxiter=10000, 
init_maxerr=1.0E-10, direct_space=True, connectivity= 'post_hoc_RE', 
N_dtram_accelerations=N_dtram_accelerations)

I have also run simulations using the default option for 'connectivity', and alternating 'N_dtram_accelerations = 5 or 10' to speed up convergence (N_dtram_accelerations = 0 I also used, but Tram would not converge within the maximum wall time allowed). The 'max_err = 5E-15' threshold is mostly reached, so I assume the calculation has converged.However, if I plot the implied timescales, I get a puzzling behavior (see plot).

I was wondering if there is a more rigorous way to look into Tram convergence, and if perhaps some combination of the optional Tram parameters is recommended over some others. I would like something to improve the connectivity of the dataset (if possible), without corrupting the result quality.

Thank you for the help, Lorenzo

[fabian-paul]

Hi Lorenzo,

to check the connectivity of the data, a starting point would be to look at the connectivity information that the TRAM object provides. They are:

• tram.csets[k] # for every thermodynamic state number k: the states in k that contribute to the stationary distribution
• tram.active_set # tram.csets[k] projected along the k axis = states that contribution to the consensus stationary distribution
• tram.model_active_set[k] # for every thermodynamic state number k: the states that contribute kinetic information

Are the end states of the transition that you want to study included in the connected set? How does the connected sets change with lag time?

As usual with MSMs connectivity information is only approximate. The connectivity between different MSM states (within the same thermodynamic state) is an approximate measure because of recrossing. Consider a trajectory where you only have one real transition from state i to state j and the system remains in j for the rest of the trajectory. If there would be no projection error this fact could be detected simply by looking at the count matrix. However because of projection error, many spurious transitions into j and out of j will appear in the count matrix, so this problem can't be detected automatically. This problem usually shows up in the implied timescale plot because the trapped trajectory gets assigned a slow process of undefined time scale. I recommend to you to carefully check the connectivity in state space first. Is it possible to build an MSM with the unbiased data? If so, I would recommand building a simple MSM, because the problems are a bit more accessible there.

Then there is the possibility that overlap between thermodynamic states is lacking. This means that the factors used to reweight between ensembles are so small that they are effective zero. In principle only part of the conformational space has to overlap between ensembles to obtain a fully connected set. So it could be enough if the connectivity between ensembles was established by a single MSM state. However if some state transitions in one ensemble are not sampled reversibly or even missing, TRAM relies on a stronger overlap between ensembles. In the best case there is overlap between ensembles for every MSM state. One could start by checking what are the irreversibly connected states in the unbiased ensemble and then check whether these are reversibly connected in the biased simuiations. Again: state transitions are affected by recrossing, so some manual checks may be needed.

Best, Fabian

[bonimba87]

Hi Fabian, thank you for your kind reply.

I did build a Markov Model out of the unbiased simulation data. It is defined on basically the same subset of clusters which were visited by regular MD, as the tram.msm.transition_matrix is. The other clusters apparently drop out even if they were visited by umbrella simulations. I am not really sure whether I can consider it reaching convergence at some point, but for sure its estimates for the timescales are way larger than anything I have been able to get with Tram so far (see previous plot).

I also checked connectivity by building a matrix where the (i, j) entry is the number of clusters that MSMs on thermodynamics state i-th and j-th share (24 umbrellas, the 25-th column/row refers to the unbiased data). It seems the central area of the matrix is really poorly connected, could that explain the huge difference between the Tram and the unbiased MSM timescales?

MSM_convergence.pdf overlap_states.pdf

Again, thank you very much, Lorenzo

[fabian-paul]

Hi Lorenzo,

I think these implied time scales are not converged at all. (Even though some publications seem to sell plots like this as converged. But I think there is quite a bit of nonsense in the published literature.) They look like shifted version of the black identity line which means that the timescales are linearly diverging. From my experience this can be a consequence of missing connectivity between the Markov states. Typically some trajectory is trapped in some conformation (and never leaves the conformation). In MSMs this can give raise to a large implied time scale that appears to increase with lag time. There is some amount of manual preprocessing that has to be done in that case. You will have to detect the trapped state and handle it appropriately, e.g. exclude it from the analysis / seed new simulations in the state / build a series of umbrellas that pull the system out of the state / ...

Here is the recipe that I use to detect traps:

• build a normal MSM at some arbitrary lag time where you already see divergence of the ITS
• do PCCA (with n PCCA states where n=1+number of large ITS)
• plot the time series of the PCCA state, i.e.

set_number = 0 # let's look only at metastable set #0 for the moment
macro = np.concatenate((pcca.metastable_assignment, [-1])) # to handle states outside of the cset
for i in range(len(msm.discrete_trajectories_active)):
    if np.any(macro[msm.discrete_trajectories_active[i]]==set_number):
        plt.figure()
        plt.plot(macro[msm.discrete_trajectories_active[i]])
        plt.title(source.filenames[i])

If there is a trap, you will see one (or a few) time series that stay in the metastable set to the end of the trajectory + some recrossing events where some other trajectory enters the metastable set and immediately leaves it. These recrossing events are what renders the automated detection of trapped states difficult.

This could allow you to answer some basic questions about your data:

• Are there trapped states?
• If there are trapped states, do your trajectories capture the event of the system entering the trapped state or not (the trajectories already start in the trapped states)?
• If you trajectories enter, can you compensate for the lack of exit events by using the umbrella sampling data? Is the macro-state visited in the US simulations?

In you answer you say: "The other clusters apparently drop out even if they were visited by umbrella simulations." Does this refer to TRAM? The thermodynamic connectivity matrix that you show certainly has a gap in the middle. I would try to improve the sampling.

Best, Fabian