"""
ER Calcium model
from Wagner et al 2004: "A wave of IP3 production..."

vserca modified

Demonstrates intracellular compartments.

snnet version 2011-03-16 by Robert A McDougal
"""

from snnet import *

cell = """
    # total cell volume (including any compartments)
    volume = 1

    # IP3 dynamics
    vprod = .075 # uM / s  (in principle depends on surface area/volume ratio)
    kprod = .4   # uM
    v1 = .001    # uM / s
    v2 = .005    # uM / s
    v3 = .02     # uM / s
    lambda = 30
    k0 = .39     # uM
    k1 = 2.5     # uM
    k2 = .5      # uM
    k3 = 30      # uM

    # declare ca and ip3 to be a concentration (necessary for using flux)
    # this means that ca depends on effective volume
    # effective volume = total volume - volume organelles
    concentration ca, ip3
    
    # IP_3 degradation and production
    theta = ca / (ca + k0)
    jkinase = (1 - theta) * v1 * ip3 / (ip3 + k1) + theta * v2 * ip3 / (ip3 + k2)
    jphos = v3 * ip3 / (ip3 + k3)
    
    # vprod includes the effects of cell surface area / volume ratio
    jprod = vprod * ca ^ 2 / (ca ^ 2 + kprod ^ 2)
    
    # still need to divide by the effective_volume
    ip3' = (jprod - lambda * (jkinase + jphos)) / effective_volume

    # initial conditions
    ca = 1
    ip3 = 1.5
"""


er = """
    # calcium regulation
    vip3r = 8.5  # /s
    kip3 = .15   # uM
    kact = .8    # uM
    vleak = .01  # /s
    vserca = 1   # uM / s (was .65 in W et al)
    kserca = .4  # uM

    # channel state
    kinh = 1.9   # uM
    tau = 2      # s

    # volume of ER
    volume = .17

    concentration ca
    
    ip3 = external_ip3
    cai = external_ca

    hinf = kinh / (kinh + cai)
    m = ip3 / (ip3 + kip3) * cai / (cai + kact)
    jip3r = vip3r * m ^ 3 * h ^ 3 * (ca - cai)
    jleak = vleak * (ca - cai)
    jserca = vserca * cai ^ 2 / (kserca ^ 2 + cai ^ 2)

    # calcium flux (outward = positive)
    ca flux jip3r - jserca + jleak
     
    # IP_3 receptor inactivation
    h' = ((hinf - h) / tau)

    # initial conditions
    h = 0
    ca = 3.94117647
"""

def simulation():
    c = neurons(cell)
    c.add_compartment(er)
    advance(200)

data = run(simulation, integrator=fast_xpp_integrator)

subplot(2, 1, 1)
plot(data.t, data.ca0)
ylabel('[Ca$^{2+}$]$_c$ in $\mu$M')

subplot(2, 1, 2)
# ca is a property of the cell, which has index 0
# h is a property of the ER, which has index 1
plot(data.ca0, data.h1)
xlabel('[Ca$^{2+}$]$_c$ in $\mu$M')
ylabel('$h$')

savefig('example3.pdf')

