# Markov State Modeling

Last updated: May 29th, 2020

## Markov State Modeling Practice¶

#### Following parameters in table are corresponding to the parameters listed in Fig 1 (Markov Map): Ignore "distribution"¶

Based on Fig 1, you should be able to construct the ordinary differential equations.

For instance,

$\frac{d[C_1]}{dt} = (L_1\times C_4 + k_1\times C_2) - (3k_2\times A\times C_1 + L_4\times C_1)$

• In here $L_4$ is not listed in both table and figure but due to the computation, it is necessary to introduce to stabilize the calculation
• There are two L3 values. You can try either one of them and decide which one works best by comparing the plots from literature given below in Step 4
In [1]:
# Function for integreating ODE and displaying results
import scipy.integrate
from scipy.integrate import odeint
import numpy as np
import matplotlib.pylab as plt
import math
%matplotlib inline


### Step 1: define "ordinary differential equations" based on Fig 1¶

In [2]:
def func(y,t,ATP,t_end):
C1, C2, C3, C4, Q1, Q2, Q3, Q4 = y

### List your parameters such as k1, k2, ....
k1 = 0.3
k2 = 40000
k3 = 2.4
k4 = 50000
k5 = 1.58
k6 = 7000
L1 = 0.0001
L2 = 0.004
L3 = 0.1
L4 = 0.00000000001

### ATP switch
rep = math.floor(t/40)
if (rep % 2) == 0:
A = ATP
elif (rep % 2) == 1:
A = 0

### ODEs: There should be 8 ordinary differential equations for each state starting from C1, C2, ...
dC1dt = L1*C4 + k1*C2       - (3*k2*A + L4)*C1
dC2dt = 3*k2*A*C1 + 2*k3*Q1 - (k1 + 2*k4*A)*C2
dC3dt = 2*k1*Q4 + 3*k2*A*C4 - (k1 + 2*k2*A)*C3
dC4dt = L4*C1 + k1*C3       - (L1 + 3*k2*A)*C4
dQ1dt = 2*k4*A*C2 + 3*k5*Q2 - (2*k3 + k6*A)*Q1
dQ2dt = k6*A*Q1 + L2*Q3     - (L3 + 3*k5)*Q2
dQ3dt = L3*Q2 + k2*A*Q4     - (L2 + 3*k1)*Q3
dQ4dt = 3*k1*Q3 + 2*k2*A*C3 - (k2*A + 2*k1)*Q4

dydt = [dC1dt, dC2dt, dC3dt, dC4dt, dQ1dt, dQ2dt, dQ3dt, dQ4dt]

return dydt


### Step 2: solve ODE¶

In [8]:
## Initial state
y0 = [1,0,0,0,0,0,0,0] ## you need 8 initial states (C1, C2, C3, C4, Q1, Q2, Q3, Q4)

## Time: Ends abount 3 mins after
## Be careful with unit
t_start = 0
t_end   = 120
t_step  = t_end*10
t = scipy.linspace(t_start,t_end,t_step)

## Additional information such as ATP concentration: Test 10, 32, 100 uM.
## Be careful with unit!
## In future, this ATP stimulation would be no longer consistent and we will try.
ATP = 32*10e-6

## Solve!
y = scipy.integrate.odeint(func,y0,t,args=(ATP,t_end,))

<ipython-input-8-247f49d6c430>:9: DeprecationWarning: scipy.linspace is deprecated and will be removed in SciPy 2.0.0, use numpy.linspace instead
t = scipy.linspace(t_start,t_end,t_step)


### Step 3: converting open probability to current¶

• I : inward current
• $Q_{1,2,3,4}$ : open probability from calculated "y"... make sure you understand the index corresponding to these dependent variables and properly extract them to calculate the inward current.
• Tip: if you wonder the structure of "y", type print(np.shape(y)) or print(y) to check what columes and rows indicate
In [9]:
## Parameters
g12 = 1.5e-8
g34 = 4.5e-8
V   = -60e-3
E   = 0

I = g12*(y[:,4] + y[:,5])*(V-E) + g34*(y[:,6] + y[:,7])*(V-E)


### Step 4: Plot and compare them against literature¶

#### Fig 2: Plot from literature¶

In [10]:
plt.figure(figsize=(6,2),dpi=150)
plt.tick_params(labelsize=10,direction='in')
plt.plot(t,I)
plt.tight_layout()

In [11]:
print(t)

[0.00000000e+00 1.00083403e-01 2.00166806e-01 ... 1.19799833e+02
1.19899917e+02 1.20000000e+02]

In [12]:
print(I)

[-0.00000000e+00 -7.00648633e-10 -7.99048803e-10 ... -2.35437883e-09
-2.35479379e-09 -2.35520733e-09]

In [ ]: