Axonal gap junctions produce fast oscillations in cerebellar Purkinje cells (Traub et al. 2008)

 Download zip file 
Help downloading and running models
Accession:114654
Examines how electrical coupling between proximal axons produces fast oscillations in cerebellar Purkinje cells. Traub RD, Middleton SJ, Knopfel T, Whittington MA (2008) Model of very fast (>75 Hz) network oscillations generated by electrical coupling between the proximal axons of cerebellar Purkinje cells. European Journal of Neuroscience.
Reference:
1 . Traub RD, Middleton SJ, Knopfel T, Whittington MA (2008) Model of very fast (greater than 75 Hz) network oscillations generated by electrical coupling between the proximal axons of cerebellar Purkinje cells. Eur J Neurosci 28:1603-16
Citations  Citation Browser
Model Information (Click on a link to find other models with that property)
Model Type: Realistic Network; Neuron or other electrically excitable cell;
Brain Region(s)/Organism:
Cell Type(s): Cerebellum Purkinje GABA cell;
Channel(s): I Na,p; I Na,t; I T low threshold; I p,q; I A; I K; I M; I h; I K,Ca; I R;
Gap Junctions: Gap junctions;
Receptor(s):
Gene(s):
Transmitter(s):
Simulation Environment: FORTRAN;
Model Concept(s): Oscillations;
Implementer(s): Traub, Roger D ; Middleton, S. J. ; Knopfel, T.; Whittington, M. A. ;
Search NeuronDB for information about:  Cerebellum Purkinje GABA cell; I Na,p; I Na,t; I T low threshold; I p,q; I A; I K; I M; I h; I K,Ca; I R;
/
axon_coupled_net
readme.txt
purkinje_bignet21.f
purkinje17.f
                            
! 21 March 2007, version of purkinje.f to run on a node.
c 13 March 2007, begin construction of Purkinje cell model, starting with bask.f
c Sources include Llinas & Sugimori; Roth & Hausser for cell structure and passive
c parameters; de Schutter & Bower & Miyasho et al. for many of the active properties;
c Roth & Hausser for anomalous rectifier; Martina et al for data on delayed rectifier.
c 1 soma cylinder (level 1), 24 smooth dendritic compartments, #2 - 25.
c compartments 2 & 3 are dendritic shaft (level 2), 4 - 25 rest of smooth dendrites (level 3.
c There are 44 "canonical treelets" for the spiny dendrites; 2 treelets attach to
c each of the level 3 smooth dendritic compartments.
c A spiny treelet has 12 compartments: 4 along the main root, then 4 in each of 2 symmetrical
c branches.  The spiny dendrites are compartments 26 - 553.
c Axon has 6 5-micron compartments, #554 - 559, level 0.

c Active conductances: gnaf (m cubed h) - axon will use shifted version. 
c                      gnap (m cubed)
c                      gcaP (m) P channels
c                      gcaT (m h) T channels
c                      gcaR (m h) R channels
c                      gar  (m)   anomolous rectifier
c                      gKDR (m fourth) delayed rectifier, leaving out inactivation
c                      gKM  (m) M current, also called persistent K
c                      gKA  (m fourth x h)  A current
c                      gKD  (m h)  D current
c                      gKC  - try my original structure
c                      gKAHP - try my original structure.
c  12 Active conductances in all

        PROGRAM PURKINJE

       INTEGER J1, I, J, K, L, O, ISEED, K1, thisno
       REAL*8  TIMTOT, Z, Z1, Z2, curr(559), c(559), DT, time
       REAL*8 TIMER, gettime, time1, time2

c CINV is 1/C, i.e. inverse capacitance
       real*8 v(559), chi(559), cinv(559), gL(559), 
     x gNaF(559), gNaP(559), gCaP(559), gCaT(559), gCaR(559),
     x gar(559), gKDR(559), gKM(559), gKA(559), gKD(559),
     x gKC(559), gKAHP(559)

        real*8 jacob(559,559), betchi(559), gam(0:559,0:559)

        real*8
     x  mnaf(559), hnaf(559),
     x  mnap(559),
     x  mcap(559),
     x  mcat(559), hcat(559),
     x  mcar(559), hcar(559),
     x  mar (559),
     x  mkdr(559),
     x  mkm(559),
     x  mka(559), hka(559),
     x  mkd(559), hkd(559),
     x  mkc(559),
     x  mkahp(559)


       real*8  gampa(559),gnmda(559),ggaba_a(559),cafor(559)

       real*8 kdr_shift

       real*8
     X alpham_naf(0:640),betam_naf(0:640),dalpham_naf(0:640),
     X   dbetam_naf(0:640),
     X alphah_naf(0:640),betah_naf(0:640),dalphah_naf(0:640),
     X   dbetah_naf(0:640),
     X alpham_nap(0:640),betam_nap(0:640),dalpham_nap(0:640),
     X   dbetam_nap(0:640),

     X alpham_cap(0:640),betam_cap(0:640),dalpham_cap(0:640),
     X   dbetam_cap(0:640),
     X alpham_cat(0:640),betam_cat(0:640),dalpham_cat(0:640),
     X   dbetam_cat(0:640),
     X alphah_cat(0:640),betah_cat(0:640),dalphah_cat(0:640),
     X   dbetah_cat(0:640),
     X alpham_car(0:640),betam_car(0:640),dalpham_car(0:640),
     X   dbetam_car(0:640),
     X alphah_car(0:640),betah_car(0:640),dalphah_car(0:640),
     X   dbetah_car(0:640),

     X alpham_ar(0:640), betam_ar(0:640), dalpham_ar(0:640),
     X   dbetam_ar(0:640),

     X alpham_kdr(0:640),betam_kdr(0:640),dalpham_kdr(0:640),
     X   dbetam_kdr(0:640),
     X alpham_km(0:640), betam_km(0:640), dalpham_km(0:640),
     X   dbetam_km(0:640),
     X alpham_ka(0:640), betam_ka(0:640),dalpham_ka(0:640) ,
     X   dbetam_ka(0:640),
     X alphah_ka(0:640), betah_ka(0:640), dalphah_ka(0:640),
     X   dbetah_ka(0:640),
     X alpham_kd(0:640), betam_kd(0:640),dalpham_kd(0:640) ,
     X   dbetam_kd(0:640),
     X alphah_kd(0:640), betah_kd(0:640), dalphah_kd(0:640),
     X   dbetah_kd(0:640),
     X alpham_kc(0:640), betam_kc(0:640), dalpham_kc(0:640),
     X   dbetam_kc(0:640) 
       real*8 vL,vk,vna,var,vca,vgaba_a

        INTEGER NEIGH(559,6), NNUM(559)

c START EXECUTION PHASE
          include 'mpif.h'
          call mpi_init (info)
          call mpi_comm_rank(mpi_comm_world, thisno, info)
          call mpi_comm_size(mpi_comm_world, nodes , info)
          time1 = gettime()

        time1 = gettime()

c       kdr_shift =  7.d0
        kdr_shift = 10.d0

       CALL   PURKINJE_SETUP
     X   (alpham_naf, betam_naf, dalpham_naf, dbetam_naf,
     X    alphah_naf, betah_naf, dalphah_naf, dbetah_naf,
     X    alpham_nap, betam_nap, dalpham_nap, dbetam_nap,

     X    alpham_cap, betam_cap, dalpham_cap, dbetam_cap,
     X    alpham_cat, betam_cat, dalpham_cat, dbetam_cat,
     X    alphah_cat, betah_cat, dalphah_cat, dbetah_cat,
     X    alpham_car, betam_car, dalpham_car, dbetam_car,
     X    alphah_car, betah_car, dalphah_car, dbetah_car,

     X    alpham_ar , betam_ar , dalpham_ar , dbetam_ar,
     X    kdr_shift,
     X    alpham_kdr, betam_kdr, dalpham_kdr, dbetam_kdr,
     X    alpham_km , betam_km , dalpham_km , dbetam_km ,
     X    alpham_ka , betam_ka , dalpham_ka , dbetam_ka ,
     X    alphah_ka , betah_ka , dalphah_ka , dbetah_ka ,
     X    alpham_kd , betam_kd , dalpham_kd , dbetam_kd ,
     X    alphah_kd , betah_kd , dalphah_kd , dbetah_kd ,
     X    alpham_kc , betam_kc , dalpham_kc , dbetam_kc) 

          

        CALL PURKINJE_MAJ (GL,GAM,
     X    gNaF, gNaP, gCaP, gCaT, gCaR, gar, gKDR,
     X    gKM, gKA, gKD, gKC, gKAHP,
     X    CAFOR,JACOB,C,BETCHI,NEIGH,NNUM)

          do i = 1, 559
             cinv(i) = 1.d0 / c(i)
          end do

c                 gnaf = 0.d0
c                 gnap = 0.d0
c                 gkdr = 0.d0
c                 gka = 0.d0
c                 gkd = 0.d0
c                 gkm = 0.d0
c                 gkc = 0.d0
c                 gkahp = 0.d0
c                 gcat = 0.d0
c                 gcaP = 0.d0
c                 gcaR = 0.d0
c                 gar = 0.d0

c Below introduces somatic shunt
c            gL(1) = 100.d0 * gL(1)
             gL(1) =   5.d0 * gL(1)

c below 4 statements disconnect axon from soma.
c            gam(1,554) = 0.d0
c            gam(554,1) = 0.d0
c            jacob(1,554) = 0.d0
c            jacob(554,1) = 0.d0

        MG = 1.0d0
C  IN MILLIMOLAR

        VL = -80.d0
        VK =  - 85.d0
        VNA = 45.d0
        VCA = 135.d0
        VAR = -30.d0
        VGABA_A = -75.d0

        DT = .00060d0
c       DT = .00200d0

        TIMTOT = 1000.00d0
c       TIMTOT = 500.d0 * dt

c ? initialize membrane state variables?
        v = VL

        k1 = idnint (4.d0 * (v(1) + 120.d0))

      hnaf = alphah_naf(k1)/(alphah_naf(k1)+betah_naf(k1))
      hka = alphah_ka(k1)/(alphah_ka(k1)+betah_ka(k1))
      hkd = alphah_kd(k1)/(alphah_kd(k1)+betah_kd(k1))
      hcat=alphah_cat(k1)/(alphah_cat(k1)+betah_cat(k1))
      hcar=alphah_car(k1)/(alphah_car(k1)+betah_car(k1))

       O = 0
       TIME = 0.d0
       TIMER = 0.d0 ! for periodic spikes

2000      TIME = TIME + DT
          O = O + 1
          IF (TIME .GT. TIMTOT) GO TO 2001

          if (mod(O,8000).eq.0) timer = timer + 4.8d0

c         curr(1) = 1.50d0
c         curr(1) = 0.75d0

c         IF  (MOD(O,100).EQ. 0 )                      THEN
c         IF  ( i.eq.i          )                      THEN
c           WRITE (6,904) TIME, v(1), v(14), chi(14),
c    X V(25), v(554), v(557), v(558), v(559), curr(559)
c    X ,curr(14)
c         ENDIF
904       FORMAT(2X,F8.3,2X,8f8.2,f8.4,f9.4)

           if (time.le. 50.d0) then
          curr(1) = -1.00d0
           else if (time.le. 950.d0) then
          curr(1) = 1.50d0
           else
          curr(1) = -1.00d0
           endif

c           if (time.le. 50.d0) then
c         curr(14) = -0.1d0
c           else if (time.le.450.d0) then
c         curr(14) = 0.15d0
c           else
c         curr(14) = -0.1d0
c           endif

c Single antidromic spike
c          if (time.le.20.d0) then
c         curr(559) =  0.00d0
c          else if (time.le.20.50d0) then
c         curr(559) = 0.40d0 * (time-20.d0) / 0.5d0
c          else if (time.le.21.d0) then
c      curr(559) = 0.4d0 - 0.40d0 * (time-20.5d0) / 0.5d0
c          else
c         curr(559) = 0.d0
c          endif

c Periodic antidromic spikes
c            if (mod(O, 8000).le.900) then
c         curr(559) = 0.4d0 * (time-timer) / 0.54d0
c            else if (mod(O, 8000).lt.1800) then
c      curr(559) = 0.4d0 -0.4d0*(time-timer-0.54d0)/0.54d0
c            else
c          curr(559) = 0.0d0
c            endif

c            if (mod(O, 8000).le.1300) then
c         curr(559) = 0.5d0
c            else
c          curr(559) = 0.15d0
c            endif

909                  CONTINUE

         CALL   PURKINJE_INT (V,CHI,CINV,
     X mnaf,hnaf,mnap,mcap,mcat,hcat,mcar,hcar,mar,
     x mkdr,mkm,mka,hka,mkd,hkd,mkc,mkahp,
     x dt,neigh,nnum,jacob,mg,
     x vL,vk,vna,var,vca,vgaba_a,betchi,gam,gL,
     x gnaf,gnap,gcaP,gcaT,gCaR,gar,
     x gKDR,gKM,gKA,gKD,gKC,gKAHP,     
     x gampa,gnmda,ggaba_a,
     x O,time,
     X    alpham_naf, betam_naf, dalpham_naf, dbetam_naf,
     X    alphah_naf, betah_naf, dalphah_naf, dbetah_naf,
     X    alpham_nap, betam_nap, dalpham_nap, dbetam_nap,

     X    alpham_cap, betam_cap, dalpham_cap, dbetam_cap,
     X    alpham_cat, betam_cat, dalpham_cat, dbetam_cat,
     X    alphah_cat, betah_cat, dalphah_cat, dbetah_cat,
     X    alpham_car, betam_car, dalpham_car, dbetam_car,
     X    alphah_car, betah_car, dalphah_car, dbetah_car,

     X    alpham_ar , betam_ar , dalpham_ar , dbetam_ar,

     X    alpham_kdr, betam_kdr, dalpham_kdr, dbetam_kdr,
     X    alpham_km , betam_km , dalpham_km , dbetam_km ,
     X    alpham_ka , betam_ka , dalpham_ka , dbetam_ka ,
     X    alphah_ka , betah_ka , dalphah_ka , dbetah_ka ,
     X    alpham_kd , betam_kd , dalpham_kd , dbetam_kd ,
     X    alphah_kd , betah_kd , dalphah_kd , dbetah_kd ,
     X    alpham_kc , betam_kc , dalpham_kc , dbetam_kc ,
     X    cafor,curr)


              GO TO 2000
2001          CONTINUE

         time2 = gettime()
         write(6,309) time2 - time1
309      format(' Elapsed time = ',f8.0,' secs')


1000    CONTINUE
        call mpi_finalize(info)
        END



C  SETS UP TABLES FOR RATE FUNCTIONS
       SUBROUTINE PURKINJE_SETUP
     X   (alpham_naf, betam_naf, dalpham_naf, dbetam_naf,
     X    alphah_naf, betah_naf, dalphah_naf, dbetah_naf,
     X    alpham_nap, betam_nap, dalpham_nap, dbetam_nap,

     X    alpham_cap, betam_cap, dalpham_cap, dbetam_cap,
     X    alpham_cat, betam_cat, dalpham_cat, dbetam_cat,
     X    alphah_cat, betah_cat, dalphah_cat, dbetah_cat,
     X    alpham_car, betam_car, dalpham_car, dbetam_car,
     X    alphah_car, betah_car, dalphah_car, dbetah_car,

     X    alpham_ar , betam_ar , dalpham_ar , dbetam_ar,
     X    kdr_shift,
     X    alpham_kdr, betam_kdr, dalpham_kdr, dbetam_kdr,
     X    alpham_km , betam_km , dalpham_km , dbetam_km ,
     X    alpham_ka , betam_ka , dalpham_ka , dbetam_ka ,
     X    alphah_ka , betah_ka , dalphah_ka , dbetah_ka ,
     X    alpham_kd , betam_kd , dalpham_kd , dbetam_kd ,
     X    alphah_kd , betah_kd , dalphah_kd , dbetah_kd ,
     X    alpham_kc , betam_kc , dalpham_kc , dbetam_kc) 

      INTEGER I,J,K
      real*8 minf, hinf, taum, tauh, V, Z,
     X alpham_naf(0:640),betam_naf(0:640),dalpham_naf(0:640),
     X   dbetam_naf(0:640),
     X alphah_naf(0:640),betah_naf(0:640),dalphah_naf(0:640),
     X   dbetah_naf(0:640),
     X alpham_nap(0:640),betam_nap(0:640),dalpham_nap(0:640),
     X   dbetam_nap(0:640),
  

     X alpham_caP(0:640),betam_caP(0:640),dalpham_caP(0:640),
     X   dbetam_caP(0:640),
     X alpham_cat(0:640),betam_cat(0:640),dalpham_cat(0:640),
     X   dbetam_cat(0:640),
     X alphah_cat(0:640),betah_cat(0:640),dalphah_cat(0:640),
     X   dbetah_cat(0:640),
     X alpham_car(0:640),betam_car(0:640),dalpham_car(0:640),
     X   dbetam_car(0:640),
     X alphah_car(0:640),betah_car(0:640),dalphah_car(0:640),
     X   dbetah_car(0:640),

     X alpham_ar(0:640), betam_ar(0:640), dalpham_ar(0:640),
     X   dbetam_ar(0:640),

     X alpham_kdr(0:640),betam_kdr(0:640),dalpham_kdr(0:640),
     X   dbetam_kdr(0:640),
     X alpham_km(0:640), betam_km(0:640), dalpham_km(0:640),
     X   dbetam_km(0:640),
     X alpham_ka(0:640), betam_ka(0:640),dalpham_ka(0:640) ,
     X   dbetam_ka(0:640),
     X alphah_ka(0:640), betah_ka(0:640), dalphah_ka(0:640),
     X   dbetah_ka(0:640),
     X alpham_kd(0:640), betam_kd(0:640),dalpham_kd(0:640) ,
     X   dbetam_kd(0:640),
     X alphah_kd(0:640), betah_kd(0:640), dalphah_kd(0:640),
     X   dbetah_kd(0:640),
     X alpham_kc(0:640), betam_kc(0:640), dalpham_kc(0:640),
     X   dbetam_kc(0:640) 

       real*8  kdr_shift

C FOR VOLTAGE, RANGE IS -120 TO +40 MV (absol.), 0.25 MV RESOLUTION

       DO 1, I = 0, 640
          V = dble (I)
          V = (V / 4.d0) - 120.d0

c gNaF
c T. Miyasho et al. 2001 
           alpham_naf(i) = 35.d0 / (0.d0 + dexp((v + 5.d0)/(-10.d0))) 
           betam_naf(i) =  7.d0 / (0.d0 + dexp((v+65.d0)/20.d0))

c T. Miyasho et al. 2001
            alphah_naf(i) = 0.225d0 / (1.d0 + dexp((v+80.d0)/10.d0)) 
            betah_naf(i) =  7.5d0 / (0.d0 + dexp((v-3.d0)/(-18.d0)))

c T. Miyasho et al. 2001
            alpham_nap(i) = 200.d0  / (1.d0 + dexp((v-18.d0)/(-16.d0))) 
            betam_nap(i) =  25.d0 / (1.d0 + dexp((v+58.d0)/8.d0))

c P channels, T. Miyasho et al. 2001
            alpham_cap(i) = 8.5d0 / (1.d0 + dexp((v-8.d0)/(-12.5d0)))  
            betam_cap(i) =  35.d0 / (1.d0 + dexp((v+74.d0)/14.5d0))

c T-current, from T. Miyasho et al., 2001
            alpham_cat(i) = 2.6d0 / (1.d0 + dexp((v+21.d0)/(-8.d0))) 
            betam_cat(i) =  0.18d0 / (1.d0 + dexp((v+40.d0)/4.d0))
            alphah_cat(i) = 0.0025d0 / (1.d0 + dexp((v+40.d0)/8.d0))
            betah_cat(i) =  0.19d0 / (1.d0 + dexp((v+50.d0)/(-10.d0)))

c R-current, from T. Miyasho et al., 2001
            alpham_car(i) = 2.6d0 / (1.d0 + dexp((v+7.d0)/(-8.d0))) 
            betam_car(i) =  0.18d0/ (1.d0 + dexp((v+26.d0)/4.d0))
            alphah_car(i) = 0.0025/ (1.d0 + dexp((v+32.d0)/8.d0))
            betah_car(i) =  0.19d0/ (1.d0 + dexp((v+42.d0)/(-10.d0)))

c h-current (anomalous rectifier), Roth and Hausser 
           alpham_ar(i) = 0.00063d0 * dexp (-0.063d0 * (v + 73.2d0))
           betam_ar(i) =  0.00063d0 * dexp(0.079d0 * (v +73.2d0))

c delayed rectifier, non-inactivating. Start with mine (e.g. bask.f)
c Perhaps modify after comparison with Martina et al., 2003
           v = v + kdr_shift
c      minf = 1.d0/(1.d0+dexp((-v-27.d0)/11.5d0))
       minf = 1.d0/(1.d0+dexp((-v-20.d0)/11.5d0)) ! elevate threshold
            if (v.le.-10.d0) then
             taum = .25d0 + 4.35d0*dexp((v+10.d0)/10.d0)
            else
             taum = .25d0 + 4.35d0*dexp((-v-10.d0)/10.d0)
            endif
              alpham_kdr(i) = minf / taum
              betam_kdr(i) = 1.d0 /taum - alpham_kdr(i)
c from Martina, Schultz et al., 1998 - for bask.f
           v = v - kdr_shift

c M-current, from plast.f, with 60 mV shift
        alpham_km(i) = .02d0/(1.d0+dexp((-v-20.d0)/5.d0))
        betam_km(i) = .01d0 * dexp((-v-43.d0)/18.d0)

c A current: T. Miyasho et al., 2001
           alpham_ka(i) = 1.4d0 / (1.d0 + dexp((v+27.d0)/(-12.d0))) 
           betam_ka(i) =  0.49d0/ (1.d0 + dexp((v+30.d0)/4.d0))
           alphah_ka(i) = 0.0175d0/(1.d0 + dexp((v+50.d0)/8.d0))
           betah_ka(i) = 1.3 / (1.d0 + dexp((v+13.d0)/(-10.d0)))

c D current: T. Miyasho et al., 2001
           alpham_kd(i) = 8.5d0 / (1.d0 + dexp((v+17.d0)/(-12.5d0))) 
           betam_kd(i) =  35.0d0/ (1.d0 + dexp((v+99.d0)/14.5d0))
           alphah_kd(i) = 0.0015d0/(1.d0 + dexp((v+89.d0)/8.d0))
           betah_kd(i) = 0.0055d0 / (1.d0 + dexp((v+83.d0)/(-8.d0)))

c voltage part of C-current, using 1994 kinetics, shift 60 mV
              if (v.le.-10.d0) then
       alpham_kc(i) = (2.d0/37.95d0)*dexp((v+50.d0)/11.d0 -
     x                                     (v+53.5)/27.d0)
       betam_kc(i) = 2.d0*dexp((-v-53.5d0)/27.d0)-alpham_kc(i)
               else
       alpham_kc(i) = 2.d0*dexp((-v-53.5d0)/27.d0)
       betam_kc(i) = 0.d0
               endif
c Speed-up of C kinetics here.
          alpham_kc(i) = 2.d0 * alpham_kc(i)
           betam_kc(i) = 2.d0 *  betam_kc(i)

1        CONTINUE

         do 2, i = 0, 639

      dalpham_naf(i) = (alpham_naf(i+1)-alpham_naf(i))/.25d0
      dbetam_naf(i) = (betam_naf(i+1)-betam_naf(i))/.25d0
      dalphah_naf(i) = (alphah_naf(i+1)-alphah_naf(i))/.25d0
      dbetah_naf(i) = (betah_naf(i+1)-betah_naf(i))/.25d0
      dalpham_nap(i) = (alpham_nap(i+1)-alpham_nap(i))/.25d0
      dbetam_nap(i) = (betam_nap(i+1)-betam_nap(i))/.25d0

      dalpham_cap(i) = (alpham_cap(i+1)-alpham_cap(i))/.25d0
      dbetam_cap(i) = (betam_cap(i+1)-betam_cap(i))/.25d0
      dalpham_cat(i) = (alpham_cat(i+1)-alpham_cat(i))/.25d0
      dbetam_cat(i) = (betam_cat(i+1)-betam_cat(i))/.25d0
      dalphah_cat(i) = (alphah_cat(i+1)-alphah_cat(i))/.25d0
      dbetah_cat(i) = (betah_cat(i+1)-betah_cat(i))/.25d0
      dalpham_car(i) = (alpham_car(i+1)-alpham_car(i))/.25d0
      dbetam_car(i) = (betam_car(i+1)-betam_car(i))/.25d0
      dalphah_car(i) = (alphah_car(i+1)-alphah_car(i))/.25d0
      dbetah_car(i) = (betah_car(i+1)-betah_car(i))/.25d0

      dalpham_ar(i) = (alpham_ar(i+1)-alpham_ar(i))/.25d0
      dbetam_ar(i) = (betam_ar(i+1)-betam_ar(i))/.25d0

      dalpham_kdr(i) = (alpham_kdr(i+1)-alpham_kdr(i))/.25d0
      dbetam_kdr(i) = (betam_kdr(i+1)-betam_kdr(i))/.25d0
      dalpham_km(i) = (alpham_km(i+1)-alpham_km(i))/.25d0
      dbetam_km(i) = (betam_km(i+1)-betam_km(i))/.25d0
      dalpham_ka(i) = (alpham_ka(i+1)-alpham_ka(i))/.25d0
      dbetam_ka(i) = (betam_ka(i+1)-betam_ka(i))/.25d0
      dalphah_ka(i) = (alphah_ka(i+1)-alphah_ka(i))/.25d0
      dbetah_ka(i) = (betah_ka(i+1)-betah_ka(i))/.25d0
      dalpham_kd(i) = (alpham_kd(i+1)-alpham_kd(i))/.25d0
      dbetam_kd(i) = (betam_kd(i+1)-betam_kd(i))/.25d0
      dalphah_kd(i) = (alphah_kd(i+1)-alphah_kd(i))/.25d0
      dbetah_kd(i) = (betah_kd(i+1)-betah_kd(i))/.25d0
      dalpham_kc(i) = (alpham_kc(i+1)-alpham_kc(i))/.25d0
      dbetam_kc(i) = (betam_kc(i+1)-betam_kc(i))/.25d0
2      CONTINUE
       END

       SUBROUTINE PURKINJE_MAJ (GL,GAM,
     X    gNaF, gNaP, gCaP, gCaT, gCaR, gar, gKDR,
     X    gKM, gKA, gKD, gKC, gKAHP,
     X    CAFOR,JACOB,C,BETCHI,NEIGH,NNUM)

c Conductances: see main program.

c Note VAR = equil. potential for anomalous rectifier.

c 1 soma cylinder (level 1), 24 smooth dendritic compartments, #2 - 25.
c compartments 2 & 3 are dendritic shaft (level 2), 4 - 25 rest of smooth dendrites (level 3.
c There are 44 "canonical treelets" for the spiny dendrites; 2 treelets attach to
c each of the level 3 smooth dendritic compartments (22 of them).
c A spiny treelet has 12 compartments: 4 along the main root, then 4 in each of 2 symmetrical
c branches.  The spiny dendrites are compartments 26 - 553, level 4.
c Axon has 6 5-micron compartments, #554 - 559, level 0.

c Drop "glc"-like terms, just using "gl"-like
c CAFOR corresponds to "phi" in Traub et al., 1994
c Consistent set of units: nF, mV, ms, nA, microS

        REAL*8 C(559),GL(559),GAM(0:559,0:559)
        REAL*8 gnaf(559), gnap(559), gcap(559), gcat(559), gcar(559)
        REAL*8 gar(559), gKDR(559), gKM(559), gKA(559), gKD(559),
     X   gKC(559), gKAHP(559) 
        REAL*8 JACOB(559,559),RI_SD,RI_AXON,RM_SD,RM_AXON,CDENS
        INTEGER LEVEL(559)
        REAL*8 GNAF_DENS(0:4), GNAP_DENS(0:4) 
        REAL*8 GCAP_DENS(0:4), GCAT_DENS(0:4), GCAR_DENS(0:4) 
        REAL*8 GAR_DENS(0:4)
        REAL*8 GKDR_DENS(0:4), GKM_DENS(0:4), GKA_DENS(0:4)
        REAL*8 GKD_DENS(0:4), GKC_DENS(0:4), GKAHP_DENS(0:4)

        REAL*8 RES, RINPUT, z1, z2, z, somaar
        REAL*8 RSOMA, PI, BETCHI(559), CAFOR(559)
        REAL*8 RAD(559), LEN(559), GAM1, GAM2, ELEN(559)
        REAL*8 RIN, D(559), AREA(559), RI
        INTEGER NEIGH(559,6), NNUM(559)
C FOR ESTABLISHING TOPOLOGY OF COMPARTMENTS

        RI_SD = 115.d0 ! Roth & Hausser
c       RM_SD = 125000.d0 ! Roth & Hausser
c       RM_SD = 100000.d0 ! close to Roth & Hausser
        RM_SD =  50000.d0 ! close to Miyasho et al.
        RI_AXON = 100.d0
c       RM_AXON = 125000.d0 ! NOTE same as for SD
        RM_AXON =   2000.d0 
        CDENS = 0.8d0 ! Roth & Hausser

        PI = 3.14159d0

c See T. Miyasho et al., Brain Research 2001, for initial densities, but
c there will be changes.  e.g. Ca spikes have voltage-dependent fast-repol.
c so I will make dendritic gKDR larger and gKC smaller.
c ALSO, since there is now axon with voltage-shifted gNaF activation, I will
c try much smaller gNaF somatic density (Miyasho et al. have 10,000 ohm-cm-sqd)
c gNaF
           gnaf_dens(0) = 3500.d0
           gnaf_dens(1) = 5000.d0
c          gnaf_dens(1) = 7000.d0
           gnaf_dens(2) =  10.d0
           gnaf_dens(3) =   0.d0
           gnaf_dens(4) =   0.d0
c gNaP
           gnap_dens(0) = 0.1d0
           gnap_dens(1) = 5.0d0
           gnap_dens(2) = 1.0d0
           gnap_dens(3) = 0.0d0
           gnap_dens(4) = 0.d0
c gCaP
           gcaP_dens(0) = 0.0d0
           gcaP_dens(1) = 0.0d0
           gcaP_dens(2) = 0.0d0
           gcaP_dens(3) = 5.0d0
           gcaP_dens(4) = 5.0d0
c gCaT
           gcaT_dens(0) = 0.0d0
           gcaT_dens(1) = 0.0d0
           gcaT_dens(2) = 0.5d0
           gcaT_dens(3) = 1.5d0
           gcaT_dens(4) = 1.5d0
c gCaR
           gcaR_dens(0) = 0.0d0
           gcaR_dens(1) = 0.0d0
           gcaR_dens(2) = 0.0d0
           gcaR_dens(3) = 8.0d0
           gcaR_dens(4) = 8.0d0

c gar - Miyasho only have Ih at soma, but many princ cells have dendritic Ih
           gar_dens(0) = 0.0d0
           gar_dens(1) = 0.005d0
           gar_dens(2) = 0.005d0
           gar_dens(3) = 0.005d0
           gar_dens(4) = 0.005d0

c gKDR
c          gKDR_dens(0) = 3500.0d0
           gKDR_dens(0) = 1000.0d0
c          gKDR_dens(1) = 7000.0d0
           gKDR_dens(1) = 1000.0d0
           gKDR_dens(2) = 0.5d0
           gKDR_dens(3) = 0.5d0
           gKDR_dens(4) = 0.5d0
c gKM
           gKM_dens(0) = 1.00d0
           gKM_dens(1) = 1.00d0
           gKM_dens(2) = 1.00d0
           gKM_dens(3) = 0.04d0
           gKM_dens(4) = 0.04d0
c gKA
           gKA_dens(0) =  1.0d0
           gKA_dens(1) = 15.0d0
           gKA_dens(2) = 80.0d0
           gKA_dens(3) = 80.0d0
           gKA_dens(4) = 80.0d0
c gKD
           gKD_dens(0) =  0.0d0
           gKD_dens(1) =  0.0d0
           gKD_dens(2) = 80.0d0
           gKD_dens(3) = 80.0d0
c          gKD_dens(3) = 200.0d0
           gKD_dens(4) = 80.0d0
c          gKD_dens(4) = 200.0d0
c gKC
           gKC_dens(0) =  0.0d0
           gKC_dens(1) =  0.0d0
           gKC_dens(2) = 25.0d0
           gKC_dens(3) = 25.0d0
           gKC_dens(4) = 25.0d0
c gKAHP
         gkahp_dens(0) =   0.d0
         gkahp_dens(1) =   0.d0
         gkahp_dens(2) =   0.d0
         gkahp_dens(3) = 1.60d0
         gkahp_dens(4) = 1.60d0

        WRITE   (6,9988)
9988    FORMAT(2X,'I',4X,'NADENS',' CADENS(P)',' KDRDEN',' KAHPDE',
     X     ' KCDENS',' KADENS')
        DO 9989, I = 0, 4
          WRITE (6,9990) I, gnaf_dens(i), gcaR_dens(i), gkdr_dens(i),
     X  gkahp_dens(i), gkc_dens(i), gka_dens(i)
9990    FORMAT(2X,I2,2X,F6.2,1X,F6.2,1X,F6.2,1X,F6.2,1X,F6.2,1X,F6.2)
9989    CONTINUE


        level(1) = 1

        level(2) = 2
        level(3) = 2

        do i = 4, 25
         level(i) = 3
        end do

        do i = 26, 553
         level(i) = 4
        end do

        do i = 554, 559
         level(i) = 0
        end do

c connectivity of axon
        nnum(554) = 2
        nnum(555) = 2
        nnum(556) = 2
        nnum(557) = 2
        nnum(558) = 2
        nnum(559) = 1
         neigh(554,1) =  1
         neigh(554,2) = 555
         neigh(555,1) = 554
         neigh(555,2) = 556
         neigh(556,1) = 555
         neigh(556,2) = 557
         neigh(557,1) = 556
         neigh(557,2) = 558
         neigh(558,1) = 557
         neigh(558,2) = 559
         neigh(559,1) = 558

c connectivity of SD part
          nnum(1) = 2  ! SOMA
          neigh(1,1) = 554
          neigh(1,2) =  2

c Now proximal smooth dendrite: no spiny treelets attached
          nnum(2) = 2
          neigh(2,1) = 1
          neigh(2,2) = 3

          nnum(3) = 3
          neigh(3,1) = 2
          neigh(3,2) = 4
          neigh(3,3) = 5

c Now distal smooth dendrite (level 3): 2 spiny treelets to each
          nnum(4) = 6
          neigh(4,1) = 3
          neigh(4,2) = 5
          neigh(4,3) = 12
          neigh(4,4) = 18
          neigh(4,5) = 26 ! base of spiny treelet
          neigh(4,6) = 38 ! base of spiny treelet

          nnum(5) = 5
          neigh(5,1) = 3
          neigh(5,2) = 4
          neigh(5,3) = 6 
          neigh(5,4) = 50 ! base of spiny treelet
          neigh(5,5) = 62 ! base of spiny treelet

          nnum(6) = 4
          neigh(6,1) = 5
          neigh(6,2) = 7
          neigh(6,3) = 74
          neigh(6,4) = 86

          nnum(7) = 4
          neigh(7,1) = 6
          neigh(7,2) = 8
          neigh(7,3) = 98
          neigh(7,4) = 110

          nnum(8) = 4
          neigh(8,1) = 7
          neigh(8,2) = 9
          neigh(8,3) = 122
          neigh(8,4) = 134

          nnum(9) = 4
          neigh(9,1) = 8
          neigh(9,2) = 10
          neigh(9,3) = 146
          neigh(9,4) = 158

          nnum(10) = 4
          neigh(10,1) = 9
          neigh(10,2) = 11
          neigh(10,3) = 170
          neigh(10,4) = 182

          nnum(11) = 3
          neigh(11,1) = 10
          neigh(11,2) = 194
          neigh(11,3) = 206

          nnum(12) = 5
          neigh(12,1) = 4 
          neigh(12,2) = 18 
          neigh(12,3) = 13  
          neigh(12,4) = 218 
          neigh(12,5) = 230 

          nnum(13) = 4
          neigh(13,1) = 12
          neigh(13,2) = 14 
          neigh(13,3) = 242 
          neigh(13,4) = 254 

          nnum(14) = 4
          neigh(14,1) = 13
          neigh(14,2) = 15 
          neigh(14,3) = 266 
          neigh(14,4) = 278 

          nnum(15) = 4
          neigh(15,1) = 14
          neigh(15,2) = 16 
          neigh(15,3) = 290 
          neigh(15,4) = 302 

          nnum(16) = 4
          neigh(16,1) = 15
          neigh(16,2) = 17 
          neigh(16,3) = 314 
          neigh(16,4) = 326 

          nnum(17) = 3
          neigh(17,1) = 16
          neigh(17,2) = 338
          neigh(17,3) = 350 

          nnum(18) = 5
          neigh(18,1) = 4 
          neigh(18,2) = 12 
          neigh(18,3) = 19  
          neigh(18,4) = 362 
          neigh(18,5) = 374 

          nnum(19) = 4
          neigh(19,1) = 18
          neigh(19,2) = 20 
          neigh(19,3) = 386 
          neigh(19,4) = 398 

          nnum(20) = 4
          neigh(20,1) = 19
          neigh(20,2) = 21 
          neigh(20,3) = 410 
          neigh(20,4) = 422 

          nnum(21) = 4
          neigh(21,1) = 20
          neigh(21,2) = 22 
          neigh(21,3) = 434 
          neigh(21,4) = 446 

          nnum(22) = 4
          neigh(22,1) = 21
          neigh(22,2) = 23 
          neigh(22,3) = 458 
          neigh(22,4) = 470 

          nnum(23) = 4
          neigh(23,1) = 22
          neigh(23,2) = 24 
          neigh(23,3) = 482 
          neigh(23,4) = 494 

          nnum(24) = 4
          neigh(24,1) = 23
          neigh(24,2) = 25 
          neigh(24,3) = 506 
          neigh(24,4) = 518 

          nnum(25) = 3
          neigh(25,1) = 24
          neigh(25,2) = 530
          neigh(25,3) = 542 

c Attachments of bases of treelets: one side to a smooth compartment, the other to next more distal part of treelet
          do i = 26, 530, 24
           nnum(i) = 2
           neigh(i,1) = (i - 26)/24 + 4
           neigh(i,2) = i + 1
          end do

          do i = 38, 542, 24
           nnum(i) = 2
           neigh(i,1) = (i - 38)/24 + 4
           neigh(i,2) = i + 1
          end do

c Now connect up the rest of the treelets
           do i = 27, 543, 12
            nnum(i) = 2
            neigh(i,1) = i-1
            neigh(i,2) = i+1
           end do

           do i = 28, 544, 12
            nnum(i) = 2
            neigh(i,1) = i-1
            neigh(i,2) = i+1
           end do

           do i = 29, 545, 12
            nnum(i) = 3
            neigh(i,1) = i-1
            neigh(i,2) = i+1
            neigh(i,3) = i+5
           end do

           do i = 30, 546, 12
            nnum(i) = 3
            neigh(i,1) = i-1
            neigh(i,2) = i+1
            neigh(i,3) = i+4
           end do

           do i = 31, 547, 12
            nnum(i) = 2
            neigh(i,1) = i-1
            neigh(i,2) = i+1
           end do

           do i = 32, 548, 12
            nnum(i) = 2
            neigh(i,1) = i-1
            neigh(i,2) = i+1
           end do

           do i = 33, 549, 12
              nnum(i) = 1
              neigh(i,1) = i-1
           end do

           do i = 34, 550, 12
              nnum(i) = 3
              neigh(i,1) = i-5
              neigh(i,2) = i-4
              neigh(i,3) = i+1
           end do

           do i = 35, 551, 12
              nnum(i) = 2
              neigh(i,1) = i-1
              neigh(i,2) = i+1
           end do

           do i = 36, 552, 12
              nnum(i) = 2
              neigh(i,1) = i-1
              neigh(i,2) = i+1
           end do

           do i = 37, 553, 12
              nnum(i) = 1
              neigh(i,1) = i-1
           end do

         DO 332, I = 1, 559
           WRITE(6,3330) I, NEIGH(I,1),NEIGH(I,2),NEIGH(I,3),NEIGH(I,4),
     X NEIGH(I,5),NEIGH(I,6)
3330     FORMAT(2X,I5,I5,I5,I5,I5,I5,I5)
332      CONTINUE
          DO 858, I = 1,559
           DO 858, J = 1, NNUM(I)
            K = NEIGH(I,J)
            IT = 0
            DO 859, L = 1, NNUM(K)
             IF (NEIGH(K,L).EQ.I) IT = 1
859         CONTINUE
             IF (IT.EQ.0) THEN
              WRITE(6,8591) I, K
8591          FORMAT(' ASYMMETRY IN NEIGH MATRIX ',I4,I4)
             ENDIF
858       CONTINUE

c length and radius of axonal compartments
          do i = 554, 559
c           len(i) = 5.d0
            len(i) = 10.d0
          end do
          rad(554) = 0.750d0
          rad(555) = 0.700d0
          rad(556) = 0.650d0
          rad(557) = 0.600d0
          rad(558) = 0.550d0
          rad(559) = 0.5d0

c  length and radius of SD compartments
          len(1) = 29.d0
          rad(1) = 9.d0 ! SOMA

c Smooth dendrites
          do i = 2, 25
           len(i) = 15.d0
          end do
c Possibly lengthen dendritic shaft, to isolate (relatively) rest of dendrites.
          len(2) = 30.d0
          len(3) = 30.d0

          rad(2) =   2.25d0  ! proximal shaft
          rad(3) =   2.25d0  ! proximal shaft
          rad(2) =   1.80d0  ! proximal shaft ! What happens if this is narrower?
          rad(3) =   1.80d0  ! proximal shaft

           do i = 4, 11
             rad(i) = 1.8d0 ! a bit smaller
           end do
           do i = 12, 25
             rad(i) = 1.42d0
           end do

c Spiny dendrites
          do i = 26, 553
            len(i) = 25.d0
          end do

          do j = 26, 29
           do i = 0, 516, 12
             rad(i + j) = 0.75d0
           end do
          end do

          do j = 30, 37
           do i = 0, 516, 12
             rad(i + j) = 0.60d0
           end do
          end do


        WRITE(6,919)
919     FORMAT('COMPART.',' LEVEL ',' RADIUS ',' LENGTH(MU)')
        DO 920, I = 1, 559
920      WRITE(6,921) I, LEVEL(I), RAD(I), LEN(I)
921     FORMAT(I3,5X,I2,3X,F6.2,1X,F6.1,2X,F4.3)

        DO 120, I = 1, 559
          AREA(I) = 2.d0 * PI * RAD(I) * LEN(I)
          K = LEVEL(I)
         if (k .eq. 4) area(i) = 3.d0 * area(i)   ! spine correction
          C(I) = CDENS * AREA(I) * (1.D-8)

           if (k.ge.1) then
          GL(I) = (1.D-2) * AREA(I) / RM_SD
           else
          GL(I) = (1.D-2) * AREA(I) / RM_AXON
           endif

          GNAF(I) = GNAF_DENS(K) * AREA(I) * (1.D-5)
          GNAP(I) = GNAP_DENS(K) * AREA(I) * (1.D-5)
          GCAP(I) = GCAP_DENS(K) * AREA(I) * (1.D-5)
          GCAT(I) = GCAT_DENS(K) * AREA(I) * (1.D-5)
          GCAR(I) = GCAR_DENS(K) * AREA(I) * (1.D-5)
          GKDR(I) = GKDR_DENS(K) * AREA(I) * (1.D-5)
          GKA(I) = GKA_DENS(K) * AREA(I) * (1.D-5)
          GKD(I) = GKD_DENS(K) * AREA(I) * (1.D-5)
          GKC(I) = GKC_DENS(K) * AREA(I) * (1.D-5)
          GKAHP(I) = GKAHP_DENS(K) * AREA(I) * (1.D-5)
          GKM(I) = GKM_DENS(K) * AREA(I) * (1.D-5)
          GAR(I) = GAR_DENS(K) * AREA(I) * (1.D-5)
c above conductances should be in microS
120           continue

         Z = 0.d0
         DO I = 2, 25
           Z = Z + AREA(I)
1019     END DO    

         z1 = 0.d0
         do i = 26, 553
           z1 = z1 + area(i)
         end do
         WRITE(6,1020) area(1),Z, z1
1020     FORMAT(2X,' SOMA AREA ',F7.0,
     X 'smooth dend ',f7.0,' spiny dends c spines ',f7.0)

        DO 140, I = 1, 559
        DO 140, K = 1, NNUM(I)
         J = NEIGH(I,K)
           if (level(i).eq.0) then
               RI = RI_AXON
           else
               RI = RI_SD
           endif
         GAM1 =100.d0 * PI * RAD(I) * RAD(I) / ( RI * LEN(I) )

           if (level(j).eq.0) then
               RI = RI_AXON
           else
               RI = RI_SD
           endif
         GAM2 =100.d0 * PI * RAD(J) * RAD(J) / ( RI * LEN(J) )

         GAM(I,J) = 2.d0/( (1.d0/GAM1) + (1.d0/GAM2) )
140     CONTINUE
c gam computed in microS

        DO 299, I = 1, 559
c299       BETCHI(I) = .20d0 ! NOTE HOW FAST
299       BETCHI(I) = .80d0 ! NOTE HOW FAST
        BETCHI( 1) =  .05d0

        DO 300, I = 1, 559
c300     D(I) = 1.D-4
c300     D(I) = 3.D-2
300     D(I) = 6.D-2
        DO 301, I = 1, 559
         IF (LEVEL(I).EQ.1) D(I) = 3.D-2
301     CONTINUE

       DO 160, I = 1, 559
160     CAFOR(I) = 5200.d0 / (AREA(I) * D(I))
C     NOTE CORRECTION

        do 200, i = 1, 559
200     C(I) = 1000.d0 * C(I)
C     TO GO FROM MICROF TO NF.

      DO 909, I = 1, 559
       JACOB(I,I) = - GL(I)
      DO 909, J = 1, NNUM(I)
         K = NEIGH(I,J)
         IF (I.EQ.K) THEN
             WRITE(6,510) I
510          FORMAT(' UNEXPECTED SYMMETRY IN NEIGH ',I4)
         ENDIF
         JACOB(I,K) = GAM(I,K)
         JACOB(I,I) = JACOB(I,I) - GAM(I,K)
909   CONTINUE

c 15 Jan. 2001: make correction for c(i)
          do i = 1, 559
          do j = 1, 559
             jacob(i,j) = jacob(i,j) / c(i)
          end do
          end do

       DO 500, I = 1, 559
        WRITE (6,501) I,C(I)
501     FORMAT(1X,I3,' C(I) = ',F8.5)
500     CONTINUE
        END


       SUBROUTINE PURKINJE_INT (V,CHI,CINV,
     X mnaf,hnaf,mnap,mcap,mcat,hcat,mcar,hcar,mar,
     x mkdr,mkm,mka,hka,mkd,hkd,mkc,mkahp,
     x dt,neigh,nnum,jacob,mg,
     x vL,vk,vna,var,vca,vgaba_a,betchi,gam,gL,
     x gnaf,gnap,gcaP,gcaT,gCaR,gar,
     x gKDR,gKM,gKA,gKD,gKC,gKAHP,     
     x gampa,gnmda,ggaba_a,
     x O,time,
     X    alpham_naf, betam_naf, dalpham_naf, dbetam_naf,
     X    alphah_naf, betah_naf, dalphah_naf, dbetah_naf,
     X    alpham_nap, betam_nap, dalpham_nap, dbetam_nap,

     X    alpham_cap, betam_cap, dalpham_cap, dbetam_cap,
     X    alpham_cat, betam_cat, dalpham_cat, dbetam_cat,
     X    alphah_cat, betah_cat, dalphah_cat, dbetah_cat,
     X    alpham_car, betam_car, dalpham_car, dbetam_car,
     X    alphah_car, betah_car, dalphah_car, dbetah_car,

     X    alpham_ar , betam_ar , dalpham_ar , dbetam_ar,

     X    alpham_kdr, betam_kdr, dalpham_kdr, dbetam_kdr,
     X    alpham_km , betam_km , dalpham_km , dbetam_km ,
     X    alpham_ka , betam_ka , dalpham_ka , dbetam_ka ,
     X    alphah_ka , betah_ka , dalphah_ka , dbetah_ka ,
     X    alpham_kd , betam_kd , dalpham_kd , dbetam_kd ,
     X    alphah_kd , betah_kd , dalphah_kd , dbetah_kd ,
     X    alpham_kc , betam_kc , dalpham_kc , dbetam_kc ,
     X    cafor,curr)

c CINV is 1/C, i.e. inverse capacitance
       real*8 v(559), chi(559), cinv(559)
       real*8 mnaf(559), hnaf(559), mnap(559), mkdr(559),
     x mka(559),hka(559),mkm(559),mkc(559),mkahp(559),
     x mkd(559),hkd(559),mcap(559),mcar(559),hcar(559),
     x mcat(559),hcat(559),mar(559),jacob(559,559),betchi(559),
     x gam(0:559,0:559),gL(559),
     x gnaf(559),gnap(559),gkdr(559),gka(559),gkd(559),
     x gkm(559),gkc(559),gkahp(559),gcat(559),gar(559),
     x gcap(559),gcar(559),
     x gampa(559),gnmda(559),ggaba_a(559),cafor(559) 

       real*8
     X alpham_naf(0:640),betam_naf(0:640),dalpham_naf(0:640),
     X   dbetam_naf(0:640),
     X alphah_naf(0:640),betah_naf(0:640),dalphah_naf(0:640),
     X   dbetah_naf(0:640),
     X alpham_nap(0:640),betam_nap(0:640),dalpham_nap(0:640),
     X   dbetam_nap(0:640),

     X alpham_cap(0:640),betam_cap(0:640),dalpham_cap(0:640),
     X   dbetam_cap(0:640),
     X alpham_cat(0:640),betam_cat(0:640),dalpham_cat(0:640),
     X   dbetam_cat(0:640),
     X alphah_cat(0:640),betah_cat(0:640),dalphah_cat(0:640),
     X   dbetah_cat(0:640),
     X alpham_car(0:640),betam_car(0:640),dalpham_car(0:640),
     X   dbetam_car(0:640),
     X alphah_car(0:640),betah_car(0:640),dalphah_car(0:640),
     X   dbetah_car(0:640),

     X alpham_ar(0:640), betam_ar(0:640), dalpham_ar(0:640),
     X   dbetam_ar(0:640),

     X alpham_kdr(0:640),betam_kdr(0:640),dalpham_kdr(0:640),
     X   dbetam_kdr(0:640),
     X alpham_km(0:640), betam_km(0:640), dalpham_km(0:640),
     X   dbetam_km(0:640),
     X alpham_ka(0:640), betam_ka(0:640),dalpham_ka(0:640) ,
     X   dbetam_ka(0:640),
     X alphah_ka(0:640), betah_ka(0:640), dalphah_ka(0:640),
     X   dbetah_ka(0:640),
     X alpham_kd(0:640), betam_kd(0:640),dalpham_kd(0:640) ,
     X   dbetam_kd(0:640),
     X alphah_kd(0:640), betah_kd(0:640), dalphah_kd(0:640),
     X   dbetah_kd(0:640),
     X alpham_kc(0:640), betam_kc(0:640), dalpham_kc(0:640),
     X   dbetam_kc(0:640) 

        real*8 fastna_shift

c the f's are the functions giving 1st derivatives for evolution of
c the differential equations for the voltages (v), calcium (chi), and
c other state variables.
       real*8 fv(559), fchi(559),fmnaf(559),fhnaf(559),fmkdr(559),
     x fmka(559),fhka(559),fmkd(559),fhkd(559),fmnap(559),
     x fmkm(559),fmkc(559),fmkahp(559),
     x fmcat(559),fhcat(559),fmar(559),
     x fmcap(559),fmcar(559),fhcar(559)

c below are for calculating the partial derivatives
       real*8 dfv_dv(559,559), dfv_dchi(559), dfv_dmnaf(559),
     x  dfv_dhnaf(559),dfv_dmnap(559),
     x  dfv_dmkdr(559),dfv_dmka(559),dfv_dhka(559),
     x  dfv_dmkd(559),dfv_dhkd(559),dfv_dmkm(559),dfv_dmkc(559),
     xdfv_dmkahp(559),dfv_dmcat(559),dfv_dhcat(559),dfv_dmcap(559),
     x  dfv_dmar(559),dfv_dmcar(559),dfv_dhcar(559)

        real*8 dfchi_dv(559), dfchi_dchi(559),
     x dfmnaf_dmnaf(559), dfmnaf_dv(559),dfhnaf_dhnaf(559),
     x dfhnaf_dv(559),
     x dfmnap_dmnap(559), dfmnap_dv(559),
     x dfmkdr_dmkdr(559),dfmkdr_dv(559),
     x dfmka_dmka(559),dfmka_dv(559),dfhka_dhka(559),dfhka_dv(559),
     x dfmkd_dmkd(559),dfmkd_dv(559),dfhkd_dhkd(559),dfhkd_dv(559),
     x dfmkm_dmkm(559),dfmkm_dv(559),dfmkc_dmkc(559),dfmkc_dv(559),
     x dfmcap_dmcap(559),dfmcap_dv(559),
     x dfmcat_dmcat(559),dfmcat_dv(559),dfhcat_dhcat(559),
     x dfhcat_dv(559),
     x dfmcar_dmcar(559),dfmcar_dv(559),dfhcar_dhcar(559),
     x dfhcar_dv(559),
     x dfmar_dmar(559),dfmar_dv(559),dfmkahp_dchi(559),
     x dfmkahp_dmkahp(559), dt2, outrcd(20), time

         REAL*8 dt,mg,vL,vk,vna,var,vca,vgaba_a,curr(559),Z,Z1,Z2
         INTEGER O, K0, K1, K2, NEIGH(559,6), NNUM(559)
       REAL*8 OPEN(559),gamma(559),gamma_prime(559)
c gamma is function of chi used in calculating KC conductance
       REAL*8 alpham_ahp(559), alpham_ahp_prime(559)
       REAL*8 gna_tot(559),gk_tot(559),gca_tot(559),gar_tot(559)

       DO 301, I = 1, 559
          FV(I) = -GL(I) * (V(I) - VL) * cinv(i)
          DO 302, J = 1, NNUM(I)
             K = NEIGH(I,J)
302     FV(I) = FV(I) + GAM(I,K) * (V(K) - V(I)) * cinv(i)
301    CONTINUE


c       CALL FNMDA (V, OPEN, MG)

      DO 421, I = 1, 559
421    FV(I) = FV(I) + ( CURR(I)
     X   - gampa(I) * V(I) 
c    X   - (gampa(I) + open(i) * gnmda(I))*V(I)
     X   - ggaba_a(I)*(V(I)-Vgaba_a) ) * cinv(i)
c above assumes equil. potential for AMPA & NMDA = 0 mV

       do i = 1, 559
        gamma(i) = dmin1 (1.d0, .004d0 * chi(i))
        if (chi(i).le.250.d0) then
          gamma_prime(i) = .004d0
        else
          gamma_prime(i) = 0.d0
        endif
       end do

      DO 88, I = 1, 559
       gna_tot(i) = gnaf(i) * (mnaf(i)**3) * hnaf(i) +
     x     gnap(i) * (mnap(i)**3)
       gk_tot(i) = gkdr(i) * (mkdr(i)**4) +
     x             gka(i)  * (mka(i)**4) * hka(i) +
     x             gkd(i)  * (mkd(i)**4) * hkd(i) +
     x             gkm(i)  * mkm(i) +
     x             gkc(i)  * mkc(i) * gamma(i) +
     x             gkahp(i)* mkahp(i)
       gca_tot(i) = gcat(i) *  mcat(i) * hcat(i) + ! NOTE CHANGE TO 1st power m
     x              gcaP(i) *  mcaP(i) +
     x              gcaR(i) *  mcaR(i) * hcaR(i) 

       gar_tot(i) = gar(i) * mar(i)


88     FV(I) = FV(I) - ( gna_tot(i) * (v(i) - vna)
     X  + gk_tot(i) * (v(i) - vK)
     X  + gca_tot(i) * (v(i) - vCa)
     X  + gar_tot(i) * (v(i) - var) ) * cinv(i)

         do i = 1, 559
         do j = 1, 559
          if (i.ne.j) then
            dfv_dv(i,j) = jacob(i,j)
          else
            dfv_dv(i,j) = jacob(i,i) - cinv(i) *
     X  (gna_tot(i) + gk_tot(i) + gca_tot(i) + gar_tot(i)
     X    + ggaba_a(i) + gampa(i)
     X   + open(i) * gnmda(I) )
          endif
         end do
         end do

          do i = 1, 559
        dfv_dchi(i)  = - cinv(i) * gkc(i) * mkc(i) *
     x                     gamma_prime(i) * (v(i)-vK)
        dfv_dmnaf(i) = -3.d0 * cinv(i) * (mnaf(i)**2) *
     X    (gnaf(i) * hnaf(i) + gnap(i)) * (v(i) - vna)
        dfv_dhnaf(i) = - cinv(i) * gnaf(i) * (mnaf(i)**3) *
     X                    (v(i) - vna)
        dfv_dmnap(i) = -3.d0 * cinv(i) * (mnap(i)**2) *
     X    (gnap(i) + gnap(i)) * (v(i) - vna)

        dfv_dmkdr(i) = -4.d0 * cinv(i) * gkdr(i) * (mkdr(i)**3)
     X                   * (v(i) - vK)
        dfv_dmka(i)  = -4.d0 * cinv(i) * gka(i) * (mka(i)**3) *
     X                   hka(i) * (v(i) - vK)
        dfv_dhka(i)  = - cinv(i) * gka(i) * (mka(i)**4) *
     X                    (v(i) - vK)
        dfv_dmkd(i)  = -4.d0 * cinv(i) * gkd(i) * (mkd(i)**3) *
     X                   hkd(i) * (v(i) - vK)
        dfv_dhkd(i)  = - cinv(i) * gkd(i) * (mkd(i)**4) *
     X                    (v(i) - vK)

        dfv_dmkm(i)  = - cinv(i) * gkm(i) * (v(i) - vK)
        dfv_dmkc(i)  = - cinv(i) * gkc(i) * gamma(i) * (v(i)-vK)
        dfv_dmkahp(i)= - cinv(i) * gkahp(i) * (v(i) - vK)

        dfv_dmcat(i)  = -1.d0 * cinv(i) * gcat(i) * 1.d0    *
     X                    hcat(i) * (v(i) - vCa)
        dfv_dhcat(i) = - cinv(i) * gcat(i) * mcat(i) *
     X                  (v(i) - vCa)
        dfv_dmcar(i)  = -1.d0 * cinv(i) * gcar(i) * 1.d0    *
     X                    hcar(i) * (v(i) - vCa)
        dfv_dhcar(i) = - cinv(i) * gcar(i) * mcar(i) *
     X                  (v(i) - vCa)
        dfv_dmcap(i) = -1.d0 * cinv(i) * gcap(i) * 1.d0    *
     X                      (v(i) - vCa)

        dfv_dmar(i) = - cinv(i) * gar(i) * (v(i) - var)
          end do

         do i = 1, 559
          fchi(i) = - cafor(i) * gca_tot(i) * (v(i) - vca)
     x       - betchi(i) * chi(i)
          dfchi_dv(i) = - cafor(i) * gca_tot(i)
          dfchi_dchi(i) = - betchi(i)
         end do

       do i = 1, 559
c       alpham_ahp(i) = dmin1(0.2d-4 * chi(i),0.01d0) ! original from baskint
        alpham_ahp(i) = dmin1(6.0d-4 * chi(i),0.30d0) ! speed this up
        if (chi(i).le.500.d0) then
          alpham_ahp_prime(i) = 6.0d-4
        else
          alpham_ahp_prime(i) = 0.d0
        endif
       end do

       do i = 1, 559
        fmkahp(i) = alpham_ahp(i) * (1.d0 - mkahp(i))
c    x                  -.001d0 * mkahp(i)
     x                  -.060d0 * mkahp(i) ! speed this up significantly
        dfmkahp_dmkahp(i) = - alpham_ahp(i) - .001d0
        dfmkahp_dchi(i) = alpham_ahp_prime(i) *
     x                     (1.d0 - mkahp(i))
       end do

          do i = 1, 559

       K1 = IDNINT ( 4.d0 * (V(I) + 120.d0) ) ! For SD
       IF (K1.GT.640) K1 = 640
       IF (K1.LT.  0) K1 =   0

c            fastNa_shift =  2.0d0 ! For axon
             fastNa_shift =  6.0d0 ! For axon
c            fastNa_shift =  0.0d0 ! For axon
       K0 = IDNINT ( 4.d0 * (V(I)+      fastNa_shift+ 120.d0) )
       IF (K0.GT.640) K0 = 640
       IF (K0.LT.  0) K0 =   0


           if (i.lt.554) then
            k2 = k1
           else
            k2 = k0
           endif
        fmnaf(i) = alpham_naf(k2) * (1.d0 - mnaf(i)) -
     X              betam_naf(k2) * mnaf(i)
        fhnaf(i) = alphah_naf(k1) * (1.d0 - hnaf(i)) -
     X              betah_naf(k1) * hnaf(i)
        fmnap(i) = alpham_nap(k1) * (1.d0 - mnap(i)) -
     X              betam_nap(k1) * mnap(i)

        fmkdr(i) = alpham_kdr(k1) * (1.d0 - mkdr(i)) -
     X              betam_kdr(k1) * mkdr(i)
        fmka(i)  = alpham_ka (k1) * (1.d0 - mka(i)) -
     X              betam_ka (k1) * mka(i)
        fhka(i)  = alphah_ka (k1) * (1.d0 - hka(i)) -
     X              betah_ka (k1) * hka(i)
        fmkd(i)  = alpham_kd (k1) * (1.d0 - mkd(i)) -
     X              betam_kd (k1) * mkd(i)
        fhkd(i)  = alphah_kd (k1) * (1.d0 - hkd(i)) -
     X              betah_kd (k1) * hkd(i)
        fmkm(i)  = alpham_km (k1) * (1.d0 - mkm(i)) -
     X              betam_km (k1) * mkm(i)
        fmkc(i)  = alpham_kc (k1) * (1.d0 - mkc(i)) -
     X              betam_kc (k1) * mkc(i)

        fmcat(i) = alpham_cat(k1) * (1.d0 - mcat(i)) -
     X              betam_cat(k1) * mcat(i)
        fhcat(i) = alphah_cat(k1) * (1.d0 - hcat(i)) -
     X              betah_cat(k1) * hcat(i)
        fmcar(i) = alpham_car(k1) * (1.d0 - mcar(i)) -
     X              betam_car(k1) * mcar(i)
        fhcar(i) = alphah_car(k1) * (1.d0 - hcar(i)) -
     X              betah_car(k1) * hcar(i)
        fmcap(i) = alpham_cap(k1) * (1.d0 - mcap(i)) -
     X              betam_cap(k1) * mcap(i)

        fmar(i)  = alpham_ar (k1) * (1.d0 - mar(i)) -
     X              betam_ar (k1) * mar(i)

       dfmnaf_dv(i) = dalpham_naf(k2) * (1.d0 - mnaf(i)) -
     X                  dbetam_naf(k2) * mnaf(i)
       dfhnaf_dv(i) = dalphah_naf(k1) * (1.d0 - hnaf(i)) -
     X                  dbetah_naf(k1) * hnaf(i)
       dfmnap_dv(i) = dalpham_nap(k1) * (1.d0 - mnap(i)) -
     X                  dbetam_nap(k1) * mnap(i)

       dfmkdr_dv(i) = dalpham_kdr(k1) * (1.d0 - mkdr(i)) -
     X                  dbetam_kdr(k1) * mkdr(i)
       dfmka_dv(i)  = dalpham_ka(k1) * (1.d0 - mka(i)) -
     X                  dbetam_ka(k1) * mka(i)
       dfhka_dv(i)  = dalphah_ka(k1) * (1.d0 - hka(i)) -
     X                  dbetah_ka(k1) * hka(i)
       dfmkd_dv(i)  = dalpham_kd(k1) * (1.d0 - mkd(i)) -
     X                  dbetam_kd(k1) * mkd(i)
       dfhkd_dv(i)  = dalphah_kd(k1) * (1.d0 - hkd(i)) -
     X                  dbetah_kd(k1) * hkd(i)
       dfmkm_dv(i)  = dalpham_km(k1) * (1.d0 - mkm(i)) -
     X                  dbetam_km(k1) * mkm(i)
       dfmkc_dv(i)  = dalpham_kc(k1) * (1.d0 - mkc(i)) -
     X                  dbetam_kc(k1) * mkc(i)
       dfmcat_dv(i) = dalpham_cat(k1) * (1.d0 - mcat(i)) -
     X                  dbetam_cat(k1) * mcat(i)
       dfhcat_dv(i) = dalphah_cat(k1) * (1.d0 - hcat(i)) -
     X                  dbetah_cat(k1) * hcat(i)
       dfmcar_dv(i) = dalpham_car(k1) * (1.d0 - mcar(i)) -
     X                  dbetam_car(k1) * mcar(i)
       dfhcar_dv(i) = dalphah_car(k1) * (1.d0 - hcar(i)) -
     X                  dbetah_car(k1) * hcar(i)
       dfmcap_dv(i) = dalpham_cap(k1) * (1.d0 - mcap(i)) -
     X                  dbetam_cap(k1) * mcap(i)

       dfmar_dv(i)  = dalpham_ar(k1) * (1.d0 - mar(i)) -
     X                  dbetam_ar(k1) * mar(i)

       dfmnaf_dmnaf(i) =  - alpham_naf(k2) - betam_naf(k2)
       dfhnaf_dhnaf(i) =  - alphah_naf(k1) - betah_naf(k1)
       dfmnap_dmnap(i) =  - alpham_nap(k1) - betam_nap(k1)

       dfmkdr_dmkdr(i) =  - alpham_kdr(k1) - betam_kdr(k1)
       dfmka_dmka(i)  =   - alpham_ka (k1) - betam_ka (k1)
       dfhka_dhka(i)  =   - alphah_ka (k1) - betah_ka (k1)
       dfmkd_dmkd(i)  =   - alpham_kd (k1) - betam_kd (k1)
       dfhkd_dhkd(i)  =   - alphah_kd (k1) - betah_kd (k1)
       dfmkm_dmkm(i)  =   - alpham_km (k1) - betam_km (k1)
       dfmkc_dmkc(i)  =   - alpham_kc (k1) - betam_kc (k1)

       dfmcat_dmcat(i) =  - alpham_cat(k1) - betam_cat(k1)
       dfhcat_dhcat(i) =  - alphah_cat(k1) - betah_cat(k1)
       dfmcar_dmcar(i) =  - alpham_car(k1) - betam_car(k1)
       dfhcar_dhcar(i) =  - alphah_car(k1) - betah_car(k1)
       dfmcap_dmcap(i) =  - alpham_cap(k1) - betam_cap(k1)

       dfmar_dmar(i)  =   - alpham_ar (k1) - betam_ar (k1)

          end do

       dt2 = 0.5d0 * dt * dt

        do i = 1, 559
          v(i) = v(i) + dt * fv(i)
           do j = 1, 559
        v(i) = v(i) + dt2 * dfv_dv(i,j) * fv(j)
           end do
        v(i) = v(i) + dt2 * ( dfv_dchi(i) * fchi(i)
     X          + dfv_dmnaf(i) * fmnaf(i)
     X          + dfv_dhnaf(i) * fhnaf(i)
     X          + dfv_dmnap(i) * fmnap(i)
     X          + dfv_dmkdr(i) * fmkdr(i)
     X          + dfv_dmka(i)  * fmka(i)
     X          + dfv_dhka(i)  * fhka(i)
     X          + dfv_dmkd(i)  * fmkd(i)
     X          + dfv_dhkd(i)  * fhkd(i)
     X          + dfv_dmkm(i)  * fmkm(i)
     X          + dfv_dmkc(i)  * fmkc(i)
     X          + dfv_dmkahp(i)* fmkahp(i)
     X          + dfv_dmcat(i)  * fmcat(i)
     X          + dfv_dhcat(i) * fhcat(i)
     X          + dfv_dmcar(i)  * fmcar(i)
     X          + dfv_dhcar(i) * fhcar(i)
     X          + dfv_dmcap(i) * fmcap(i)
     X          + dfv_dmar(i)  * fmar(i) )

        chi(i) = chi(i) + dt * fchi(i) + dt2 *
     X   (dfchi_dchi(i) * fchi(i) + dfchi_dv(i) * fv(i))

        mnaf(i) = mnaf(i) + dt * fmnaf(i) + dt2 *
     X   (dfmnaf_dmnaf(i) * fmnaf(i) + dfmnaf_dv(i)*fv(i))
        hnaf(i) = hnaf(i) + dt * fhnaf(i) + dt2 *
     X   (dfhnaf_dhnaf(i) * fhnaf(i) + dfhnaf_dv(i)*fv(i))
        mnap(i) = mnap(i) + dt * fmnap(i) + dt2 *
     X   (dfmnap_dmnap(i) * fmnap(i) + dfmnap_dv(i)*fv(i))

        mkdr(i) = mkdr(i) + dt * fmkdr(i) + dt2 *
     X   (dfmkdr_dmkdr(i) * fmkdr(i) + dfmkdr_dv(i)*fv(i))
        mka(i) =  mka(i) + dt * fmka(i) + dt2 *
     X   (dfmka_dmka(i) * fmka(i) + dfmka_dv(i) * fv(i))
        hka(i) =  hka(i) + dt * fhka(i) + dt2 *
     X   (dfhka_dhka(i) * fhka(i) + dfhka_dv(i) * fv(i))
        mkd(i) =  mkd(i) + dt * fmkd(i) + dt2 *
     X   (dfmkd_dmkd(i) * fmkd(i) + dfmkd_dv(i) * fv(i))
        hkd(i) =  hkd(i) + dt * fhkd(i) + dt2 *
     X   (dfhkd_dhkd(i) * fhkd(i) + dfhkd_dv(i) * fv(i))
        mkm(i) =  mkm(i) + dt * fmkm(i) + dt2 *
     X   (dfmkm_dmkm(i) * fmkm(i) + dfmkm_dv(i) * fv(i))
        mkc(i) =  mkc(i) + dt * fmkc(i) + dt2 *
     X   (dfmkc_dmkc(i) * fmkc(i) + dfmkc_dv(i) * fv(i))
        mkahp(i) = mkahp(i) + dt * fmkahp(i) + dt2 *
     X (dfmkahp_dmkahp(i)*fmkahp(i) + dfmkahp_dchi(i)*fchi(i))
        mcat(i) =  mcat(i) + dt * fmcat(i) + dt2 *
     X   (dfmcat_dmcat(i) * fmcat(i) + dfmcat_dv(i) * fv(i))
        hcat(i) =  hcat(i) + dt * fhcat(i) + dt2 *
     X   (dfhcat_dhcat(i) * fhcat(i) + dfhcat_dv(i) * fv(i))
        mcar(i) =  mcar(i) + dt * fmcar(i) + dt2 *
     X   (dfmcar_dmcar(i) * fmcar(i) + dfmcar_dv(i) * fv(i))
        hcar(i) =  hcar(i) + dt * fhcar(i) + dt2 *
     X   (dfhcar_dhcar(i) * fhcar(i) + dfhcar_dv(i) * fv(i))
        mcap(i) =  mcap(i) + dt * fmcap(i) + dt2 *
     X   (dfmcap_dmcap(i) * fmcap(i) + dfmcap_dv(i) * fv(i))

        mar(i) =   mar(i) + dt * fmar(i) + dt2 *
     X   (dfmar_dmar(i) * fmar(i) + dfmar_dv(i) * fv(i))
         end do


c      IF (MOD(O,75).EQ.0) THEN
       IF (MOD(O,150).EQ.0) THEN
           OUTRCD(1) = TIME
           OUTRCD(2) = v(1)
           outrcd(3) = v(14)
           outrcd(4) = V(17)    
           outrcd(5) = v( 22)
           outrcd(6) = v(555)
           outrcd(7) = v(559)
           outrcd(8) = chi(1)
       outrcd(9) =  chi(14)
       outrcd(10)=  chi(22)
       outrcd(11)=  chi(17)
           outrcd(12) = curr(1)
           outrcd(13) = curr(559)
           outrcd(14) = curr( 14)
           outrcd(15) = gna_tot(1) * (v(1)-vna) * cinv(1)
           outrcd(16) = gk_tot(1) * (v(1)-vk) * cinv(1)
           outrcd(17) = gca_tot( 14) * (v( 14)-vca)*cinv( 14)
           outrcd(18) = gk_tot( 14) * (v( 14)-vk)*cinv( 14)
         OPEN(11,FILE='PURK17.OU')
         WRITE (11,FMT='(18F10.3)') (OUTRCD(I),I=1,18)
C900      FORMAT (100A4)
       ENDIF


              END

               SUBROUTINE FNMDA (VSTOR,OPEN,MG)
               REAL*8 VSTOR(559), OPEN(559)
               REAL*8 A, BB1, BB2, VM, A1, A2, B1, B2, MG
c modify so that potential is absolute and not relative to
c  "rest"
         A = DEXP(-2.847d0)
         BB1 = DEXP(-.693d0)
         BB2 = DEXP(-3.101d0)
C  TO DETERMINE VOLTAGE-DEPENDENCE OF NMDA CHANNELS
           DO 1, I = 1, 559
c          VM = VSTOR(I) - 60.
           VM = VSTOR(I)
           A1 = DEXP(-.016d0*VM - 2.91d0)
           A2 = 1000.d0 * MG * DEXP (-.045d0 * VM - 6.97d0)
           B1 = DEXP(.009d0*VM + 1.22d0)
           B2 = DEXP(.017d0*VM + 0.96d0)
        OPEN(I)     = 1.d0/(1.d0 + (A1+A2)*(A1*BB1 + A2*BB2) /
     X   (A*A1*(B1+BB1) + A*A2*(B2+BB2))  )
C  FROM JAHR & STEVENS, EQ. 4A
C               DO 124, J = 1, 19
C          OPEN(J) = 1./(1.+.667* EXP(-0.07*(VSTOR(J)-60.)) )
C  FROM CHUCK STEVENS
1               CONTINUE
                        END