pro Oscill_t

;Simple three-state ring model. violates microreversibility and generates some oscillation

;time in ms, concentrations in micromolar

wdelete,0 & wdelete,1

rl=400 ;record length

Po=dblarr(rl) ;open state probability

Pc=dblarr(rl) ;closed state

Pi=dblarr(rl) ;inactivated

pout=fltarr(4,rl)



g=fltarr(7,rl) ;output array



dummy=''

Euler=50. ;number of Euler increments per point

tint=0.25 ;milliseconds per point

t=tint*findgen(rl) ;time in ms,

Pout(0,*)=t(*) ;pass time to output array

dt=tint/Euler ;Euler increment

;************************************Parameters

kco=0.2 ;inverse ms

koc=0.0001 ;inverse ms

koi=0.2

kio=0.0001

kic=0.1

kci=0.0001



;************************************equations:

;

;************************************initial conditions



;for l=0,5 do begin ;(loop in off voltages)



Po0=0 & x1= 1 ;this is just a big jump in activation, from 0 to 1. Must include effect of voltage later.

Po(0) = 0

Pc(0)=0.99

Pi(0)=1-Po(0)-Pc(0)



;********************************** these values determined by voltage and duration of pulse

dur= 100. ;pulse duration in milliseconds

;********************************** computation of solution

;will go over time in two loops, one for milliseconds, one for Euler

tim=0.0 ;time

Pos=Po(0) & Pcs=Pc(0) & Pis=Pi(0)



for i=1,dur/tint <rl-1 do begin

for j=1,Euler do begin

tim=tim+dt

;xt is a scalar

dpo=(kco*Pcs+kio*Pis-(koc+koi)*Pos)*dt

dpc=(kic*Pis+koc*Pos-(kco+kci)*Pcs)*dt



Pos=Pos+dpo & Pcs=Pcs+dpc & Pis=1.0-Pos-Pcs

endfor

;now update arrays

t(i)=tim

Po(i)=Pos & Pc(i)=Pcs & Pi(i)=Pis

endfor

i2=i



kco=0.001



t0=tim

for i=i2,rl-1 do begin

for j=1,Euler do begin

tim=tim+dt

dpo=(kco*Pcs+kio*Pis-(koc+koi)*Pos)*dt

dpc=(kic*Pis+koc*Pos-(kco+kci)*Pcs)*dt



Pos=Pos+dpo & Pcs=Pcs+dpc & Pis=1.0-Pos-Pcs

endfor

;now update arrays

t(i)=tim

Po(i)=Pos & Pc(i)=Pcs & Pi(i)=Pis

endfor

;******************************display

window,0, xpos=0,ypos=0, title='thin, a, Po; red, Pi'

plot, t, Po, xrange=[0,tint*rl-1],xtitle='time, ms'

oplot, t, Pi, thick=2,color=200

oplot,t,Pc, thick=2,color=200,line=2

;pout(0,*)=t

Pout(1,*)=Pc & Pout(2,*)=Po & Pout(3,*)=Pi

fname='oscill3.asc'

openw,un,fname,/get_lun

printf,un,format='(f7.2,3f8.5)',pout

close,un & free_lun,un

stop

end