<jupyter_start><jupyter_text>HIV-1 protease Markov State Model, Gating Factor Calculation and Conformational AnalysisAuthor: <NAME>Correspondence: <EMAIL>, Affiliation: 1. Heidelberg Institute for Theoretical Studies, HITS gGmbH 2. European Moelcular Biology LaboratoryThis notebook contains the computational workflow for analysing (MD) simulations and Markov state model analyses of apo HIV-1 protease conformational gating for the manuscript:S.Kashif. Sadiq‡, <NAME>, <NAME>, <NAME> (2021)   A multiscale approach for computing gated ligand binding from molecular dynamics and Brownian dynamics simulations  Import modules<jupyter_code>%pylab inline
import hivprgating
from hivprgating import *<jupyter_output><empty_output><jupyter_text>Markov State Model Analysis  Load coordinates and trajectory list<jupyter_code>#Specify the details of where the trajs are and the top file
indir = '../traj/'
topfile =  '../reference/protein.pdb'
#Specify list of trajs from a matrix of: batch sim frames
#Load data on number of frame for each trajectory
fname = '../traj/sims_array.dat'
sims_array = read_int_matrix(fname)
#Choose list of trajs from array: 
traj_list = traj_list_from_sims_array_xtc(sims_array, indir)

#Featurizer
feat = coor.featurizer(topfile)<jupyter_output><empty_output><jupyter_text>Load Pre-Computed Multidimensional Feature Data<jupyter_code>#Number of snapshot per trajectory
nsnaps=1300

# Get lambda, Alpha/beta and firemans grip coordinates as a pyEMMA data object
dir_array=['../features/lambda']

Y = multidir_obj(dir_array, sims_array,nsnaps)

# Convert to 2D matrix of all features and snapshots
Z=np.reshape(Y,(np.shape(Y)[0]*np.shape(Y)[1],np.shape(Y)[2]))

#Adds lambda modulus as a column
Z=np.transpose(np.vstack((np.transpose(Z),np.linalg.norm(Z[:,0:3],axis=1))))

# 3column lambda projection data in cylindrical coordinates
#Introduce 50 deg offset to zero point of angle
Z_cyl=xyz_to_cyl_coords(Z[:,:3],50)
#Adds cylindrical coords to Z as 3 columns at end
Z=np.hstack((Z,Z_cyl))

#Picks out the first three columns - later could be used for MSM
Yred=[Y[x][:,0:3] for x in range(np.shape(Y)[0])]

#Picks out the first three columns of Y but in cylindrical coordinates
Yred_cyl=[xyz_to_cyl_coords(Y[x][:,0:3],50) for x in range(np.shape(Y)[0])]

#Here we select a 461 subset of data up to 50 ns representing size of earlier data set
#Randomly pick 461 trajs from all trajs started from semi-open (500)
trajlist=np.random.permutation(np.where(sims_array[:,1]<=500)[0])[:461]
#trajlist=[x for x in range(461)]
t_end=500
#Z2 are the indices in Z marking the location of the subset data
Z2_inds=[j for i in [[trajid*nsnaps + i for i in range(t_end)] for trajid in trajlist] for j in i]

#Labelling of each column in Z matrix data
Zlabels=["$\lambda_{x}$ (\AA)","$\lambda_{y}$ (\AA)","$\lambda_{z}$ (\AA)",
        "$\lambda$ (\AA)","$\lambda_{r}$ (\AA)","$\lambda_{\Theta}$ ($^{\circ}$)","$\lambda_{z}$ (\AA)"]<jupyter_output><empty_output><jupyter_text>Cluster Feature Data into Microstates<jupyter_code>n_clusters = 100   # number of k-means clusters

#To compute clusters based on the 3D lambda metric in cylindircal polar coords - Note, this will compute a new clustering
#lambda_cl_obj = coor.cluster_kmeans(Yred_cyl, k=n_clusters, tolerance=1e-05, max_iter=100)

### Save and Load Previous Cluster object
#lambda_cl_obj.save('./MSM_obj/cl_obj100_lambda_cyl')
#To load the precomputed clustering used in the manuscript:
lambda_cl_obj=pyemma.load('../MSM_obj/cl_obj100_lambda_cyl')

#Create discretized trajectories - where each snapshot is assigned to belong to a cluster e.g. microstate
lambda_dtrajs = lambda_cl_obj.dtrajs<jupyter_output><empty_output><jupyter_text>Relaxation Timescale analysis<jupyter_code>#Compute relaxation timescales
#lambda_its = msm.timescales_msm(lambda_dtrajs, lags=500, nits=12, errors='bayes')

### Save and Load Previous Implied Timescales Object
#lambda_its.save('./MSM_obj/its_cl100_lambda_cyl')
##To load the precomputed relaxation timescales obtained for the manuscript
lambda_its=pyemma.load('../MSM_obj/its_cl100_lambda_cyl')
#Plot Timescale curves
plot_its(mplt,lambda_its,50,1000)

#Plot Ratio of Timescale curves    
plot_timescale_ratios(lambda_its)
#save_figure('its.png')

# Plot Kinetic Variance
plot_kinetic_variance(lambda_its)
#save_figure('kin_var.png')

#save_current_fig2(plt, '../figures/its.png')<jupyter_output><empty_output><jupyter_text>MSM Analysis at chosen lag time<jupyter_code>msm_lag = 300
#To compute MSM based on defined lag time
#M = msm.bayesian_markov_model(lambda_dtrajs, msm_lag)

### Save and Load Previous MSM
#M.save('./MSM_obj/M_cl100_tau30')
#To load precomputed MSM used in the manuscript
M=pyemma.load('../MSM_obj/M_cl100_tau30')
M_stat=M.stationary_distribution
PI_all = np.hstack([M_stat[dtraj] for dtraj in M.discrete_trajectories_full])

#np.set_printoptions(precision=3)
#print(PI_all)
#print(M_stat)<jupyter_output><empty_output><jupyter_text>PCCA Observation probability and Membership matrices<jupyter_code>n_sets = 5
M.pcca(n_sets)
pccaX = M.metastable_distributions
pccaM = M.metastable_memberships  # get PCCA memberships
pcca_sets = M.metastable_sets
pcca_assign = M.metastable_assignments<jupyter_output><empty_output><jupyter_text>Hidden Markov Model (HMM) approach - Baum-Welch Optimized Coarse-Graining<jupyter_code>#This takes a LONG TIME!
#hmm = M.coarse_grain(5)

# Save and Load previous HMM - This file is large ~ 2Gb 
#hmm.save('./MSM_obj/hmm_cl100_tau30_meta5')
hmm=pyemma.load('../MSM_obj/hmm_cl100_tau30_meta5')
hmm_dist = hmm.metastable_distributions
hmm_membership = hmm.metastable_memberships  # get PCCA memberships
# memberships over trajectory
hmm_dist_all = [np.hstack([hmm_dist[i,:][dtraj] for dtraj in hmm.discrete_trajectories_full]) for i in range(n_sets)]
hmm_mem_all = [np.hstack([hmm_membership[:,i][dtraj] for dtraj in hmm.discrete_trajectories_full]) for i in range(n_sets)]
hmm_sets = hmm.metastable_sets<jupyter_output><empty_output><jupyter_text>Plot projected observation probability distribution of microstates across collective variable landscape<jupyter_code>cl_obj=lambda_cl_obj
MSM_dims=np.array([4,5,6])
MSM_ticks=[[0,10,20,30,40], [0,100,200,300], [-10,0,10]]

dim=[5,4]
fig, axes = subplots(1, 5, figsize=(16, 3))
matplotlib.rcParams.update({'font.size': 12})
axes = axes.flatten()

np.seterr(invalid='warn') 
for k in range(n_sets):
        plot_sampled_density(k, Z[:,dim[0]], Z[:,dim[1]], hmm_dist_all[k],  dim, MSM_dims, MSM_ticks, Zlabels, 
                             ax=axes[k], cmap="Spectral", cbar=False)     
        ax=axes[k]
        mstate_color=['yellow','red','magenta','blue','black']
        ax.plot(cl_obj.clustercenters[hmm_sets[k],
np.where(MSM_dims==dim[0])[0][0]],cl_obj.clustercenters[hmm_sets[k],
np.where(MSM_dims==dim[1])[0][0]], linewidth=0, marker='o', markersize=2, color='black')
        
#fig.set_size_inches(15, 3)
#save_current_fig2(plt, '../figures/hmm_distributions_theta_r.png',pad=0.3)

plt.show()
dim=[5,6]
fig, axes = subplots(1, 5, figsize=(16, 3))
matplotlib.rcParams.update({'font.size': 12})
axes = axes.flatten()

np.seterr(invalid='warn') 
for k in range(n_sets):
        plot_sampled_density(k, Z[:,dim[0]], Z[:,dim[1]], hmm_dist_all[k],  dim, MSM_dims, MSM_ticks, Zlabels, 
                             ax=axes[k], cmap="Spectral", cbar=False)     
        ax=axes[k]
        mstate_color=['yellow','red','magenta','blue','black']
        ax.plot(cl_obj.clustercenters[hmm_sets[k],
np.where(MSM_dims==dim[0])[0][0]],cl_obj.clustercenters[hmm_sets[k],
np.where(MSM_dims==dim[1])[0][0]], linewidth=0, marker='o', markersize=2, color='black')
        
#save_current_fig2(plt, '../figures/hmm_distributions_theta_z.png',pad=0.1)<jupyter_output><empty_output><jupyter_text>Chapman-Kolmogorov test<jupyter_code>m=5
ck = M.cktest(m, memberships=hmm_membership, mlags=4, err_est=False)
#ck = M.cktest(5)
#matplotlib.rcParams.update({'font.size': 14})
plt.rc('text', usetex=False)
mplt.plot_cktest(ck, diag=True, figsize=(15,9), layout=(m,m), padding_top=0.1, y01=True, padding_between=0.3, dt=0.1, units='ns')


#save_current_fig2(plt, '../figures/cktest.png')
plt.show()<jupyter_output><empty_output><jupyter_text>Free Energy Plot Based on MSM Distribution with Microstate Cluster Centers<jupyter_code>plt.rc('text', usetex=True)
cl_obj=lambda_cl_obj
MSM_dims=np.array([4,5,6])
dim=[5,4]
mstate_color=['yellow','red','magenta','blue','black']
mpl.rcParams.update({'font.size': 12})
plt=plot_weighted_free_energy_landscape(Z,plt,dim[0],dim[1],Zlabels,cmap='Spectral',cbar=True,wg=PI_all,lev_max=0,shallow=True)
#Unweighted PMF
#plot_free_energy_landscape(Z,plt,dim[0],dim[1],Zlabels,cmap="Spectral")

#Plot All Cluster centers
#plot(cl_obj.clustercenters[:,np.where(MSM_dims==dim[0])[0][0]],cl_obj.clustercenters[:,np.where(MSM_dims==dim[1])[0][0]], linewidth=0, marker='o', markersize=8, color='k')

#Once HMM sets has been calculated - Calculate clusters belonging to each HMM metastable state
plot_metastable_sets(plt,cl_obj,hmm_sets,MSM_dims,dim,mstate_color)

#save_current_fig2(plt, '../figures/pmf_theta_r.png')
plt.rc('text', usetex=True)
cl_obj=lambda_cl_obj
MSM_dims=np.array([4,5,6])
dim=[5,6]
mstate_color=['yellow','red','magenta','blue','black']
mpl.rcParams.update({'font.size': 12})
plt=plot_weighted_free_energy_landscape(Z,plt,dim[0],dim[1],Zlabels,cmap='Spectral',cbar=True,wg=PI_all,lev_max=0,shallow=True)
#Unweighted PMF
#plot_free_energy_landscape(Z,plt,dim[0],dim[1],Zlabels,cmap="Spectral")

#Plot All Cluster centers
#plot(cl_obj.clustercenters[:,np.where(MSM_dims==dim[0])[0][0]],cl_obj.clustercenters[:,np.where(MSM_dims==dim[1])[0][0]], linewidth=0, marker='o', markersize=8, color='k')

#Once HMM sets has been calculated - Calculate clusters belonging to each HMM metastable state
plot_metastable_sets(plt,cl_obj,hmm_sets,MSM_dims,dim,mstate_color)

#save_current_fig2(plt, '../figures/pmf_theta_z.png')<jupyter_output><empty_output><jupyter_text>Macrostate thermodynamics, kinetics and gating factor calculation<jupyter_code>np.set_printoptions(precision=3,suppress=True)    

tau=30.
rho=hmm.pi
P=hmm.transition_matrix
kon=1000.*P/tau
init=[0,2]
fin=[1,3,4]
nd_kon=nondiag_rates(kon)
gam_kon= gamma_factor(kon,init,fin)
tau_c= tau_c(kon,init,fin)

#Macrostate thermodynamics 
#stationary distribution
print(rho)
#Free energy of macrostates
print(-0.596*np.log(rho))
#Zero-maximumed free energy of macrostates
print(-0.596*np.log(rho)-np.max(-0.596*np.log(rho)))
#Macrostate Transition kinetics matrix 
print(nd_kon)
#Conformation gating factor
print(gam_kon)
### Conformational Transition timescale
print(tau_c)<jupyter_output><empty_output><jupyter_text>Kinetic network of states<jupyter_code>#states 0 to 5 correspond to C1 - C5 respectively

pos=np.array([[-5, 1], [2, 3], [5, 0], [-1, 0], [0, 8]])    
#pos=np.array([[0, -5], [5, 0], [5, -5]])    
#pos=np.array([[-3, -1], [1, 3]])    
mplt.plot_markov_model(hmm, pos=pos, minflux=3e-4, arrow_labels=kon, arrow_label_format='%.3f')
gca().set_frame_on(False)<jupyter_output><empty_output><jupyter_text>Representative conformations  All Microstates Analysis<jupyter_code>#mac0=micro_order(Xsets,Xdist,0)

Xdist=hmm_dist
Xsets=hmm_sets

all_mac=np.empty((0,4))
for i in range(len(Xsets)):
    mac=micro_order(Xsets,Xdist,i)
    mac=np.hstack(( i*np.ones((np.shape(mac)[0],1)).astype(int), mac ))
    all_mac=np.vstack((all_mac,mac))


#Count frames of each microstate and append microstate stationary distribution and frame count to all_mac
all_mac_sets=all_mac[:,1].astype(int)
frame_count=np.empty((0,2))
for i in all_mac_sets:
    macro_state=all_mac[i,0].astype(int)
    micro_state=i
    micro_confs, all_snaps=micro_snap_extractor(macro_state, Xdist, [micro_state], lambda_dtrajs)
    frame_count=np.vstack((frame_count,micro_confs[0])).astype(int)

    
all_mac=np.hstack((all_mac, M.pi[all_mac[:,1].astype(int)].reshape(-1,1), frame_count[:,1].reshape(-1,1)))

# Table of macrostate, top micorstates, observation probability of the corresponding microstate
#print(all_mac[:,:3])<jupyter_output><empty_output><jupyter_text>Microstates plot<jupyter_code>#Select specific set of microstates that belong to a given metastable state
cl_obj=lambda_cl_obj
MSM_dims=np.array([4,5,6])
dim=[5,4]
mstate_color=['y','red','magenta','blue','black']
matplotlib.rcParams.update({'font.size': 12})
plt=plot_weighted_free_energy_landscape(Z,plt,dim[0],dim[1],Zlabels,fill=False,contour_label=False, 
                                        contour_color='0.70',wg=PI_all,lev_max=0,shallow=True,fsize=(16,16))

#Plot chi_conf_sets specific microstates
#cl_x=cl_obj.clustercenters[sets,np.where(MSM_dims==dim[0])[0][0]]
#cl_y=cl_obj.clustercenters[sets,np.where(MSM_dims==dim[1])[0][0]]
#plot(cl_x,cl_y, linewidth=0, marker='o', markersize=4, color='red')
#plt=annotate_microstates(plt,sets,cl_x,cl_y)

#plot_metastable_sets(plt,cl_obj,best_micro,MSM_dims,dim,mstate_color,msize=10,annotate=True, textsize=30)
#plot_metastable_sets(plt,cl_obj,top_sets,MSM_dims,dim,mstate_color,msize=8,annotate=True,textsize=16)
plot_metastable_sets(plt,cl_obj,hmm_sets,MSM_dims,dim,mstate_color,msize=4,annotate=True)
#plot_metastable_sets(plt,cl_obj,pcca_sets,MSM_dims,dim,mstate_color,msize=4,annotate=True)


#save_current_fig2(plt, '../figures/labelled_microstates_theta_r.png')<jupyter_output><empty_output><jupyter_text>Identify representative microstates and snapshots from X distributions<jupyter_code>energy_factor=1
top_sets=[]
top_sets_all=[]
top_b_set=[]
for i in range(len(Xsets)):
    chi_conf_sets, b_set = top_microstates(i, Xdist, Xsets, energy_factor)
    top_sets.append(chi_conf_sets[:5])
    top_sets_all.append(chi_conf_sets)
    top_b_set.append(b_set)
    
best_micro=[np.array(x[0]).reshape(1) for x in top_sets]    

#Number of microstates in each macrostate
print([len(top_sets_all[x]) for x in range(len(top_sets_all))])

print(top_sets_all)
#print(top_b_set)
print(best_micro)


bm1=[np.array(top_sets[0][:2]).reshape(2)]
bm2=[np.array(x[0]).reshape(1) for x in top_sets[1:]]
bm_states=bm1+bm2
print(bm_states)

print(top_sets)<jupyter_output><empty_output><jupyter_text>Best and Top microstates<jupyter_code>#Select specific set of microstates that belong to a given metastable state
cl_obj=lambda_cl_obj
MSM_dims=np.array([4,5,6])


mstate_color=['y','red','magenta','blue','black']
matplotlib.rcParams.update({'font.size': 12})

fig = plt.figure(figsize=(15, 6)) 
fig.subplots_adjust(hspace=0.4, wspace=0.4)
gs = gridspec.GridSpec(1, 2)

ax=plt.subplot(gs[0])
dim=[5,4]
plt=plot_weighted_free_energy_landscape(Z,plt,dim[0],dim[1],Zlabels,fill=False,contour_label=False, 
                                        contour_color='0.70',wg=PI_all,lev_max=0,shallow=True,standalone=False)
plot_metastable_sets(plt,cl_obj,best_micro,MSM_dims,dim,mstate_color,msize=10,annotate=True, textsize=25)
#plot_metastable_sets(plt,cl_obj,bm_states,MSM_dims,dim,mstate_color,msize=10,annotate=True, textsize=25)
#plot_metastable_sets(plt,cl_obj,top_sets_all,MSM_dims,dim,mstate_color,msize=8,annotate=True,textsize=16)

ax=plt.subplot(gs[1])
dim=[5,6]
plt=plot_weighted_free_energy_landscape(Z,plt,dim[0],dim[1],Zlabels,fill=False,contour_label=False, 
                                        contour_color='0.70',wg=PI_all,lev_max=0,shallow=True,standalone=False)

#plot_metastable_sets(plt,cl_obj,best_micro,MSM_dims,dim,mstate_color,msize=10,annotate=True,textsize=25)

plot_metastable_sets(plt,cl_obj,bm_states,MSM_dims,dim,mstate_color,msize=10,annotate=True, textsize=25)
#plot_metastable_sets(plt,cl_obj,top_sets_all,MSM_dims,dim,mstate_color,msize=8,annotate=True,textsize=16)

#save_current_fig2(plt, '../figures/best_microstates.png',pad=0.1)

plt.show()<jupyter_output><empty_output><jupyter_text>Extract and save trajectories of representative Macrostate conformations - from representative microstates and snapshots  Lambda sorted bestmicrostate snapshots<jupyter_code>#Trajectories are required for this step!

# Pre-extracted conformations and their associated data are already present in ../traj/macro_bestmicro_lamsort1000_xtcs - 
# so this step can be skipped.
# To redo this step - download all xtc trajectories from zenodo first and uncomment following commands within this cell

######################################
#inp = coor.source(traj_list, feat)
#systems=['c1','c1b','c2','c3','c4','c5']
#mac_mic=np.array([[0,0],[0,1],[1,0],[2,0],[3,0],[4,0]])
#macro_dir='../traj/macro_bestmicro_lamsort1000_xtcs'
#
#for i in range(len(systems)):
#    file_prefix=macro_dir+'/hmm_'+str(systems[i])
#    macro_state=mac_mic[i,0]
#    micro_state=top_sets[macro_state][mac_mic[i,1]]
#    mic_centroid=cl_obj.clustercenters[micro_state]
#    micro_confs, all_snaps=micro_snap_extractor(macro_state, Xdist, [micro_state], lambda_dtrajs)
#    
#    sort_all_snaps, sort_all_lam=lamsort_all_snaps(mic_centroid,all_snaps)
#    snaps, trajcounts, snaps_m=msm_conformation_selection(inp,
#                                                     sims_array,sort_all_snaps[:1000,:],sample=False,
#                                                      datfile=str(file_prefix)+'.dat',
#                                                      xtcfile=str(file_prefix)+'.xtc',
#                                                      mlabelfile=str(file_prefix)+'.mlabel.dat')<jupyter_output><empty_output><jupyter_text>Analyzing extracted/representative trajectories<jupyter_code># Analysis types:
# 0. Visual analysis of best 100 representative snapshots
# 1. lambda_x - comparison with previous paper
# 2. Curl analysis a-b analysis
# 3. RMSD w.r.t. semi-open and closed structures
# 4. distance-distance map
# 5. Volumetric density plots<jupyter_output><empty_output><jupyter_text>Extract Concatenated Feature data for representative snapshots of each macrostate<jupyter_code>#sys_id=0
#Zdata,snap_inds, snaps=conformation_data(sys_id, systems, macro_dir)        
#print(np.shape(Zdata))
#print(top_sets)
#print(snap_inds.reshape(-1,1))
#plot([(x+1)/10 for x in range(100)],Z[Zind,13],'r')


systems=['c1','c1b','c2','c3','c4','c5']
macro_dir='../traj/macro_bestmicro_lamsort1000_xtcs'

ZD=np.empty((1000,8,0))
SI=np.empty((1000,0))
S=np.empty((1000,2,0))

microstates=np.array([79,12,98,2,1,6])

for i in range(len(systems)):
    Zdata,snap_inds, snaps=conformation_data(Z, i, systems, macro_dir,sims_array)  
    
    #Proximity of cylindrical coords to centroid of a defined microstate
    #mic_centroid=cl_obj.clustercenters[microstates[i]]
    #dist=np.linalg.norm(Zdata[:,13:]-mic_centroid,axis=1).reshape(-1,1)

    #Proximity of lam cartesian coords to centroid of a defined microstate
    th_off=50
    mic_centroid=lambda_cl_obj.clustercenters[microstates[i]]
    mic_centroid_xyz=np.array([ mic_centroid[0]*np.cos(pi*(mic_centroid[1]+th_off)/180), 
                               mic_centroid[0]*np.sin(pi*(mic_centroid[1]+th_off)/180), mic_centroid[2] ])
    
    dist=np.linalg.norm(Zdata[:,4:]-mic_centroid,axis=1).reshape(-1,1)
    dist_xyz=np.linalg.norm(Zdata[:,:3]-mic_centroid_xyz,axis=1).reshape(-1,1)

    
    Zdata=np.hstack((Zdata, dist_xyz))
    
    ZD=np.dstack((ZD,Zdata))
    SI=np.hstack((SI,snap_inds.reshape(-1,1)))
    S=np.dstack((S,snaps))
    
#print(np.shape(ZD[:,:,0]))
#print(ZD[:,:,0])<jupyter_output><empty_output><jupyter_text>Visual Analysis<jupyter_code>fig = plt.figure(figsize=(18, 4)) 
fig.subplots_adjust(hspace=0.4, wspace=0.2)
gs = gridspec.GridSpec(2, 6)


gs_id=0
for i in range(len(systems)):
    ax=plt.subplot(gs[gs_id])
    ax.axis("off")
    img = mpimg.imread('../figures/conf_figs/png/'+str(systems[i])+'_top.png')
    imgplot = plt.imshow(img)
    
    plt.xticks([],fontsize=10, rotation=0)
    plt.yticks([],fontsize=10, rotation=0)
    
    gs_id+=1
    
for i in range(len(systems)):
    ax=plt.subplot(gs[gs_id])
    ax.axis("off")
    img = mpimg.imread('../figures/conf_figs/png/'+str(systems[i])+'_side.png')
    imgplot = plt.imshow(img)
    fig.set_facecolor("white")
    
    plt.xticks([],fontsize=10, rotation=0)
    plt.yticks([],fontsize=10, rotation=0)
    
    gs_id+=1

#save_figure('best_micro_extracted_snapshots.png')    
    
plt.show()<jupyter_output><empty_output><jupyter_text>Plot 1D Feature Data histogram distribution for extracted snapshots of each conformation  Lambda distribution<jupyter_code># Metric id follows Zlabels
#Zlabels=["$\lambda_{x}$ (\AA)","$\lambda_{y}$ (\AA)","$\lambda_{z}$ (\AA)",
#        "$\\alpha$ (\AA)","$\\beta$ (\AA)", "$\delta$ (\AA)",
#        "$\epsilon_{1}$ (\AA)","$\epsilon_{2}$ (\AA)", 
#        "$\kappa_{12}$ (\AA)", "$\kappa_{11}$ (\AA)", "$\kappa_{21}$ (\AA)", "$\kappa_{22}$ (\AA)",
#        "$\lambda$ (\AA)","$\lambda_{r}$ (\AA)","$\lambda_{\Theta}$ ($^{\circ}$)","$\lambda_{z}$ (\AA)"]

fig = plt.figure(figsize=(18, 2)) 
fig.subplots_adjust(hspace=0.6, wspace=0.4)
gs = gridspec.GridSpec(1, 6)


metric=3
metric_bounds=[0,40,10]
y_range=[0,1.2]
y_ticks=[0.5,1]

metric_range=metric_bounds[:2]
metric_ticks=[x for x in range(metric_bounds[0],metric_bounds[1],metric_bounds[2])]
D_space=np.array([y/10 for y in range(metric_bounds[0]*10,metric_bounds[1]*10,1)])
h=0.3
for i in range(len(systems)):
    ax=distribution_plot(gs[i],ZD,metric,i,D_space,h,'k') 
    ax=timeseries_axes(ax,metric_range,y_range,metric_ticks,y_ticks,None,None)

#print(D)
#ax1.hist(D, bins=D_plot,orientation='vertical',density=True,color='maroon',alpha=1,histtype='step')
#ax1.fill_betweenx(yrange,epan, fc='maroon',alpha=0.2)<jupyter_output><empty_output><jupyter_text>RMSD Analysis<jupyter_code>rmsd_so=np.empty((1000,5,0))
rmsd_cl=np.empty((1000,5,0))
ref_prmtop='../reference/protein.prmtop'
ref_pdb='../reference/protein.pdb'
cl_prmtop='../reference/p2nc/protein.prmtop'
cl_pdb='../reference/p2nc/protein.pdb'

for i in range(len(systems)):
    r_so=rmsd_analysis(i,systems, macro_dir, ref_prmtop, ref_pdb, ref_prmtop)
    r_cl=rmsd_analysis(i,systems, macro_dir, cl_prmtop, cl_pdb, ref_prmtop)
    
    rmsd_so=np.dstack((rmsd_so,r_so))
    rmsd_cl=np.dstack((rmsd_cl,r_cl))
    
    
#print(rmsd_cl[:,:,3])
fig = plt.figure(figsize=(18, 2)) 
fig.subplots_adjust(hspace=0.6, wspace=0.4)
gs = gridspec.GridSpec(1, 6)

D_space=np.array([y/10 for y in range(0,150,1)])
h=0.3
metric_range=[0,15]
metric_ticks=[x for x in range(0,15,5)]
y_range=[0,1.2]
y_ticks=[0.5,1]



for i in range(len(systems)):
    ax=distribution_plot(gs[i],rmsd_cl,4,i,D_space,h,'k',alph=0.5)
    ax=distribution_plot(gs[i],rmsd_so,4,i,D_space,h,'r',alph=0.5) 
    ax=timeseries_axes(ax,metric_range,y_range,metric_ticks,y_ticks,None,None)



#save_figure('best_micro_rmsd.png')

#Print mean and std of each conformation
print(np.mean(rmsd_so[:,:,4],axis=0),np.std(rmsd_so[:,:,4],axis=0))<jupyter_output><empty_output><jupyter_text>Analysis of Volmap calculated dx files for all macrostates<jupyter_code>#### Load grid files in a given vomap directory - for all macrostates

Dlist=[]
for i in range(len(systems)):
    dxfile="../volmap/"+systems[i]+".dx"
    D = density.Density(dxfile)
    Dlist.append(D)


#### Load grid files for reference structure - closed bound complex p2nc avg of 1000 snapshots in traj 1

Rlist=[]
dxpath="../volmap/"
dxfile=dxpath+"r.dx"
D = density.Density(dxfile)
Rlist.append(D)
dxfile=dxpath+"p.dx"
D = density.Density(dxfile)
Rlist.append(D)
dxfile=dxpath+"rcom.dx"
D = density.Density(dxfile)
Rlist.append(D)<jupyter_output><empty_output><jupyter_text>Volmap of closed flap HIV-1 protease with bound p2-NC substrate averaged over MD simulation<jupyter_code>D1=Rlist[0]
D2=Rlist[1]
#dim=0:x axis, y-z cross section - Lateral side
#dim=1:y axis, x-z cross section - Lateral face
#dim=2:z axis, x-y cross section - Top down
#dim=0
#extent=0
fig = plt.figure(figsize=(18, 6))
fig.subplots_adjust(hspace=0.2, wspace=0.3)
gs = gridspec.GridSpec(1, 3)
ax = plt.subplot(gs[0])
plt=plot_cross_section_overlay(plt, D1, D2, 2, 8,cmp1='binary',cmp2='Reds',overlay=True,flip=True,fsize=20)
ax = plt.subplot(gs[1])
plt=plot_cross_section_overlay(plt, D1, D2, 1, 0,cmp1='binary',cmp2='Reds',overlay=True,flip=True,fsize=20)
ax = plt.subplot(gs[2])
plt=plot_cross_section_overlay(plt, D1, D2, 0, 0,cmp1='binary',cmp2='Reds',overlay=True,flip=True,fsize=20)

#plt=plot_cross_section_overlay(plt, D1, D2, 0, 0,cmp1='binary',cmp2='Reds',overlay=False,flip=True)

fig.text(0.14,0.72,'A',fontsize=40,rotation=0)
fig.text(0.42,0.72,'B',fontsize=40,rotation=0)
fig.text(0.7,0.72,'C',fontsize=40,rotation=0)

fig.text(0.215,0.25,'z = 8 \AA',fontsize=20,rotation=0)
fig.text(0.495,0.25,'y = 0 \AA',fontsize=20,rotation=0)
fig.text(0.775,0.25,'x = 0 \AA',fontsize=20,rotation=0)

#plt.colorbar()

#Xtic,Xlab=make_ticks(D,0,20)
#Ytic,Ylab=make_ticks(D,1,20,-1)
#plt.xticks(Xtic,Xlab,fontsize=10, rotation=0)
#plt.yticks(Ytic,Ylab,fontsize=10, rotation=0)


#save_current_fig2(plt, '../figures/volplots_closed_bound_p2nc.png')


plt.show()<jupyter_output><empty_output><jupyter_text>Volmap of superposition of of p2-NC reference peptide with average of extracted macrostate conformations   Single cross-sectional plane<jupyter_code>fig = plt.figure(figsize=(24, 18)) 
fig.subplots_adjust(hspace=0.4, wspace=0.4)
gs = gridspec.GridSpec(5, 6)


dim=1
h_array=[[8,4,0,-4], [4,2,0,-2],[16,12,8,4]]
h_list=h_array[dim]

gs_id=0
for j in h_list:

    extent=j
    for i in range(len(systems)):

        ax = plt.subplot(gs[gs_id])        
        D1=Dlist[i]
        D2=Rlist[1]
        plt=plot_cross_section_overlay(plt, D1, D2, dim, extent,cmp1='binary',cmp2='Reds',overlay=True,fsize=20)

        gs_id+=1

plt.show()<jupyter_output><empty_output><jupyter_text>All three cross-sectional planes - xy,xz,yz<jupyter_code>systems=['c1','c1b','c2','c3','c4','c5']
#macro_dir='macro_bestmicro_lamsort100_xtcs'

fig = plt.figure(figsize=(24, 18)) 
fig.subplots_adjust(hspace=0.4, wspace=0.4)
gs = gridspec.GridSpec(5, 6)


#h_array=[[8,4,0,-4], [4,2,0,-2],[16,12,8,4]]
#h_list=h_array[dim]

fs=20

gs_id=0

#First line - Top view z=8
dim=2
extent=8
plt, gs_id = plot_cross_section_systems(plt, systems, gs, gs_id, Dlist, Rlist[1], dim, extent,fs)

#Second Line - Top view z=12
dim=2
extent=12
plt, gs_id = plot_cross_section_systems(plt, systems, gs, gs_id, Dlist, Rlist[1], dim, extent,fs)

#Third line - Lateral face view y=0
dim=1
extent=0
plt, gs_id = plot_cross_section_systems(plt, systems, gs, gs_id, Dlist, Rlist[1], dim, extent,fs)

#Fourth line - Lateral side view x=0
dim=0
extent=0
plt, gs_id = plot_cross_section_systems(plt, systems, gs, gs_id, Dlist, Rlist[1], dim, extent,fs)


#save_current_fig2(plt, '../figures/volplots.png')

plt.show()<jupyter_output><empty_output><jupyter_text>Structural Characterization - Compiled Master Figure<jupyter_code>fig = plt.figure(figsize=(24, 30)) 
fig.subplots_adjust(hspace=0.4, wspace=0.4)
gs = gridspec.GridSpec(8, 6)

gs_id=0

# Line 1
orient="top"
plt,gs_id = plot_structures_systems(plt, systems, gs, gs_id,orient)

#Line 2
orient="side"
plt,gs_id = plot_structures_systems(plt, systems, gs, gs_id,orient)

#Line 3 - RMSD w.r.t. semi-open and closed
D_space=np.array([y/10 for y in range(0,150,1)])
h=0.3
metric_range=[0,15]
metric_ticks=[x for x in range(0,15,5)]
y_range=[0,1.2]
y_ticks=[0.5,1]

for i in range(len(systems)):
    ax=distribution_plot(gs[gs_id],rmsd_cl,4,i,D_space,h,'m',alph=0.5)
    ax=distribution_plot(gs[gs_id],rmsd_so,4,i,D_space,h,'b',alph=0.8) 
    ax=timeseries_axes(ax,metric_range,y_range,metric_ticks,y_ticks,r'$d_{RMS}$ (\AA)',r'$\rho(d)$')
    
    gs_id+=1

    
#Line 4 - lambda

metric=3
metric_bounds=[0,40,10]
y_range=[0,1.2]
y_ticks=[0.5,1]

metric_range=metric_bounds[:2]
metric_ticks=[x for x in range(metric_bounds[0],metric_bounds[1],metric_bounds[2])]
D_space=np.array([y/10 for y in range(metric_bounds[0]*10,metric_bounds[1]*10,1)])
h=0.3
for i in range(len(systems)):
    ax=distribution_plot(gs[gs_id],ZD,metric,i,D_space,h,'k',alph=0.5) 
    ax=timeseries_axes(ax,metric_range,y_range,metric_ticks,y_ticks,r'$\lambda $ (\AA)',r'$\rho(\lambda)$')

    gs_id+=1
    
    
fs=20

#Line 5 - Top view z=8
dim=2
extent=8
plt, gs_id = plot_cross_section_systems(plt, systems, gs, gs_id, Dlist, Rlist[1], dim, extent,fs)

fig.text(0.91,0.44,'z = 8 \AA',fontsize=20,rotation=90)

#Line 6 - Top view z=12
dim=2
extent=12
plt, gs_id = plot_cross_section_systems(plt, systems, gs, gs_id, Dlist, Rlist[1], dim, extent,fs)

fig.text(0.91,0.34,'z = 12 \AA',fontsize=20,rotation=90)

#Line 7 - Lateral face view y=0
dim=1
extent=0
plt, gs_id = plot_cross_section_systems(plt, systems, gs, gs_id, Dlist, Rlist[1], dim, extent,fs)

fig.text(0.91,0.24,'y = 0 \AA',fontsize=20,rotation=90)


#Line 8 - Lateral side view x=0
dim=0
extent=0
plt, gs_id = plot_cross_section_systems(plt, systems, gs, gs_id, Dlist, Rlist[1], dim, extent,fs)

fig.text(0.91,0.14,'x = 0 \AA',fontsize=20,rotation=90)

fig.text(0.15,0.90,'C1a',fontsize=40)
fig.text(0.30,0.90,'C1b',fontsize=40)
fig.text(0.43,0.90,'C2',fontsize=40)
fig.text(0.56,0.90,'C3',fontsize=40)
fig.text(0.70,0.90,'C4',fontsize=40)
fig.text(0.84,0.90,'C5',fontsize=40)

fig.text(0.05,0.85,'A',fontsize=50)
fig.text(0.05,0.75,'B',fontsize=50)
fig.text(0.05,0.65,'C',fontsize=50)
fig.text(0.05,0.55,'D',fontsize=50)
fig.text(0.05,0.45,'E',fontsize=50)
fig.text(0.05,0.35,'F',fontsize=50)
fig.text(0.05,0.25,'G',fontsize=50)
fig.text(0.05,0.15,'H',fontsize=50)


#save_current_fig2(plt,'../figures/structural_characterization.png')

plt.show()<jupyter_output><empty_output><jupyter_text>Residue-Residue Interaction Analysis  Flap-Opposite Monomer Sidechain-Sidechain Analysis for a given conformation<jupyter_code>#Macrostate id: 0-5 corresponds to C1 (C1b)-C5 respectively
mac_id=0
ref_prmtop='../reference/protein.prmtop'
u=load_universe(mac_id,systems,macro_dir, ref_prmtop)
reslim_a=[43, 58]
reslim_b=[100, 198]
sidechain="not name H* and not name N CA O OX*"
min_c=general_min_dist_matrix(u, u,reslim_a,reslim_b,sidechain,
                              mfa=True,ax_points1=HIVPR_ax_points(u),ax_points2=HIVPR_ax_points(u))
resid_array=close_interactions(min_c,reslim_a,reslim_b,10)
C_mind,C_stats=res_res_array_min_sidechain_distance_trajectory(u, u, resid_array,sidechain)
#Reconvert monB
C_stats[:,1]-=99
#print(C_stats[C_stats[:,2]<=5,:])

reslim_a=[1, 99]
reslim_b=[142, 157]
sidechain="not name H* and not name N CA O OX*"
min_c_rev=general_min_dist_matrix(u, u,reslim_a,reslim_b,sidechain,
                              mfa=True,ax_points1=HIVPR_ax_points(u),ax_points2=HIVPR_ax_points(u))
resid_array=close_interactions(min_c_rev,reslim_a,reslim_b,10)
C_mind_rev,C_stats_rev=res_res_array_min_sidechain_distance_trajectory(u, u, resid_array,sidechain)
#Reconvert monB
C_stats_rev[:,1]-=99
#print(C_stats[C_stats[:,2]<=5,:])

C_stats_all=np.unique(np.vstack((C_stats,C_stats_rev)),axis=0)
#print(np.unique(C_stats_all,axis=0))

print(C_stats_all[C_stats_all[:,2]<=5.5,:])<jupyter_output><empty_output><jupyter_text>Flap-Opposite Monomer Hydrogen bond analysis for a given conformation<jupyter_code>#Macrostate id: 0-5 corresponds to C1 (C1b)-C5 respectively
mac_id=0
ref_prmtop='../reference/protein.prmtop'
u=load_universe(mac_id,systems,macro_dir,ref_prmtop)
hbonds = HBA(universe=u)
hbonds.d_a_cutoff=3.5
hbonds.d_h_a_angle_cutoff=150
hbonds.donors_sel='name N'
hbonds.hydrogens_sel='name H'
hbonds.acceptors_sel='name O'
#hbonds.donors_sel = hbonds.guess_donors("protein and resid 43:58")
#hbonds.hydrogens_sel = hbonds.guess_hydrogens("protein and resid 43:58")
#hbonds.acceptors_sel = hbonds.guess_acceptors("protein and resid 100:198")

hbonds.run()
hbdat=hbonds.hbonds
atmlist=hbond_atom_definitions(u,hbdat)
hbdat_atms=join_atom_defs_to_hbond_data(hbdat,atmlist)

#Flap A with Monomer B
reslim1=[43, 58]
reslim2=[100, 198]
unique_hb_list,hb_combo=unique_subset_hbond_data(hbdat_atms,reslim1,reslim2)
#print(unique_hb_list)

#Flap B with Monomer A
reslim1=[1, 99]
reslim2=[142, 157]
unique_hb_list_rev,hb_combo_rev=unique_subset_hbond_data(hbdat_atms,reslim1,reslim2)
#print(unique_hb_list_rev)

#Combine lists
un_hb_list_all=np.unique(np.vstack((unique_hb_list,unique_hb_list_rev)).astype('U21'),axis=0)
print(un_hb_list_all)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15993532941578834, 'fraction_numerical': 0.03954854957559929, 'mean_word_length': 6.669051025274201, 'pattern_counts': {'":': 0, '<': 130, '<?xml version=': 0, '>': 126, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 13, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '30', 'dir': 'jupyter-structured-clean-dedup', 'id': '274257', 'license': 'MIT', 'n_tokens_mistral': 13897, 'n_tokens_neox': 13014, 'n_words': 2036, 'path': 'data/MD/notebook/hivpr-conf-gating.ipynb', 'repo_name': 'kashifsadiq/hiv1pr-msm', 'size': '32219'}