Based on: Optimal Dimensionality Reduction of Multistate Kinetic and Markov-State Models

(Hummer and Szabo, 2014, JAPC B) http://pubs.acs.org/doi/pdfplus/10.1021/jp508375q



In [2]:

    
import numpy as np
import numpy.random as npr
import pylab as pl
%matplotlib inline

Abstract

Goal: model reduction: optimally describe the dynamics of aggregated supersates formed by combining microstates

Tactic: match the time-dependent occupancy-number correlation functions of the superstates in full and aggregated systems

Alternate tactic: projection operator formalism

Introduction

Consider a high-dimensional dynamical system partitioned into two regions ($1$ and $2$), described by a two-state kinetic system (state space $\Omega=\{ 1,2 \}$, transition matrix $T$ contains transition rate to $2$ from $1$ $\equiv$ $R_{21}$, and reverse $\equiv$ $R_{12}$)
How should we choose the rate constants $T$?

Toy example: biased random walk on a bimodal energy landscape



In [22]:

    
def bimodal_dist(x,dist=1.8):
    comp1 = np.exp((-sum((x+dist*0.5)**2)))
    comp2 = np.exp((-sum((x-dist*0.5)**2)))
    return comp1 + comp2

def metropolis_hastings(density,n_samples=100,step_size=0.1,ndim=2):
    """ generates samples given a density function which maps from R^ndim
    to R using the defined step_size.
    returns (samples, acceptance ratio)"""
    x0 = npr.randn(ndim)
    d0 = density(x0)
    samples = np.zeros((n_samples,ndim))
    acceptances = 0
    for i in range(n_samples):
        x1 = x0 + npr.randn(ndim)*step_size
        d1 = density(x1)
        if npr.rand() <= (d1/d0):
            samples[i] = x1
            x0 = x1
            d0 = d1
            acceptances+=1
        else:
            samples[i] = x0
    
    return samples,1.0*acceptances/n_samples



In [23]:

    
step_size = np.sqrt(0.1)
samples, ar = metropolis_hastings(bimodal_dist,100000,step_size)
print("Acceptance ratio: {0:.4}".format(ar))
hist = pl.hist2d(samples[:,0],samples[:,1],bins=100);









    



Acceptance ratio: 0.8046



In [34]:

    
def state_decomp(x):
    return x[:,1] - x[:,0] < 0



In [36]:

    
discrete_states = np.array(state_decomp(samples),dtype=int)
1.0*sum(discrete_states)/len(samples)









    Out[36]:





0.49453



In [38]:

    
pl.plot(discrete_states[20000:20000+100])









    Out[38]:





[<matplotlib.lines.Line2D at 0x108c9a5d0>]



In [105]:

    
def transition_count_matrix(sequence):
    n = len(set(sequence)) # number of states
    C = np.zeros((n,n),dtype=int)
    for t in range(1,len(sequence)):
        i,j = sequence[t-1:t+1]
        C[i,j] += 1
    return C

def row_normalize(matrix):
    for i in xrange(len(matrix)):
        matrix[i] /= sum(matrix[i])
    return matrix

def trans_mat(sequence):
    return row_normalize(transition_count_matrix(sequence))



In [108]:

    
discrete_states









    Out[108]:





array([1, 1, 1, ..., 0, 0, 0])



In [109]:

    
T = trans_mat(discrete_states)



In [99]:

    
T[0] = T[0] / sum(T[0])
T[1] = T[1] / sum(T[1])



In [93]:

    
np.sum(T,axis=1)









    Out[93]:





array([ 1.,  1.])



In [94]:

    
T = T / np.sum(T,axis=1)



In [102]:

    
counts









    Out[102]:





array([[ 44713.,   5833.],
       [  5834.,  43619.]])



In [ ]:

    
np.sum(T,axis=1)



In [68]:

    
sum(T[0])









    Out[68]:





1.0025505433110127



In [58]:

    
sum(T)









    Out[58]:





array([ 1.,  1.])



In [59]:

    
T[0]









    Out[59]:





array([ 0.88458267,  0.11795276])



In [60]:

    
T[:,0]









    Out[60]:





array([ 0.88458267,  0.11541733])



In [ ]: