<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 MSM of the alanine dipeptide
In this example we use different methods for estimating the transition and rate matrices for a simple two-state MSM of the alanine dipeptide. Again, first one must download files from the following link.
Step1:  First we read the trajectory data, here corresponding to the Gromacs xtc files using the MDtraj library.
Step2:  Then we discretize the data using Ramachandran states. Only alpha (A), extended (E) and left-handed (L) configurations are accepted.
Step3:  In order to enforce detailed balance we revert the assigned time series and use both trajectories to build the Markov model.
Step4:  Using the discretized data we build the MSM invoking the SuperMSM class.
Step5:  We then calculate the transition matrix for a set of lag times.
Step6:  The rate matrix can be calculated using a simple Taylor expansion. The  method produces acceptable solutions for short lag times although the result rapidly diverges from the transition matrix relaxation time at long lag times. The rate matrix elements have the following values
Step7:  Alternatively, we can use the lifetime based method, where we use branching probabilities and average lifetimes. The rate is almost identical to that from the Taylor expansion.
Step8:  Clearly, both methods yield approximately the same result. Finally, we use Buchete and Hummer's maximum likelihood propagator-based method, which is also implemented in MasterMSM.
Step9:  Finally we compare the eigenvalues and equilibrium probabilities from the three methods and those from the transition probability matrix.
<ASSISTANT_TASK:>
Python Code:
%load_ext autoreload
%autoreload 2
%matplotlib inline
%pylab inline
import math
import numpy as np
import copy
import seaborn as sns
sns.set(style="ticks", color_codes=True, font_scale=1.5)
sns.set_style({"xtick.direction": "in", "ytick.direction": "in"})
import mdtraj as md
from mastermsm.trajectory import traj
Explanation: MSM of the alanine dipeptide
In this example we use different methods for estimating the transition and rate matrices for a simple two-state MSM of the alanine dipeptide. Again, first one must download files from the following link.
End of explanation
tr = traj.TimeSeries(top='data/alaTB.gro', traj=['data/protein_only.xtc'])
Explanation: First we read the trajectory data, here corresponding to the Gromacs xtc files using the MDtraj library.
End of explanation
tr.discretize(states=['A', 'E', 'L'])
tr.find_keys()
Explanation: Then we discretize the data using Ramachandran states. Only alpha (A), extended (E) and left-handed (L) configurations are accepted.
End of explanation
tr_rev = copy.deepcopy(tr)
tr_rev.distraj = [x for x in reversed(tr.distraj)]
Explanation: In order to enforce detailed balance we revert the assigned time series and use both trajectories to build the Markov model.
End of explanation
from mastermsm.msm import msm
msm_alaTB = msm.SuperMSM([tr, tr_rev], sym=False)
Explanation: Using the discretized data we build the MSM invoking the SuperMSM class.
End of explanation
lagts = [1, 2, 5, 7, 10, 20, 50, 100]
for lt in lagts:
    msm_alaTB.do_msm(lt)
    msm_alaTB.msms[lt].do_trans()
    #msm_alaTB.msms[i].boots(plot=False)
Explanation: We then calculate the transition matrix for a set of lag times.
End of explanation
msm_alaTB.msms[1].do_rate()
print (msm_alaTB.msms[1].rate)
rate_Taylor = msm_alaTB.msms[1].rate
tau_Taylor = msm_alaTB.msms[1].tauK
peq_Taylor = msm_alaTB.msms[1].peqK
Explanation: The rate matrix can be calculated using a simple Taylor expansion. The  method produces acceptable solutions for short lag times although the result rapidly diverges from the transition matrix relaxation time at long lag times. The rate matrix elements have the following values:
End of explanation
msm_alaTB.do_lbrate()
print (msm_alaTB.lbrate)
Explanation: Alternatively, we can use the lifetime based method, where we use branching probabilities and average lifetimes. The rate is almost identical to that from the Taylor expansion.
End of explanation
msm_alaTB.msms[1].do_rate(method='MLPB', init=msm_alaTB.lbrate, \
                          report=True)
print (msm_alaTB.msms[1].rate)
Explanation: Clearly, both methods yield approximately the same result. Finally, we use Buchete and Hummer's maximum likelihood propagator-based method, which is also implemented in MasterMSM.
End of explanation
fig, ax = subplots(1,2, figsize=(10,5))
ax[0].plot(msm_alaTB.msms[1].peqT, 'x-', markersize=16, \
           markeredgewidth=2, color='g', markeredgecolor='g', \
           markerfacecolor='None', label='T')
ax[0].plot(peq_Taylor, 'o-', markersize=16, markeredgewidth=3, \
           color='b', markeredgecolor='b', markerfacecolor='None',\
           label='Taylor')
ax[0].plot(msm_alaTB.peqK, 's-', markersize=16,markeredgewidth=2, \
           color='r', markeredgecolor='r', markerfacecolor='None', \
           label='LB')
ax[0].plot(msm_alaTB.msms[1].peqK, '^-', markersize=16, \
           markeredgewidth=2, color='c', markeredgecolor='c', \
           markerfacecolor='None', label='MLPB')
ax[0].set_xlim(-0.2,2.2)
ax[0].set_xlabel('states', fontsize=20)
ax[0].set_ylim(5e-3,1)
ax[0].set_xticks([0,1,2])
ax[0].set_xticklabels(["A","E", "L"])
ax[0].set_ylabel('P$_\mathrm{eq}$', fontsize=20)
ax[0].set_yscale('log')
ax[1].plot(msm_alaTB.msms[1].tauT, 'x-', markersize=16, \
           markeredgewidth=2, color='g', markeredgecolor='g', \
           markerfacecolor='None', label='T')
ax[1].plot(tau_Taylor, 'o-', markersize=16, markeredgewidth=3, \
           color='b', markeredgecolor='b', markerfacecolor='None', \
           label='Taylor')
ax[1].plot(msm_alaTB.tauK, 's-', markersize=16,markeredgewidth=2, \
           color='r', markeredgecolor='r', markerfacecolor='None', \
           label='LB')
ax[1].plot(msm_alaTB.msms[1].tauK, '^-', markersize=16, \
           markeredgewidth=2, color='c', markeredgecolor='c', \
           markerfacecolor='None', label='MLPB')
ax[1].set_ylabel(r'$\tau$ [ps]', fontsize=20)
ax[1].set_xlabel('eigenvalues', fontsize=20)
ax[1].set_xlim(-0.2,1.2) 
ax[1].set_ylim(0,100) 
ax[1].set_xticks([0,1])
ax[0].legend(loc=3)
ax[1].legend(loc=1)
plt.tight_layout()
Explanation: Finally we compare the eigenvalues and equilibrium probabilities from the three methods and those from the transition probability matrix.
End of explanation
<END_TASK>