#!/usr/bin/env python
# coding: utf-8

# # Dual Phosphorylation Cycle

# In[1]:


get_ipython().run_line_magic('matplotlib', 'inline')
from ecell4.prelude import *


# In[2]:


citation(20133748)


# In[3]:


@species_attributes
def attrgen(radius, D):
    K | Kp | Kpp | KK | PP | K_KK | Kp_KK | Kpp_PP | Kp_PP | {"radius": radius, "D": D}

@reaction_rules
def rulegen(kon1, koff1, kcat1, kon2, koff2, kcat2):
    (K + KK == K_KK | (kon1, koff1)
        > Kp + KK | kcat1
        == Kp_KK | (kon2, koff2)
        > Kpp + KK | kcat2)

    (Kpp + PP == Kpp_PP | (kon1, koff1)
        > Kp + PP | kcat1
        == Kp_PP | (kon2, koff2)
        > K + PP | kcat2)


# In[4]:


radius, D = 0.0025, 1.0
ka1, kd1, kcat1 = 0.04483455086786913, 1.35, 1.5
ka2, kd2, kcat2 = 0.09299017957780264, 1.73, 15.0

m = NetworkModel()
m.add_species_attributes(attrgen(radius, D))
m.add_reaction_rules(rulegen(ka1, kd2, kcat1, ka2, kd2, kcat2))


# In[5]:


show(m)


# In[6]:


session = Session(model=m, y0={"K": 120, "KK": 30, "PP": 30})


# In[7]:


ret = session.run(60.0)


# In[8]:


ret.plot(y=["K+K_KK", "Kp+Kp_KK+Kp_PP", "Kpp+Kpp_PP"], legend=True)


# In[9]:


ret = session.run(60.0, ndiv=100, solver='gillespie')


# In[10]:


ret.plot(y=["K+K_KK", "Kp+Kp_KK+Kp_PP", "Kpp+Kpp_PP"], legend=True, step=True)