#Readme: The author recommend to run the code in silent mode. The data will be output in test.txt. 
#When to do parameter sensitivity analysis, use command 'RANGE=1, RANGEOVER=step, RANGESTEP=1100, RANGELOW=0, RANGEHIGH=1100, RANGERESET=yes, RANGEOLDIC=yes, output=test1'. The data will be output in 1100 data file named test1.1 (step value), test1.2........

#A is BDNF activity, B is BDNF mRNA, CEBP is CEBP

dA/dt=(stim + k_1*B^2/(B^2+k_trs^2)*(1-PSI*turn)-kdega*A/(A+K_A2))

dB/dt=((K_min_B +k_2*feed3) - kdegb*B/(B+K_B2))

feed3=(CEBP^2)/(K_BCEBP^2+CEBP^2)*(heav(t2+3600*endh-t))+kbcedbnf*(heav(t-t2-3600*endh))

endh=18

BDNF=A


CREB1=CREBt-pCREB1

CREBt=0.085


pCREB1'=((K_phosC1+feed4)*CREB1-K_dephC1*pCREB1)

feed4=k_Cbdnf*BDNF^2/(BDNF^2+K_bpcreb^2)*(1-PSI*0)




#synthesis OF CEBP protein

CEBP'=((feed2*v_cebp+K_min_C)*(CEBPmax-CEBP)*(1-PSI*turn)-K_dxCBP*CEBP/(CEBP+K_dCEBP))

feed2 =(pCREB1^2)/(k_crece^2+pCREB1^2)*heav(t-t2-3600*starth)+kbpcrebce*heav(t2+3600*starth-t)

starth=6

CEBPmax=0.228


#W1 is pathway 1, W2 is pathway 2, W3 is synaptic weight
W1'=((kbasw1+kfw1*(pCREB1-0.0066)/kw1creb*(CEBP-0.076)/kw1cebp*w1^2/(w1^2+kw1^2))*(1-PSI*turn)-kdw1*w1)

w2'=((kbasw2+kfw2*(CEBP-0.076)/kw2cebp*w2^2/(w2^2+kw2^2))*(1-PSI*turn)-kdw2*w2)

w3'=(kfw3*(W1+w2)*(1-PSI*turn)-kdw3*w3)


# stimulus

stim = test1*HEAV(T-t2)*HEAV(dur+t2-T)+0*(1-HEAV(T-t2)*HEAV(dur+t2-T))
t2=10000 

# inhibitors

PSI=amp*heav(t-t1)*heav(t1+duration-t)

p amp=0

p turn=0


#ADD ODN
#p amp=0.95

#ADD anti-BDNF 
#p amp=0.8

#ADD PSI
#p amp=0.8

#duration for anti-BDNF
#duration=3600*6

#duration for ODN
#duration=3600*24

#duration for PSI
duration=3600*6


# time to start inhbit
t1=t2+3600*startIn
startIn=24

init A=0.0722924
init B=0.00320135
init pCREB1=0.0066
init CEBP=0.076
init w1=100
init w2=100
init w3=100



#parameter sensitivity analysis
K_BCEBP=1*(1+0.03*a4*heav(step+1)*heav(31-step)*HEAV(T-t2+7200))
k_trs=1*(1+0.03*a4*heav(step-30)*heav(62-step)*HEAV(T-t2+7200))
kdega=0.002*(1+0.03*a4*heav(step-61)*heav(93-step)*HEAV(T-t2+7200))
kdegb=0.0000226*(1+0.03*a4*heav(step-92)*heav(124-step)*HEAV(T-t2+7200))
K_bpcreb=5*(1+0.03*a4*heav(step-123)*heav(155-step)*HEAV(T-t2+7200))
K_B2=0.6*(1-0.03*a4*heav(step-154)*heav(186-step)*HEAV(T-t2+7200))
k_crece=2.24*(1+0.03*a4*heav(step-185)*heav(217-step)*HEAV(T-t2+7200))
K_min_B=0.000000078*(1-0.03*a4*heav(step-216)*heav(248-step)*HEAV(T-t2+7200))
k_1=21*(1-0.03*a4*heav(step-247)*heav(279-step)*HEAV(T-t2+7200))
K_dephC1=0.000024*(1+0.03*a4*heav(step-278)*heav(310-step)*HEAV(T-t2+7200))
K_dxCBP=0.00000565*(1+0.03*a4*heav(step-309)*heav(341-step)*HEAV(T-t2+7200))
kbcedbnf=0.00574*(1-0.03*a4*heav(step-340)*heav(372-step)*HEAV(T-t2+7200))
K_phosC1=0.0000008*(1-0.03*a4*heav(step-371)*heav(403-step)*HEAV(T-t2+7200))
k_2=0.00000726*(1-0.03*a4*heav(step-402)*heav(434-step)*HEAV(T-t2+7200))
k_Cbdnf=0.00592*(1-0.03*a4*heav(step-433)*heav(465-step)*HEAV(T-t2+7200))
K_A2=0.6*(1-0.03*a4*heav(step-464)*heav(496-step)*HEAV(T-t2+7200))
v_cebp=0.31*(1-0.03*a4*heav(step-495)*heav(527-step)*HEAV(T-t2+7200))
K_min_C=0.0000158*(1-0.03*a4*heav(step-526)*heav(558-step)*HEAV(T-t2+7200))
K_dCEBP=0.076*(1-0.03*a4*heav(step-557)*heav(589-step)*HEAV(T-t2+7200))
test1=0.03*(1-0.03*a4*heav(step-588)*heav(620-step)*HEAV(T-t2+7200))
dur=60*(1-0.03*a4*heav(step-619)*heav(651-step)*HEAV(T-t2+7200))
kbpcrebce=0.0000089*(1-0.03*a4*heav(step-650)*heav(682-step)*HEAV(T-t2+7200))
kbasw1=0.0004*(1-0.03*a4*heav(step-681)*heav(713-step)*HEAV(T-t2+7200))
kw1creb=0.198*(1+0.03*a4*heav(step-712)*heav(744-step)*HEAV(T-t2+7200))
kw1cebp=0.228*(1+0.03*a4*heav(step-743)*heav(775-step)*HEAV(T-t2+7200))
kfw1=1.2*(1-0.03*a4*heav(step-774)*heav(806-step)*HEAV(T-t2+7200))
kw1=150*(1+0.03*a4*heav(step-805)*heav(837-step)*HEAV(T-t2+7200))
kdw1=1/250000*(1+0.03*a4*heav(step-836)*heav(868-step)*HEAV(T-t2+7200))
kbasw2=1/100000*(1-0.03*a4*heav(step-867)*heav(899-step)*HEAV(T-t2+7200))
kw2cebp=0.228*(1+0.03*a4*heav(step-898)*heav(930-step)*HEAV(T-t2+7200))
kfw2=0.032*(1-0.03*a4*heav(step-929)*heav(961-step)*HEAV(T-t2+7200))
kw2=150*(1+0.03*a4*heav(step-960)*heav(992-step)*HEAV(T-t2+7200))
kdw2=1/10000000*(1+0.03*a4*heav(step-991)*heav(1023-step)*HEAV(T-t2+7200))
kdw3=0.00002*(1+0.03*a4*heav(step-1022)*heav(1054-step)*HEAV(T-t2+7200))
kfw3=0.00001*(1-0.03*a4*heav(step-1053)*heav(1085-step)*HEAV(T-t2+7200))

a4=(step-31*flr(step/31))

aux 0=(T-10000)/(3600*24)
aux 1=stim
aux 3=BDNF/0.072*100
aux 4=B/0.0032*100
aux 5=pCREB1/0.0066*100
aux 6=CEBP/0.076*100
aux 7=a4

p step=0


@ total=260000, xlo=0, xhi=5000, ylo=0, yhi=10, bounds=10e30, MAXSTOR=1300000,xp=vs, yp=mk,nout=2500, dt=0.05,

@ output=test.txt
#@ RANGE=1, RANGEOVER=step, RANGESTEP=1100, RANGELOW=0, RANGEHIGH=1100, RANGERESET=yes, RANGEOLDIC=yes, output=test1


d