METHOD RK4

STARTTIME = 0
STOPTIME=60
DT = 0.001
init x2=9e-13
init x3=0
init x4=0
init x5=0
init x6=1e-13
init x7=0
init x8=0
init x9=1e-12
init x10=0
init x10a=0
init x11=1e-13
init x12=0
init x13=5e-8
init x14=1.6e-5
init x15=4.68e-8
init x16=1e-8
init x17=0
init x18=100
init x19=0
init x20= 1.5e-8
init x21= 0
init x22= 0
init x23=0.5e-6
init x24=0
init x25=12e-8
init x26=8.4e-9
init x27=13e-9
init x28=0
init x29=0.5e-6
init x30=0
init x31=0
init x32=0.5e-6
init x33=0
init x34=0
init x35=0
init x36=96
init x37=4
init Z=0

insstep=1e-7
stepdur=15.0

k1=6e7
km1=0.20
k2=k1
k3=2500.00
km2=100*km1
km3=km1
k4=km4/9
km4=0.003
k4p=2.1e-3
km4p=2.1e-4
k5=if((x6+x7+x8)>1e-13)then(10*km5)else(60*km5)
km5=1.66667e-18
k6=0.461
k7=4.16
km7=(2.5/7.45)*k7
k7p = if(negfeed)then(.346574)else(0)
km7p= if(negfeed)then(k7p*((2.5/7.45)*(3.7e-13))/((6.27e-13)-(2.5/7.45)*(3.7e-13)))else(0)
km8=10.0
k8=km8*5.00e12/70.775
k9stim=1.39
km9=(1.52e-5/5e-7)*k9stim
k9basal=(5e-8/1.6e-5)*km9
k9=(k9stim - k9basal)*(x12/PI3K) + k9basal
k10=km10*(5e-7/4.68e-7) 
km10=2.77
k11=0.1*km11*(x13-5e-8)/(5e-7-5e-8)
k12=0.1*km12*(x13-5e-8)/(5e-7-5e-8)
km11=6.93
km12=6.93
km13= .34
KMI13= 9e-9
k13=2 
k14=80
kMI14= 2e-6
k15= .6
k16= 1.386
km16= k16*(3.75/1.25)
V17= 2.3e-4
kMI17= 7.3e-6
k18= 13.2
kMI18=.18e-9
k19=360
kMI19= 1.5e-6
k20=60
kMI20 = 0.5e-9
k21=360
kMI21=15.7e-6
k22=50
kMI22= 1.4e-8
k23= 360
kMI23= 15.7e-6
k24 = 50
kMI24= 1.4e-8
k25=60
kMI25=.167e-6
k26=120
kMI26=1.5e-7
k27=60
kMI27= .167e-6
k28= 120
kMI28= 1.5e-7
k29=0.005
k30= (.693/15)


ARG= 1e-4
MKP3= 2.4e-9 
PI3K=  if(negfeed)then((k8*(3.7e-13)*(1e-13))/(k8*(3.7e-13)+km8))else(5.0e-15)
PKC= Vmax*delay(x19,tau)^n/(delay(x19,tau)^n + kd^n)
PTP=1 
PTEN= 1
SHIP= 1
PP2A=2.7e-6
MKP1= .0032e-6
IRp=8.97e-13
Vmax= 20
tau=1.5
kd=12
n=4
n2=1
ki= 1e-9
negfeed=1


x1 = IF (Z-stepdur)< 0 THEN insstep ELSE 0
d/dt(Z)=1
d/dt(x2) =  km1*x3 + km3*PTP*x5 - k1*x1*x2 + km4*x6 - k4*x2
d/dt(x3) =  k1*x1*x2 - km1*x3 - k3*x3
d/dt(x4) =  k2*x1*x5 - km2*x4 + km4p*x7 - k4p*x4
d/dt(x5) =  k3*x3 + km2*x4 - k2*x1*x5 - km3*PTP*x5 + km4p*x8 - k4p*x5
d/dt(x6) =  k5 - km5*x6 + k6*PTP*(x7 + x8) + k4*x2 - km4*x6
d/dt(x7) =  k4p*x4 - km4p*x7 - k6*PTP*x7
d/dt(x8) =  k4p*x5 - km4p*x8 - k6*PTP*x8
d/dt(x9) =  km7*PTP*x10 - k7*x9*(x4 + x5)/IRp + km7p*x10a - k7p*PKC*x9
d/dt(x10) = k7*x9*(x4 + x5)/IRp + km8*x12 - (km7*PTP + k8*x11)*x10
d/dt(x10a) = k7p*PKC*x9 - km7p*x10a
d/dt(x11) = km8*x12 - k8*x10*x11
d/dt(x12) = k8*x10*x11 - km8*x12
d/dt(x13) = k9*x14 + k10*x15 - (km9*PTEN + km10*SHIP)*x13
d/dt(x14) = km9*PTEN*x13 - k9*x14
d/dt(x15) = km10*SHIP*x13 - k10*x15
d/dt(x16) = km11*x17 - k11*x16
d/dt(x17) = k11*x16 - km11*x17
d/dt(x18) = km12*x19 - k12*x18
d/dt(x19) = k12*x18 - km12*x19



d/dt(x20)= km13*x21-k13*x17*x20/(kMI13+(x20))
d/dt(x21)= k13*x17*x20/(kMI13+(x20)) - km13*x21
d/dt(x22)= ((k14*x21*ARG)/(kMI14 + ARG))- k15*x22




d/dt(x23)= km16*x24 - k16*x23*(x4+x5)/IRp
d/dt(x24)= k16*(x23)*(x4+x5)/IRp - km16*(x24)
d/dt(x25)= V17*(x26)/(kMI17+(x26))-k18*(x24)*(x25)/(kMI18+(x25))
d/dt(x26)= k18*(x24)*(x25)/(kMI18+(x25))-V17*(x26)/(kMI17+(x26))
d/dt(x27)= k19*(x28)*(PP2A)/(kMI19+(x28))-k20*(x26)*(x27)/((kMI20+(x27))*(1 + (x17/Ki)^n2)) 
d/dt(x28)= k20*(x26)*(x27)/((kMI20+(x27))*(1 + (x17/Ki)^n2)) - k19*(x28)*(PP2A)/(kMI19+(x28))
d/dt(x29)= k21*(PP2A)*(x30)/(kMI21+(x30))-k22*(x28)*(x29)/(kMI22+(x29))
d/dt(x30)= k22*(x28)*(x29)/(kMI22+(x29))-k21*(PP2A)*(x30)/(kMI21+(x30))+k23*(PP2A)*(x31)/(kMI23+(x31))-k24*(x28)*(x30)/(kMI24+(x30))
d/dt(x31)= k24*(x28)*(x30)/(kMI24+(x30))-k23*(PP2A)*(x31)/(kMI23+(x31))
d/dt(x32)= k25*(MKP1)*(x33)/(kMI25+(x33))-k26*(x31)*(x32)/(kMI26+(x32))
d/dt(x33)= k26*(x31)*(x32)/(kMI26+(x32))+k27*(MKP1)*(x34)/(kMI27+(x34)) -k28*(x31)*(x33)/(kMI28+(x33))-k25*(MKP1)*(x33)/(kMI25+(x33))
d/dt(x34)= k28*(x31)*(x33)/(kMI28+(x33))-k27*(MKP1)*(x34)/(kMI27+(x34))
d/dt(x35)= k29*(x34+x33) - k30*x35



effect= 100* (.6*x22- 66.667*.4*x35)/(3.95e-7)
Vasctone= 100-effect