# Fig1.ode
#
# This XPPAUT file contains the mathematical model used in Figure 1 of
# Bertram, Tabak, Toporikova, and Freeman, Science STKE, in press, 2006.
# At time 1 minute a nanomolar concentration of endothelin is added.

# Variables (dimensionless): v=voltage
#                            n=delayed rectifier activation
#                            c=cytosolic calcium concentration 
#                            cer=ER calcium concentration
#                            inh=variable for direct inhibition of secretion
#                            cAMP=cAMP concentration

v(0)=-60
n(0)=0
c(0)=0.3
cer(0)=260
inh(0)=1
cAMP(0)=1

# Conductances in pS; currents in fA; Ca concentrations in uM; time in ms

# conductances
num vh=-20, sh=70, tauh=20
num kserca=0.4
# Ca parameters  (sigmav=cyt volume/ER volume)
num sigmav=10
par kc=0.15

# Ik
number vn=-16, vk=-75, taun=20
num gk=3500
num sn=5

# Ica
number vca=25
num gca=2000
number vm=-20, sm=12

# Miscellaneous
number lambda=1.25
number cm=5300

# Calcium Handling
number f=0.01, fer=0.01
number alpha=4.50e-6

# Calcium Handling: ER
num perl=0.0005, dact=0.35, dip3=0.5, dinh=0.4

# Functions

ninf = 1/(1+exp((vn-v)/sn))
minf = 1/(1+exp((vm-v)/sm))
hinf = 1/(1+exp((vh-v)/sh))
ica = gca*minf*(v-vca)
igirk = girk*h*(v-vk)
ik = gk*n*(v-vk)

% ET-1 switch
par tET=60000
ETswitch = heav(t-tet)

% IP3 concentration
par IP3mag=0.3
IP3 = IP3mag*ETswitch

% GIRK conductance
par girk_noET=1000, girk_ET=3000
girk=(1-ETswitch)*girk_noET + ETswitch*girk_ET

% ER functions
ainf = 1/(1 + dact/c)
hinfer = 1/(1 + c/dinh)
jerp = kserca*c
binf = IP3/(IP3 + dip3)
o = ainf^3*binf^3*hinfer^3

% Ca fluxes
jmemtot = -(alpha*ica + kc*c)
jerleak = perl*(cer - c)
jerip3 = o*(cer - c)
jertot = jerleak + jerip3 - jerp

# Prl secretion
par ki=0.5
num pow=4
prl_inf = inh*cAMP*c^pow/(ki^pow + c^pow)

# direct effect of ET on secretion
par taudir=20000

# Set cAMPlow=0.2 for nM ET, and cAMPlow=0.5 for pM ET
par cAMPlow=0.2

# Equations

v' = (-ica - ik - igirk)/cm
n' =  lambda*(ninf - n)/taun
h' = (hinf-h)/tauh
c' = f*(jertot + jmemtot)
cer' = -fer*sigmav*jertot
inh' = ETswitch*((0.2-inh)/taudir)
cAMP' = ETswitch*((cAMPlow-cAMP)/taudir)

aux ip3=IP3
aux girk=girk

### average of last qqq ms
par qqq=600
@ delay=6000
storP'=Prl_inf/qqq
aux smoothP=(storP-delay(storP,qqq))

@ meth=cvode, toler=1.0e-10, atoler=1.0e-10, dt=20.0, total=180000 
@ maxstor=20000,bounds=10000000, xp=tmin, yp=c
@ xlo=0, xhi=3, ylo=0, yhi=0.8, bell=off

aux tsec=t/1000
aux tmin=t/60000
aux PRL=Prl_inf

done