Markov State Modeling

Last updated: May 29th, 20202020-05-29Project preview

Markov State Modeling Practice

Purinergic Receptor, P2X7, and its gating mechanism based on ATP concentration

Link to the original paper

Fig 1: Markov State Model of P2X7

image

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

image

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

image

  • 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

image

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 [ ]:
 
Notebooks AI
Notebooks AI Profile20060