Network model with dynamic ion concentrations (Ullah et al. 2009)

 Download zip file 
Help downloading and running models
Accession:144590
This is a network model composed of 100 excitatory and 100 inhibitory neurons with dynamic ion concentrations as described in "The Influence of Sodium and Potassium Dynamics on Excitability, Seizures, and the Stability of Persistent States: II. Network and Glia Dynamics (2009) Journal of Computational Neuroscience, 26:171-183".
Reference:
1 . Ullah G, Cressman JR, Barreto E, Schiff SJ (2009) The influence of sodium and potassium dynamics on excitability, seizures, and the stability of persistent states. II. Network and glial dynamics. J Comput Neurosci 26:171-83 [PubMed]
Citations  Citation Browser
Model Information (Click on a link to find other models with that property)
Model Type: Realistic Network; Extracellular;
Brain Region(s)/Organism:
Cell Type(s):
Channel(s):
Gap Junctions:
Receptor(s):
Gene(s):
Transmitter(s):
Simulation Environment: FORTRAN;
Model Concept(s): Epilepsy;
Implementer(s): Ullah, Ghanim [ghanim.phy at gmail.com];
/
Ullah_etal_2009
Readme.rtf
Readme.txt
Activity.f90
Network.f90
                            
!      Ghanim Ullah, John R Cressman Jr, Ernest Barreto, and Steven J Schiff,
!                        July 29, 2008,
!                        Penn State University,
!                        Email: ghanim.phy@gmail.com
!      Archived with 
!      "The Influence of Sodium and Potassium Dynamics on Excitability, Seizures, 
!      and the Stability of Persistent States: II. Network and Glia Dynamics" (2009) 
!      Journal of Computational Neuroscience, 26:171-183
!      This program couples 100 inhibitory neurons and 100 excitatory neurons where the membrane potential dynamics of 
!      these neurons is taken from Gutkin et al. model, 2001,
!      J. Computational Neuroscience, 11, 121-134. 
!      The synaptic currents here are modified from that given in Gutkin et al. model, 2001 
!      The model also includes dynamic potassium and sodium concentrations
!      that builds on the model from companion paper "The Influence of Sodium and Potassium Dynamics on Excitability, Seizures, 
!      and the Stability of Persistent States: I. Single Neuron Dynamics" by
!      John R Cressman Jr, Ghanim Ullah, Jokubas Ziburkus, Steven J Schiff , and Ernest Barreto (2009),
!      Journal of Computational Neuroscience, 26:159-170.
!      NOTE: It is important to start the simulation with steady state dynamics in order to reproduce 
!      various results for given parameters set. The single cell version 
!      of this code is also archived with Cressman et al., (2009) J. Comp. Neurosci. 26:159-170. 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 

       implicit none
       integer i, j
	   integer N, N1, dim                                   !number of cells from each cell type, dimensions of each cell, number of data points
       Parameter(N=100, dim=15, N1=10000000)
       double precision time, time_step                     !total time, time step for integration
	   double precision derivs(N*dim), x(N*dim), thresh     !thresh = spiking threshold
       double precision activity_ex, activity_in            !activity of excitatory cells, activity of inhibitory cells
       integer ee, bb      
       EXTERNAL rk4                                         !rk4 algorithm is used for integration            
       EXTERNAl SystemDerives                               !function for the dynamics of the network
       time_step = 0.01d0                                   !time step (in ms)
       thresh = 0.0                                         !spiking threshold
       time = 0.0                                           !total time (in ms)

       open(unit=1, file='raster_excitatory.dat')           !saves data for raster plot of excitatory network (huge file) 
       open(unit=2, file='raster_inhibitory.dat')           !saves data for raster plot of inhibitory network (huge file)
       open(unit=3, file='activity_excitatory.dat')         !saves data for activity of excitatory network
       open(unit=4, file='activity_inhibitory.dat')         !saves data for activity of inhibitory network
       open(unit=5, file='v_cell50_excitatory.dat')         !saves membrane potential of cell 50 in excitatory network 
       open(unit=6, file='ko_cell50_excitatory.dat')        !saves extracellular potassium for cell 50 in excitatory network
       open(unit=7, file='Nai_cell50_excitatory.dat')       !saves intracellular sodium for cell 50 in in excitatory network
	   open(unit=8,file='time.dat')                         !saves total time 

!      initial conditions
       do i = 1, N
          x((i-1)*dim+1) = -65.1034452264476897             !membrane potential of excitatory network
          x((i-1)*dim+2) = 0.645021617064742286d-01         !potassium channel activation variable n for excitatory network
          x((i-1)*dim+3) = 0.981306641820766101             !sodium channel inactivation variable h for excitatory network
          x((i-1)*dim+4) = 0.217441399671948571d-05         !variable s for temporal evolution of synaptic efficacy emanating from excitatory network 
          x((i-1)*dim+5) = 0.978772044450795857d-07         !calcium concentration for excitatory network
          x((i-1)*dim+6) = -65.1034452264476897             !membrane potential of inhibitory network
          x((i-1)*dim+7) = 0.645021617064742286d-01         !potassium channel activation variable n for excitatory network
          x((i-1)*dim+8) = 0.981306641820766101             !sodium channel inactivation variable h for inhibitory network
          x((i-1)*dim+9) = 0.217441399671948571d-05         !variable s for temporal evolution of synaptic efficacy emanating from inhibitory network   
          x((i-1)*dim+10) = 7.22                            !extracellular potassium concentration for excitatory network
          x((i-1)*dim+11) = 18.54                           !intracellular sodium concentration for excitatory network
          x((i-1)*dim+12) = 7.22                            !extracellular potassium concentration for inhibitory network
          x((i-1)*dim+13) = 18.54                           !intracellular sodium concentration for inhibitory network
          x((i-1)*dim+14) = 0.0								!variable eta, modeling the synaptic block due to the depolarization for excitatory network
          x((i-1)*dim+15) = 0.0								!variable eta, modeling the synaptic block due to the depolarization for inhibitory network  
       end do 
 
!     main program, calls the subroutines SystemDerives and rk4 to integrate and tracks the activity of the network
 
       do i = 1, N1/10                                      !begin time loop
          activity_ex = 0.0                                 !activity of the excitatory cells
          activity_in = 0.0                                 !activity of the inhibitory cells
          do j = 1, 10                                      !loops the rk4 method 10 times before recording a data point, can be changed
             time = time + time_step
             call SystemDerives(time,x,derivs)
             call rk4(x,derivs,N*dim,time,time_step,x,SystemDerives)
          end do
          do j = 1, N
             ee = 0
             bb = 0
             if(x((j-1)*dim+1) .gt. thresh)then
               ee = j
               write(1,*)time,ee
               activity_ex = activity_ex + 1.0
             end if
             if(x((j-1)*dim+6) .gt. thresh)then
               bb = j+N
               write(2,*)time,bb
               activity_in = activity_in + 1.0
             end if
          end do 
          write(3,*)activity_ex
          write(4,*)activity_in
		  write(5,*)x(49*dim+1)
		  write(6,*)x(49*dim+10)
		  write(7,*)x(49*dim+11)
		  write(8,*)time
       end do                                               !end time loop
       end
             
!      This subroutine gives the network dynamics

       SUBROUTINE SystemDerives(x,y,dydx)
       implicit none
       double precision PI, PIOVERTWO, TWOPI
       Parameter (PI = 3.14159265358979323846,PIOVERTWO = PI/2.0d0, TWOPI = 2.0d0 * PI)
       double precision ran2
       integer N, dim, ii, jj
       parameter(N=100, dim=15)

!	   Various parameters and variables used in the model are declared here

       double precision C, phi, tau_e, A, g_ca, V_ca, g_l, V_l, g_na, V_na_e,V_na_i, g_k, g_ahp,V_na
       double precision V_k, V_k_e, V_k_i, tau_i, V_ee, V_ie, V_ei, V_ii, alpha_ee, alpha_ie
       double precision alpha_ei, alpha_ii,alpha_g, g_ee, g_ie, g_ei, g_ii,g_g
       double precision dx1,dx2, specnum1, specnum2, synsum1, synsum2,gapsum
       double precision alpha_n_Ve, beta_n_Ve, alpha_n_Vi, beta_n_Vi
       double precision alpha_m_Ve, beta_m_Ve, alpha_m_Vi, beta_m_Vi
       double precision alpha_h_Ve, beta_h_Ve, alpha_h_Vi, beta_h_Vi
       double precision sigma_Ve, sigma_Vi, m_inf_Ve, m_inf_Vi
       double precision Ie_mem, Ie_syn, Ie_ext, Ie_rand,V_sp
       double precision Ii_mem, Ii_syn, Ii_ext, Ii_rand,vau
       double precision x, y(N*dim), dydx(N*dim)
       double precision thresh,see,sie,sei,sii
       double precision Ik_inf_e, tau_Ik_e, A_IK_e, iksmall_inf_e, p1_e, p2_e
       double precision gamma_Ik_e, alpha_Ik_e, c1, c2, c3, c4, min_e
       double precision beta, I_pump_e, I_diff_e, I_ki_e, I_ko_e, n_k0
       double precision Ik_inf_i, tau_Ik_i, A_IK_i, iksmall_inf_i, p1_i, p2_i
       double precision gamma_Ik_i, alpha_Ik_i, min_i
       double precision I_pump_i, I_diff_i, I_ki_i, I_ko_i, diffusion, deltax
       double precision diffusion_e, diffusion_i,epsilon,G_glia
       double precision gamma_e, gamma_i, gamma_telda, v_b, beta_telda_e, beta_telda_i
       integer seed1, seed2, seed3
	   double precision expotassium_current,exsodium_current,inpotassium_current,insodium_current
	   double precision sodium_oe, sodium_oi, current_conversion
	   double precision Kin_e, Kin_i
       double precision V_l_e, V_l_i

       beta = 7.0         
       n_k0 = 14.0				            	           !n_k0 is represented as k_o,infinity in the paper
       diffusion = 250.0			                       !diffusion coefficient D
       epsilon = 0.07/0.3                                  !rest of the variables as given in the paper
	   G_glia = 5.0/0.3
       deltax = 10.0       
	   v_b = -50.0
       gamma_telda = 0.4			                       
       vau=5
       C = 1.0
       phi = 3.0
       tau_e = 4.0
       A = 20.0
       g_ca = 0.1   
       V_ca = 120.0
       g_l = 0.05
       V_l = -65.0
       g_na = 100.0 
       V_na = 55.0
       g_k = 40.0
       g_ahp = 0.01
       V_k = -80.0
       tau_i = 8.0      
       V_ee = 0.0
       V_ie = -80.0
       V_ei = 0.0
       V_ii = -80.0
       V_sp=40
       specnum1 = sqrt(100.0/PI)
       specnum2 = sqrt(30.0/PI)
       thresh=0.0
       alpha_ee = 0.12 
	   alpha_ie = 0.06     
       alpha_ei = 0.2     
       alpha_ii = 0.02

       alpha_g=0.0					                      !coupling term for gap junctions NOT INCLUDED IN THE MODEL
	
       do  ii=1, N                                        !Loop for 100 excitatory and 100 inhibitory cells
      
		 seed1 = (-1)*int(ii*(x*100))
         seed2=(-1)*int((N-ii)*(x*100))

!        rates for various ions gatting variables, as in paper. The subscripts e and i respectively 
!        represent variables for excitatory and inhibitory network.

         alpha_n_Ve = 0.01 * (y(((ii-1)*dim)+1)+34.0)/( 1.0 - exp(-0.1 * (y(((ii-1)*dim)+1)+34.0)) )
         beta_n_Ve = 0.125 * exp(-1.0*(y(((ii-1)*dim)+1)+44.0)/80.0)
         alpha_n_Vi = 0.01 * (y(((ii-1)*dim)+6)+34.0)/( 1.0 - exp(-0.1 * (y(((ii-1)*dim)+6)+34.0)) )
         beta_n_Vi = 0.125 * exp(-1.0*(y(((ii-1)*dim)+6)+44.0)/80.0)
         alpha_m_Ve = 0.1 * (y(((ii-1)*dim)+1)+30.0)/( 1.0 - exp(-0.1 * (y(((ii-1)*dim)+1)+30.0)) )
         beta_m_Ve = 4.0 * exp(-1.0*(y(((ii-1)*dim)+1)+55.0)/18.0)
         alpha_m_Vi = 0.1 * (y(((ii-1)*dim)+6)+30.0)/( 1.0 - exp(-0.1 * (y(((ii-1)*dim)+6)+30.0)) )
         beta_m_Vi = 4.0 * exp(-1.0*(y(((ii-1)*dim)+6)+55.0)/18.0)
         alpha_h_Ve = 0.07 * exp(-1.0*(y(((ii-1)*dim)+1)+44.0)/20.0)
         beta_h_Ve = 1.0/( 1.0 + exp(-0.1 * (y(((ii-1)*dim)+1)+14.0)) )
         alpha_h_Vi = 0.07 * exp(-1.0*(y(((ii-1)*dim)+6)+44.0)/20.0)
         beta_h_Vi = 1.0/( 1.0 + exp(-0.1 * (y(((ii-1)*dim)+6)+14.0)) )
         sigma_Ve = 1.0/( 1.0 + exp(-1.0*(y(((ii-1)*dim)+1)+20.0)/4.0) )
         sigma_Vi = 1.0/( 1.0 + exp(-1.0*(y(((ii-1)*dim)+6)+20.0)/4.0) )
         m_inf_Ve = alpha_m_Ve/(alpha_m_Ve + beta_m_Ve)
         m_inf_Vi = alpha_m_Vi/(alpha_m_Vi + beta_m_Vi)

!        Potassium and Sodium dynamics, I_diff and I_glia from the paper are combined into I_diff here

         I_pump_e = (1.25/(1.0+exp((25.0-y(((ii-1)*dim)+11))/3.0)))*(1.0/(1.0+exp(8.0-y(((ii-1)*dim)+10)))) 
         I_diff_e = epsilon*(y(((ii-1)*dim)+10)-n_k0)+G_glia/(1.0d0 + exp((18.0-y(((ii-1)*dim)+10))/2.5d0))
         I_pump_i = (1.25/(1.0+exp((25.0-y(((ii-1)*dim)+13))/3.0)))*(1.0/(1.0+exp(8.0-y(((ii-1)*dim)+12)))) 
         I_diff_i = epsilon*(y(((ii-1)*dim)+12)-n_k0)+G_glia/(1.0d0 + exp((18.0-y(((ii-1)*dim)+12))/2.5d0))
		 Kin_e=140.0+(18.0-y(((ii-1)*dim)+11))    
		 Kin_i=140.0+(18.0-y(((ii-1)*dim)+13))         
         sodium_oe = 144.0 - beta * (y(((ii-1)*dim)+11) - 18.0)  
		 sodium_oi = 144.0 - beta * (y(((ii-1)*dim)+13) - 18.0)

!        reversal potentials are updated based on instantaneous ion concentrations

         V_k_e = 26.64 * log(y(((ii-1)*dim)+10)/Kin_e) 
         V_k_i = 26.64 * log(y(((ii-1)*dim)+12)/Kin_i) 
         V_na_e = 26.64 * log(sodium_oe/y(((ii-1)*dim)+11)) 
         V_na_i = 26.64 * log(sodium_oi/y(((ii-1)*dim)+13)) 
         V_l_e = 26.64 * log ( (y(((ii-1)*dim)+10) + 0.065*sodium_oe + 0.6*6.0) / (Kin_e + 0.065*y(((ii-1)*dim)+11) + 0.6*130.0))
         V_l_i = 26.64 * log ( (y(((ii-1)*dim)+12) + 0.065*sodium_oi + 0.6*6.0) / (Kin_i + 0.065*y(((ii-1)*dim)+13) + 0.6*130.0))

!        diffusion of potassium from cell's extracellular space to the nearest neighbours         

         if(ii .eq. 1)then
            diffusion_e = diffusion*(y(((ii)*dim)+10)+y(((ii-1)*dim)+12)-2.0*y(((ii-1)*dim)+10))/(deltax*deltax)
            diffusion_i = diffusion*(y(((ii)*dim)+12)+y(((ii-1)*dim)+10)-2.0*y(((ii-1)*dim)+12))/(deltax*deltax)
         else if(ii .eq. N)then
            diffusion_e = diffusion*(y(((ii-2)*dim)+10)+y(((ii-1)*dim)+12)-2.0*y(((ii-1)*dim)+10))/(deltax*deltax)
            diffusion_i = diffusion*(y(((ii-2)*dim)+12)+y(((ii-1)*dim)+10)-2.0*y(((ii-1)*dim)+12))/(deltax*deltax)
         else
            diffusion_e = diffusion*(y(((ii-2)*dim)+10)+y(((ii)*dim)+10)+y(((ii-1)*dim)+12)-3.0*y(((ii-1)*dim)+10))/(deltax*deltax)
            diffusion_i = diffusion*(y(((ii-2)*dim)+12)+y(((ii)*dim)+12)+y(((ii-1)*dim)+10)-3.0*y(((ii-1)*dim)+12))/(deltax*deltax)
         end if

!        membrane currents for excitatory cells

         Ie_mem = -g_l * (y(((ii-1)*dim)+1) - V_l_e)            &
		 &       - g_na * (m_inf_Ve * m_inf_Ve * m_inf_Ve) * y(((ii-1)*dim)+3) * (y(((ii-1)*dim)+1)-V_na_e)        &
		 &       - (g_k * (y(((ii-1)*dim)+2) * y(((ii-1)*dim)+2) * y(((ii-1)*dim)+2) * y(((ii-1)*dim)+2))          &
		 &               + g_ahp * y(((ii-1)*dim)+5)/(1.0+y(((ii-1)*dim)+5)) )*(y(((ii-1)*dim)+1)-V_k_e)         
		 expotassium_current = (g_k * (y(((ii-1)*dim)+2) * y(((ii-1)*dim)+2) * y(((ii-1)*dim)+2) * y(((ii-1)*dim)+2))    &
		 &               + g_ahp * y(((ii-1)*dim)+5)/(1.0+y(((ii-1)*dim)+5)) )*(y(((ii-1)*dim)+1)-V_k_e)    
		 exsodium_current = g_na * (m_inf_Ve * m_inf_Ve * m_inf_Ve) * y(((ii-1)*dim)+3) * (y(((ii-1)*dim)+1)-V_na_e)   
		 
!******* synaptic current for excitatory cells; if uncoupled, use Ie_syn = 0.0; ***************/

         thresh=0
		 synsum1 = 0.0
         synsum2 = 0.0
         do jj=1, N
!                if((x .gt. 400.0).and.(x .lt. 430.0))THEN                   !use this loop for figure 4
!                    alpha_ee = 0.25                           
!                else
!                    alpha_ee = 0.2165
!                end if
    			 see=0
	      		 sie=0 
			     dx1=float(ii-jj)/float(N)
				 dx2=float(ii+(N-jj))/float(N)
                 g_ee = alpha_ee * specnum1 * (exp(-100.0*(dx1*dx1))+exp(-100.0*(dx2*dx2)))  
                 g_ie = alpha_ie * specnum2 * (exp(-30.0*(dx1*dx1))+exp(-30.0*(dx2*dx2))) 
			     if(ii .eq. jj)g_ee=0
			     if(y(((jj-1)*dim)+1) .gt. thresh)see=1
				 if(y(((jj-1)*dim)+6) .gt. thresh)sie=1
				 if(y(((jj-1)*dim)+14) .gt. 5.0)THEN
				   beta_telda_e = y(((jj-1)*dim)+14)
				 else
				   beta_telda_e = 0.0
				 end if  
				 if(y(((jj-1)*dim)+15) .gt. 5.0)THEN
				   beta_telda_i = y(((jj-1)*dim)+15)
				 else
				   beta_telda_i = 0.0
				 end if  
                 synsum1 = synsum1 + g_ee * y(((jj-1)*dim)+4) * exp(-beta_telda_e/vau)
                 synsum2 = synsum2 + g_ie * y(((jj-1)*dim)+9) * exp(-beta_telda_i/vau)
		 end do		 
         Ie_syn = - (y(((ii-1)*dim)+1) - V_ee) * synsum1 / float(N)         &
         &        - (y(((ii-1)*dim)+1) - V_ie) * synsum2 / float(N)

!!!**********END synaptic currents for excitatory cells**************!!!!!!!!

!		 if ( ((x .gt. 12.0) .and. (x .lt. 32.0)).and.((ii .ge. 21).and.(ii .le. 79)))then                  !use this loop for figure 2,3, 4
!             Ie_ext = 1.5*exp(-60.0*( (float(ii)-float(N)/2.0)/float(N) )*( (float(ii)-float(N)/2.0)/float(N) ) )
!        else if(((x .gt. 400.0).and.(x .le. 420.0)).and.((ii .ge. 1).and.(ii .le. 100)))then
!             Ie_ext = 0.85d0        !1.49,.99        !2nd stim 0.5d0
!        else
		      Ie_ext = 0.0
!		 end if

         Ie_rand = 20*(0.5-ran2(seed1))

!         /* currents into the inhibitory neurons */

         Ii_mem = - g_l * (y(((ii-1)*dim)+6) - V_l_i)           &
         &        - g_na * (m_inf_Vi * m_inf_Vi * m_inf_Vi) * y(((ii-1)*dim)+8) * (y(((ii-1)*dim)+6)-V_na_i)       &
         &        - g_k * (y(((ii-1)*dim)+7) * y(((ii-1)*dim)+7) * y(((ii-1)*dim)+7) * y(((ii-1)*dim)+7))          &
         &          * (y(((ii-1)*dim)+6)-V_k_i)
	     inpotassium_current = g_k * (y(((ii-1)*dim)+7) * y(((ii-1)*dim)+7) * y(((ii-1)*dim)+7) * y(((ii-1)*dim)+7))    &
         &                     * (y(((ii-1)*dim)+6)-V_k_i)  
		 insodium_current = g_na * (m_inf_Vi * m_inf_Vi * m_inf_Vi) * y(((ii-1)*dim)+8) * (y(((ii-1)*dim)+6)-V_na_i)  
         
!******** synaptic current into inhibitory neurons; if uncoupled, use Ii_syn = 0.0; *********************/

         thresh=0
		 synsum1 = 0.0
         synsum2 = 0.0
		 gapsum=0.0
		 do jj = 1, N
			     sei=0
			     sii=0
			     dx1=float(ii-jj)/float(N)
				 dx2=float(ii+(N-jj))/float(N)
				 
				 if(y(((jj-1)*dim)+1) .gt. thresh) sei=1
				 if(y(((jj-1)*dim)+6) .gt. thresh) sii=1
                 g_ei = alpha_ei * specnum2 * (exp( -30.0*(dx1*dx1) )+exp( -30.0*(dx2*dx2) )) 
                 g_ii = alpha_ii * specnum2 * (exp( -30.0*(dx1*dx1) )+exp( -30.0*(dx2*dx2) ))
!				 g_g=alpha_g*specnum2 * (exp( -30.0*(dx1*dx1) )+exp( -30.0*(dx2*dx2) ))
				 if(jj .eq. ii)g_ii=0
				 if(y(((jj-1)*dim)+14) .gt. 5.0)THEN
				   beta_telda_e = y(((jj-1)*dim)+14)
				 else
				   beta_telda_e = 0.0
				 end if  
				 if(y(((jj-1)*dim)+15) .gt. 5.0)THEN
				   beta_telda_i = y(((jj-1)*dim)+15)
				 else
				   beta_telda_i = 0.0
				 end if  
                 synsum1 = synsum1 + g_ei * y(((jj-1)*dim)+4) * exp(-beta_telda_e/vau)
                 synsum2 = synsum2 + g_ii * y(((jj-1)*dim)+9) * exp(-beta_telda_i/vau)
!                gapsum=gapsum+g_g*(y(((jj-1)*dim)+6)-y(((ii-1)*dim)+6))
		 end do
         Ii_syn = - (y(((ii-1)*dim)+6) - V_ei) * synsum1 / float(N)       &
         &        - (y(((ii-1)*dim)+6) - V_ii) * synsum2 / float(N)         ! &
!		 &		 +gapsum/float(N)

! *********END synaptic currents into inhibitory cells*****************!!!!!!!!!!!

		 Ii_ext = 0.5                  
		 Ii_rand =20.0*(0.5-ran2(seed2)) 

!!!!!!!!!variables modeling the interplay between neurons (see paper)
        
		 gamma_e = 0.0
         gamma_i = 0.0
         if((y(((ii-1)*dim)+1) .gt. -30.0).and.(y(((ii-1)*dim)+1) .lt. -10.0))THEN
            gamma_e = 0.4
         end if   
         if((y(((ii-1)*dim)+6) .gt. -30.0).and.(y(((ii-1)*dim)+6) .lt. -10.0))THEN
            gamma_i = 0.4
         end if         
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!         

!         /* excitatory neurons */

         dydx(((ii-1)*dim)+1) = (1/C) * (Ie_mem + Ie_ext + Ie_rand + Ie_syn)
         dydx(((ii-1)*dim)+2) = phi * ( alpha_n_Ve * (1.0-y(((ii-1)*dim)+2)) - beta_n_Ve * y(((ii-1)*dim)+2) )
         dydx(((ii-1)*dim)+3) = phi * ( alpha_h_Ve * (1.0-y(((ii-1)*dim)+3)) - beta_h_Ve * y(((ii-1)*dim)+3) )
         dydx(((ii-1)*dim)+4) = (1/tau_e) * (A * sigma_Ve * (1.0 - y(((ii-1)*dim)+4)) - y(((ii-1)*dim)+4))
         dydx(((ii-1)*dim)+5) = -0.002 * g_ca * (y(((ii-1)*dim)+1) - V_ca)/( 1.0 + exp( -1.0*(y(((ii-1)*dim)+1)+25.0)/2.5 ) ) - y(((ii-1)*dim)+5)/80.0

!         /* inhibitory neurons */

         dydx(((ii-1)*dim)+6) = (1/C) * (Ii_mem + Ii_ext + Ii_rand + Ii_syn)
         dydx(((ii-1)*dim)+7) = phi * ( alpha_n_Vi * (1.0-y(((ii-1)*dim)+7)) - beta_n_Vi * y(((ii-1)*dim)+7) )
         dydx(((ii-1)*dim)+8) = phi * ( alpha_h_Vi * (1.0-y(((ii-1)*dim)+8)) - beta_h_Vi * y(((ii-1)*dim)+8) )
         dydx(((ii-1)*dim)+9) = (1/tau_i) * (A * sigma_Vi * (1.0 - y(((ii-1)*dim)+9)) - y(((ii-1)*dim)+9))

!        Extracellular Potassium, Intracellular Sodium equations
!        Excitatory Neuron

         dydx(((ii-1)*dim)+10) = 0.001*(expotassium_current/3.0 - 7.0* 2.0*I_pump_e - I_diff_e + diffusion_e)        !factor 0.001 converts from seconds to milliseconds  
         dydx(((ii-1)*dim)+11) = 0.001*(-1.0*exsodium_current/(7.0*3.0) - 3.0*I_pump_e)  

!        Extracellular Potassium, Intracellular Sodium equations         
!        Inhibitory Neuron         

         dydx(((ii-1)*dim)+12) = 0.001*(inpotassium_current/3.0 - 7.0* 2.0*I_pump_i - I_diff_i + diffusion_i)
         dydx(((ii-1)*dim)+13) = 0.001*(-1.0*insodium_current/(7.0*3.0) - 3.0*I_pump_i)

!        The depolarization variable
!        Excitatory Neuron

         dydx(((ii-1)*dim)+14) = gamma_e * (y(((ii-1)*dim)+1) - v_b) - gamma_telda * y(((ii-1)*dim)+14)

!        The depolarization variable         
!        Inhibitory Neuron

         dydx(((ii-1)*dim)+15) = gamma_i * (y(((ii-1)*dim)+6) - v_b) - gamma_telda * y(((ii-1)*dim)+15)

       end do                             !end LOOP for 100 excitatory and 100 inhibitory cells
       end

       SUBROUTINE rk4(y,dydx,n,x,h,yout,derivs) 
       INTEGER n
       double precision h,x,dydx(n),y(n),yout(n) 
       EXTERNAL derivs 
       INTEGER i 
       double precision h6,hh,xh,dym(n),dyt(n),yt(n) 
       hh=h*0.5 
       h6=h/6. 
       xh=x+hh

       do i = 1, n     ! First step. 
          yt(i)=y(i)+hh*dydx(i) 
       end do  
       call derivs(xh,yt,dyt)    !  Secondstep. 
       do i = 1, n 
          yt(i)=y(i)+hh*dyt(i) 
       end do  
       call derivs(xh,yt,dym)    !  Thirdstep. 
       do i = 1, n 
          yt(i)=y(i)+h*dym(i) 
          dym(i)=dyt(i)+dym(i) 
       end do  
       call derivs(x+h,yt,dyt)   !  Fourthstep. 
       do i = 1, n          !  Accumulate increments with proper weights. 
          yout(i)=y(i)+h6*(dydx(i)+dyt(i)+2.*dym(i)) 
       end do 
       END 

       FUNCTION ran2(idum)
       INTEGER idum,IM1,IM2,IMM1,IA1,IA2,IQ1,IQ2,IR1,IR2,NTAB,NDIV
       REAL ran2,AM,EPS,RNMX
       PARAMETER (IM1=2147483563,IM2=2147483399,AM=1./IM1,      & 
       &     IMM1=IM1-1,IA1=40014,IA2=40692,IQ1=53668,          &
       &     IQ2=52774,IR1=12211,IR2=3791,NTAB=32,              &
       &     NDIV=1+IMM1/NTAB,EPS=1.2e-7,RNMX=1.-EPS)
       INTEGER idum2,j,k,iv(NTAB),iy
       SAVE iv,iy,idum2
       DATA idum2/123456789/, iv/NTAB*0/, iy/0/
          if (idum.le.0)then
             idum=max(-idum,1) 
             idum2=idum
             do  j=NTAB+8,1,-1 
                k=idum/IQ1
                idum=IA1*(idum-k*IQ1)-k*IR1
                   if (idum.lt.0) idum=idum+IM1
                     if (j .le. NTAB) iv(j)=idum
             end do 
             iy=iv(1)
          end if
          k=idum/IQ1
          idum=IA1*(idum-k*IQ1)-k*IR1 
          if (idum.lt.0) idum=idum+IM1
          k=idum2/IQ2
          idum2=IA2*(idum2-k*IQ2)-k*IR2
          if (idum2.lt.0) idum2=idum2+IM2
          j=1+iy/NDIV
          iy=iv(j)-idum2 
          iv(j)=idum
          if(iy.lt.1)iy=iy+IMM1
          ran2=min(AM*iy,RNMX) 
       END