Cortical Interneuron & Pyramidal Cell Model of Cortical Spreading Depression (Stein & Harris 2022)

 Download zip file 
Help downloading and running models
Accession:267033
This 2-cell cortical circuit model consists of a negative feedback loop between a single compartment pyramidal cell and a single compartment interneuron. Ion concentrations in the extra- and intracellular spaces are included in the model. The model is used to test the contribution of cortical inhibitory interneurons to the initiation of cortical spreading depression, as characterized by spike block in the pyramidal cell. Results show that interneuronal inhibition provides a wider dynamic range to the circuit and generally improves stability against spike block. Despite these beneficial effects, strong interneuronal firing contributed to rapidly changing extracellular ion concentrations, which facilitated hyperexcitation and led to spike block first in the interneuron and then in the pyramidal cell. The model results demonstrate that while the role of interneurons in cortical microcircuits is complex, they are critical to the initiation of pyramidal cell spike block and CSD. See reference below for more details.
Reference:
1 . Stein W, Harris AL (2022) Interneuronal dynamics facilitate the initiation of spike block in cortical microcircuits Journal of Computational Neuroscience [PubMed]
Model Information (Click on a link to find other models with that property)
Model Type: Neuron or other electrically excitable cell; Realistic Network;
Brain Region(s)/Organism: Neocortex;
Cell Type(s):
Channel(s):
Gap Junctions:
Receptor(s): GabaA; Glutamate;
Gene(s):
Transmitter(s): Gaba; Glutamate;
Simulation Environment: FORTRAN;
Model Concept(s): Spreading depression;
Implementer(s):
Search NeuronDB for information about:  GabaA; Glutamate; Gaba; Glutamate;
/
2pi
readme.txt
2pi.model.f
2pi.in
const.f
ode.par
                            
! This program is based on the model by A. L. Harris and W. Stein
! See details in BioRXIV 
!	https://biorxiv.org/cgi/content/short/2021.04.25.441350v1
! It is a 2 cell model with connected pyramidal neuron and interneuron
!  through gaba and glutamate synapses
! The pyramidal cell and interneuron reversal potentials vary
!  with ion concentrations
! Ion concentrations are updated in real-time

! This code needs the following files: ode.par, 2pi.in, const.f

!Follow these steps to compile the code
!1) compile const.f module
! 	gfortran -fdefault-real-8 -c const.f 
!2) compile 2pi.model.f code
!	mpif90 -fdefault-real-8 -c 2pi.model.f
!3) link module and code
!	mpif90 const.o 2pi.model.o -o 2pi.exe
!4) code can be executed using MPI fortran. Max processors is (Jeend - Jestart)/Jestep + 1
!	mpirun -np xx 2pi.exe
!  (xx = number of processors)

!The code produces the following output files:
!1) frequency.out.int.Jex.x.Jiy.y.gabaz.z
!	contains the firing frequency (in Hz) of the interneuron as a function of 
!	pyramidal cell injected current (Je), interneuron injected current (Ji),
! 	gaba maximum conductance (ggaba), and time
!	Files are labeled as x.x = Je value, y.y = Ji value, z.z=ggaba value
!2) frequency.out.pyr.Jex.x.Jiy.y.gabaz.z
!	same as (1), but for the interneuron
!3) if writeyn=1, membrane.pot.out.Jex.x.Jiy.y.gabaz.z
!	contains the pyramidal cell and interneuron membrane potential (mV)
!  	as a function of time - these are large files



	program twocell
	implicit none

	include 'ode.par'

!MPI	
	include 'mpif.h'
	integer num_proc,my_rank,ierr
	integer status(MPI_STATUS_SIZE)
	integer begNJe,finNJe,aveNJe,extra
!end MPI


! Declare variables
	real ystart(nvar)
	real ak2(NMAX),ak3(NMAX),ak4(NMAX),ak5(NMAX),ak6(NMAX),xp,yp
	real dxsav,hdid,hnext,htry,x,dydx,y,yscal,errmax,htemp,xnew
	real yerr,ytemp,SAFETY,PGROW,PSHRINK,ERRCON,dnm(6)
	integer Nt,Nx,Ncl,kmax,Nstep
	integer kount,nok,nbad,n
	character*20 real_to_char
	real ti,tf,eps,h1,dxsav_den
	external derivs
	real timei,timef,runtime,gamim,ggaba
	real Jestart,Jeend,Jestep,Je
	integer cntJe,NJe,cntJi_int,NJi_int,i,cntggaba,Nggaba
	character*3 dum,dumJi,dumgaba
	real Ji_intstart,Ji_intend,Ji_intstep,Ji_int
	real dum1,dum2,dum3,dum4
	real ggabastart,ggabaend,ggabastep
	integer writeyn
	real V0,n0,h0,Ca0,K_i0,K_o0
	real Na_i0,Na_o0,Cl_i0,Cl_o0,lils0
	real lilse0,lilsi0,V0_int,n0_int,h0_int
	real K_i0_int,Na_i0_int


	call cpu_time (timei)

!MPI
	call MPI_INIT(ierr)
	call MPI_COMM_SIZE(MPI_COMM_WORLD,num_proc,ierr)
	call MPI_COMM_RANK(MPI_COMM_WORLD,my_rank,ierr)

!end MPI


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Read input file								
										
	call readinput(ti,tf,ystart,Nstep,gamim,ggabastart
     &	,ggabaend,ggabastep,Jestart,Jeend,Jestep,Ji_intstart,
     &	Ji_intend,Ji_intstep,writeyn,n0,h0,Ca0,K_i0,K_o0,
     &	Na_i0,Na_o0,Cl_i0,Cl_o0,lils0,lilse0,lilsi0,V0_int,
     &	n0_int,h0_int,V0,K_i0_int,Na_i0_int)				
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


! setup loop over pyramidal cell injected current Je
	NJe = nint((Jeend - Jestart)/Jestep) + 1

	aveNJe = NJe/num_proc
	extra = mod(NJe,num_proc)

	begNJe = my_rank*aveNJe+1
	if(num_proc.eq.1)then
	  finNJe = (my_rank+1)*aveNJe
	else
	  finNJe = (my_rank+1)*aveNJe+(my_rank/(num_proc-1))*extra
	endif

	
! loop over Je values - this loop is parallel	
	do cntJe = begNJe,finNJe
	  Je = Jestart + (cntJe-1)*Jestep

	  write(dum,"(F3.1)")Je

! loop over interneuron injecte current values Ji 
	  NJi_int = nint((Ji_intend - Ji_intstart)/Ji_intstep) + 1
	  do cntJi_int = 1,NJi_int
	    Ji_int = Ji_intstart + (cntJi_int -1)*Ji_intstep

	    write(dumji,"(F3.1)")Ji_int

! loop over gaba conductance ggaba
	    Nggaba = nint((ggabaend - ggabastart)/ggabastep) + 1
	    do cntggaba = 1,Nggaba
	      ggaba = ggabastart + (cntggaba - 1)*ggabastep

	      write(dumgaba,"(F3.1)")ggaba

! open/create membrane potential files
	      if(writeyn.eq.1)then
	        open(1000*cntJi_int+100*cntggaba+cntJe,file=
     &'membrane.pot.out.Je'//dum//'.Ji'//dumji//'.gaba'//dumgaba)
	        write(1000*cntJi_int+100*cntggaba+cntJe,*)           
     &'time (s)     V_pyr (mV)    V_int (mV)'
	      endif

! open/create frequency file
	      write(dum,"(F3.1)")Je
	      open(10000*cntJi_int+1000*cntggaba+10*cntJe,file=
     &'frequency.out.pyr.Je'//dum//'.Ji'//dumji//'.gaba'//dumgaba)
	      write(10000*cntJi_int+1000*cntggaba+10*cntJe,*)
     &		'Je     Ji     ggaba     time (s)      freq (Hz)'

	      open(100000*cntJi_int+10000*cntggaba+100*cntJe,file=
     &'frequency.out.int.Je'//dum//'.Ji'//dumji//'.gaba'//dumgaba)
	      write(100000*cntJi_int+10000*cntggaba+100*cntJe,*)
     &		'Je     Ji     ggaba     time (s)      freq_p (Hz)'


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Solve ODEs with Runge-Kutta 	 	
!  See Numerical recipes												

	      call rkdumb(ystart,ti,tf,nstep,derivs,gamim,ggaba,Je
     &	    ,cntJe,Ji_int,cntJi_int,cntggaba,writeyn)										

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

! if writing to membrane potential files
	      if(writeyn.eq.1)then
! write input parameters to output file
		write(1000*cntJi_int+100*cntggaba+cntJe,*)
     &		'#V0 = ',V0
		write(1000*cntJi_int+100*cntggaba+cntJe,*)
     &		'#n0,h0,Ca0 = ',n0,h0,Ca0
		write(1000*cntJi_int+100*cntggaba+cntJe,*)
     &		'#K_i0,K_o0 = ',K_i0,K_o0
		write(1000*cntJi_int+100*cntggaba+cntJe,*)
     &		'#Na_i0,Na_o0 = ',Na_i0,Na_o0
		write(1000*cntJi_int+100*cntggaba+cntJe,*)
     &		'#Cl_i0,Cl_o0 = ',Cl_i0,Cl_o0
		write(1000*cntJi_int+100*cntggaba+cntJe,*)
     &		'#lils0,lilse0,lilsi0 = ',
     &		lils0,lilse0,lilsi0
		write(1000*cntJi_int+100*cntggaba+cntJe,*)
     &		'#V0_int = ',V0_int
		write(1000*cntJi_int+100*cntggaba+cntJe,*)
     &		'#n0_int,h0_int, = ',
     &		n0_int,h0_int
		write(1000*cntJi_int+100*cntggaba+cntJe,*)
     &		'#K_i0_int,Na_i0_int = ',K_i0_int,Na_i0_int
		write(1000*cntJi_int+100*cntggaba+cntJe,*)
     &		'#Nstep = ',Nstep
		write(1000*cntJi_int+100*cntggaba+cntJe,*)
     &		'#gami = ',gamim 

	        close(1000*cntJi_int+100*cntggaba+cntJe)
	      endif

! write input parameters to output file pyramidal
		write(10000*cntJi_int+1000*cntggaba+10*cntJe,*)
     &		'#V0 = ',V0
		write(10000*cntJi_int+1000*cntggaba+10*cntJe,*)
     &		'#n0,h0,Ca0 = ',n0,h0,Ca0
		write(10000*cntJi_int+1000*cntggaba+10*cntJe,*)
     &		'#K_i0,K_o0 = ',K_i0,K_o0
		write(10000*cntJi_int+1000*cntggaba+10*cntJe,*)
     &		'#Na_i0,Na_o0 = ',Na_i0,Na_o0
		write(10000*cntJi_int+1000*cntggaba+10*cntJe,*)
     &		'#Cl_i0,Cl_o0 = ',Cl_i0,Cl_o0
		write(10000*cntJi_int+1000*cntggaba+10*cntJe,*)
     &		'#lils0,lilse0,lilsi0 = ',
     &		lils0,lilse0,lilsi0
		write(10000*cntJi_int+1000*cntggaba+10*cntJe,*)
     &		'#V0_int = ',V0_int
		write(10000*cntJi_int+1000*cntggaba+10*cntJe,*)
     &		'#n0_int,h0_int = ',
     &		n0_int,h0_int
		write(10000*cntJi_int+1000*cntggaba+10*cntJe,*)
     &		'#K_i0_int,Na_i0_int = ',K_i0_int,Na_i0_int
		write(10000*cntJi_int+1000*cntggaba+10*cntJe,*)
     &		'#Nstep = ',Nstep
		write(10000*cntJi_int+1000*cntggaba+10*cntJe,*)
     &		'#gami = ',gamim 

	      close(10000*cntJi_int+1000*cntggaba+10*cntJe)

! write input parameters to output file interneuron
		write(100000*cntJi_int+10000*cntggaba+100*cntJe,*)
     &		'#V0 = ',V0
		write(100000*cntJi_int+10000*cntggaba+100*cntJe,*)
     &		'#n0,h0,Ca0 = ',n0,h0,Ca0
		write(100000*cntJi_int+10000*cntggaba+100*cntJe,*)
     &		'#K_i0,K_o0 = ',K_i0,K_o0
		write(100000*cntJi_int+10000*cntggaba+100*cntJe,*)
     &		'#Na_i0,Na_o0 = ',Na_i0,Na_o0
		write(100000*cntJi_int+10000*cntggaba+100*cntJe,*)
     &		'#Cl_i0,Cl_o0 = ',Cl_i0,Cl_o0
		write(100000*cntJi_int+10000*cntggaba+100*cntJe,*)
     &		'#lils0,lilse0,lilsi0 = ',
     &		lils0,lilse0,lilsi0
		write(100000*cntJi_int+10000*cntggaba+100*cntJe,*)
     &		'#V0_int = ',V0_int
		write(100000*cntJi_int+10000*cntggaba+100*cntJe,*)
     &		'#n0_int,h0_int ',
     &		n0_int,h0_int
		write(100000*cntJi_int+10000*cntggaba+100*cntJe,*)
     &		'#Nstep = ',Nstep
		write(100000*cntJi_int+10000*cntggaba+100*cntJe,*)
     &		'#gami = ',gamim 

	      close(100000*cntJi_int+10000*cntggaba+100*cntJe)

	    enddo	!ggaba loop

	  enddo		!Ji_int loop

	enddo	!Je loop

	call cpu_time (timef)
      	runtime=(timef-timei)
	print*,'total runtime was ',runtime,' seconds'

! MPI
	call MPI_FINALIZE(ierr)
!end MPI

	stop
	end program twocell



! Subroutines (mostly) in order that they are called

! This routine contains the differential equations to be integrated
! modified from Numerical Recipes for the model here
*=====================================================================*
	subroutine derivs(x,y,dydx,gamim,ggaba,Je,Ji_int)
*	These are the derivatives of the original functions, which are
*	ystart(i)
*=====================================================================*

	use const
	implicit none

	include 'ode.par'

	dimension y(nmax),dydx(nmax)	
	real x,y,dydx
	real V,n,h,Ca,m
	real IL,Ik,INa,INap,IAHP
	real an,bn,ah,bh,am,bm
	real ma1,minf
	real EL,Ek,ENa,ECl
	real c_k_i,c_k_o,c_Na_i,c_Na_o,c_Cl_i,c_Cl_o
	real IkL,Ipump
	real IKcc,INKcc,IKi,Isink
	real IClL,Igaba,lils
	real INaL
	real se,si,Iglute,Igluti

	real V_int,n_int,h_int,Ca_int
	real IL_int,Ik_int,INa_int,IAHP_int
	real an_int,bn_int,ah_int,bh_int
	real ma1_int
	real ECa_int,ECl_int
	real m_int,INaL_int,IKL_int,IClL_int,minf_int
	real am_int,bm_int
	real gamim,gami,ggaba,Je,Ji_int
	real Ek_int,ENa_int
	real c_k_i_int,c_Na_i_int,Ipump_int
	real ICa

	gami = gamim*gam

!=======USER SUPPLIED==================================
! these are the coupled DE's for the pyramidal neuron
!  and the concentrations

	V = y(1)
	n = y(2)
	h = y(3)
	Ca = y(4)

! concentrations are labeled with c_xxx_i or c_xxx_o
! the i/o refers to inside or outside, xxx is the ion

	c_k_i = y(5)
	c_k_o = y(6)
	c_Na_i = y(7)
	c_Na_o = y(8)
	c_Cl_i = y(9)
	c_Cl_o = y(10)
	lils = y(11)
	se = y(12)
	si = y(13)


	Ek = RTF*log(c_k_o/c_k_i)
	ENa = RTF*log(c_Na_o/c_Na_i)
	ECl = -RTF*log(c_Cl_o/c_Cl_i)

	dydx(11) = -lils/taugaba
	Igaba = ggaba*lils*(V - ECl)

	dydx(12) = -se/tauglut
	dydx(13) = -si/tauglut

	IkL = gkL*(V-Ek)
	IClL = gClL*(V - ECl)
	INaL = gNaL*(V - ENa)

	IL = IkL + IClL + INaL

	an = 0.032*(V+52.)/(1. - exp(-(V+52.)/5.))
	bn = 0.5*exp(-(V+57.)/40.)

	ah = 0.128*exp(-(V+50.)/18.)
	bh = 4./(1. + exp(-(V+27.)/5.))
	
	am = 0.32*(V+ 54.)/(1. - exp(-(V+ 54.)/4.))
	bm = .28*(V+27)/(exp((V+ 27.)/5.)-1.)

	ma1 = 1./(1. + exp(-(V+25.)/2.5))

	minf = am/(am+bm)
	m = minf

	Ik = gk*n**4.*(V - Ek)
	INa = gNa*m**3.*h*(V - ENa)
	INap = gp*m**3.*(V - ENa)
	IAHP = gAHP*(Ca/(Ca + 1.))*(V - Ek)
	ICa = gCa*ma1*(V - ECa)

	Ipump = (rhopump/gam)/( (1.+exp((Nasat-c_Na_i)/3.))
     &	  *(1.+exp(Ksat - c_K_o)) )
	IKcc = rhokcc*log(c_K_i*c_Cl_i/(c_K_o*c_Cl_o))
	INKcc = rhoNKcc*( log(c_K_i*c_Cl_i/(c_K_o*c_Cl_o))
     &	  + log(c_Na_i*c_Cl_i/(c_Na_o*c_Cl_o)) )
     &	  /(1.+exp(16. - c_K_o))
	Isink = epsK*(c_K_o - Kbath)
	Iglute = gglut1*se*(V - Eglut)


	dydx(5) = -tauinv*(gam*(Ik + IAHP + IKL - 2.*Ipump) 
     &	  + IKcc + INKcc)
	dydx(7) = tauinv*( -gam*(INa + INap + INaL + 3.*Ipump) 
     &	  - INKcc )
	dydx(9) = tauinv*( gam*(Igaba + IClL) - IKcc - 2.*INKcc )
	dydx(10) = -beta*dydx(9)


	dydx(1) = (Je - IL - Ik - INa - INap - Igaba
     &	  - Iglute - IAHP - Ipump)/C
	dydx(2) = phi*(an*(1.-n) - bn*n)
	dydx(3) = phi*(ah*(1.-h) - bh*h)
	dydx(4) = -epsi*gCa*ma1*(V - ECa) - Ca/tauCa

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! these are the coupled DE's for the interneuron and concentrations
	V_int = y(14)
	n_int = y(15)
	h_int = y(16)
	c_k_i_int = y(17)
	c_Na_i_int = y(18)

	Ek_int = RTF*log(c_k_o/c_k_i_int)
	ENa_int = RTF*log(c_Na_o/c_Na_i_int)

! USER - if you want to write the ion concentrations as a function of time
!	uncomment the next line - the files are large

!	write(27,*)x/1000.,c_k_i,c_Na_i,c_k_o,c_Na_o,c_k_i_int,c_Na_i_int

	am_int = 0.1*(V_int + 35.)/(1. - exp(-(V_int + 35.)/10.))
	bm_int = 4.*exp(-(V_int + 60.)/18.)
	ah_int = 0.07*exp(-(V_int + 58.)/20.)
	bh_int = 1./(1. + exp(-(V_int + 28.)/10.))
	an_int = 0.01*(V_int + 34.)/(1. - exp(-(V_int + 34.)/10.))
	bn_int = 0.125*(exp(-(V_int + 44.)/80.))

	minf_int = am_int/(am_int + bm_int)
	m_int = minf_int

	IkL_int = gkL_int*(V_int - Ek_int)
	INaL_int = gNaL_int*(V_int - ENa_int)
	IL_int = IkL_int + INaL_int
	
	Ipump_int = (rhopump/gami)/( (1.+exp((Nasat-c_Na_i_int)/3.))
     &	  *(1.+exp(Ksat - c_K_o)) )
	Ik_int = gk_int*n_int**4.*(V_int - Ek_int)
	INa_int = gNa_int*m_int**3.*h_int*(V_int - ENa_int)

	Igluti = gglut2*si*(V_int - Eglut)

	dydx(14) = (Ji_int - IL_int - Ik_int - INa_int 
     &	  - Igluti - Ipump_int)/C_int

	dydx(15) = phi_int*(an_int*(1.-n_int) - bn_int*n_int)
	dydx(16) = phi_int*(ah_int*(1.-h_int) - bh_int*h_int)

	dydx(17) = -tauinv*(gami*(Ik_int + IKL_int - 2.*Ipump_int) )
	dydx(18) = tauinv*(-gami*(INa_int + INaL_int + 3.*Ipump_int))

	dydx(6) = tauinv*(gam*beta*(Ik + IAHP + IkL - 2.*Ipump)
     &	  + beta*(IKcc + INKcc) - Isink) - beta*dydx(17) 		
	dydx(8) = -beta*(dydx(7) + dydx(18))


!=======END USER SUPPLIED===============================


	return
	end



!-----------------------------------------------------------------------
! This function converts a real number into a character, so that it can be
! used for naming files, etc.

	function real_to_char(angle)
!-----------------------------------------------------------------------
		
	real angle
	character*20 cangle
	character*10 frmt
	character*20 real_to_char
	

	if(angle.gt.9999)then
	   print*,'angle is too big for function real_to_char'
	   stop	   
	else if(angle.lt.-9999)then
	   print*,'angle is too small for function real_to_char'
	   stop	
	else if(angle.eq.0)then
	   frmt='(F7.5)'
	else if(abs(angle).lt.10)then
	   frmt='(F7.5)'
	else
	   frmt='(F7.5)'
	endif

	write(cangle,frmt)angle

	real_to_char=cangle

	end


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!Reads input file
	subroutine readinput(ti,tf,ystart,Nstep,gamim,ggabastart,
     &	ggabaend,ggabastep,Jestart,Jeend,Jestep,Ji_intstart,
     &	Ji_intend,Ji_intstep,writeyn,n0,h0,Ca0,K_i0,K_o0,
     &	Na_i0,Na_o0,Cl_i0,Cl_o0,lils0,lilse0,lilsi0,V0_int,
     &	n0_int,h0_int,V0,K_i0_int,Na_i0_int)

	implicit none
	include 'ode.par'

	real ti,tf,V0,n0,h0,Ca0,K_i0,K_o0
	real Na_i0,Na_o0,Cl_i0,Cl_o0,lils0
	integer Nstep,writeyn
	real ystart(nvar)
	real lilse0,lilsi0,V0_int,n0_int,h0_int
	real gamim,Jestart,Jeend,Jestep
	real Ji_intstart,Ji_intend,Ji_intstep
	real ggabastart,ggabaend,ggabastep
	real K_i0_int,Na_i0_int

	open(2,file='2pi.in')
	
	read(2,*)ti
	read(2,*)tf
	read(2,*)V0
	read(2,*)n0,h0,Ca0
	read(2,*)K_i0,K_o0
	read(2,*)Na_i0,Na_o0
	read(2,*)Cl_i0,Cl_o0
	read(2,*)lils0,lilse0,lilsi0
	read(2,*)V0_int
	read(2,*)n0_int,h0_int
	read(2,*)K_i0_int
	read(2,*)Na_i0_int
	read(2,*)Nstep
	read(2,*)gamim
	read(2,*)ggabastart,ggabaend,ggabastep
	read(2,*)Jestart,Jeend,Jestep
	read(2,*)Ji_intstart,Ji_intend,Ji_intstep
	read(2,*)writeyn

	if(nstep.gt.nstpmx)then
	  print*,'nstep > nstpmax'
	  print*,'increase nstpmax'
	  stop
	endif

! put initial conditions into array ystart
	ystart(1) = V0
	ystart(2) = n0
	ystart(3) = h0
	ystart(4) = Ca0

	ystart(5) = K_i0
	ystart(6) = K_o0
	ystart(7) = Na_i0
	ystart(8) = Na_o0
	ystart(9) = Cl_i0
	ystart(10) = Cl_o0

	ystart(11) = lils0
	ystart(12) = lilse0	
	ystart(13) = lilsi0

	ystart(14) = V0_int
	ystart(15) = n0_int
	ystart(16) = h0_int

	ystart(17) = K_i0_int
	ystart(18) = Na_i0_int

	return
	end

! From Numerical Recipes
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	SUBROUTINE rk4(y,dydx,n,x,h,yout,derivs,gamim,ggaba,Je,
     &	  Ji_int)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

	INTEGER n,NMAX
	REAL h,x,dydx(n),y(n),yout(n)
	EXTERNAL derivs
	PARAMETER (NMAX=50) 	!Set to the maximum number of functions.
!Given values for the variables y(1:n) and their derivatives dydx(1:n) known at x, use
!the fourth-order Runge-Kutta method to advance the solution over an interval h and return
!the incremented variables as yout(1:n), which need not be a distinct array from y. The
!user supplies the subroutine derivs(x,y,dydx), which returns derivatives dydx at x.
	INTEGER i
	REAL h6,hh,xh,dym(NMAX),dyt(NMAX),yt(NMAX)
	real gamim,ggaba,Je,Ji_int

	hh=h*0.5
	h6=h/6.
	xh=x+hh

	do i=1,n 	!First step.
	  yt(i)=y(i)+hh*dydx(i)
	enddo 

	call derivs(xh,yt,dyt,gamim,ggaba,Je,Ji_int) 	!Second step.

	do i=1,n
	  yt(i)=y(i)+hh*dyt(i)
	enddo 

	call derivs(xh,yt,dym,gamim,ggaba,Je,Ji_int) 	!Third step.

	do i=1,n
	  yt(i)=y(i)+h*dym(i)
	  dym(i)=dyt(i)+dym(i)
	enddo 

	call derivs(x+h,yt,dyt,gamim,ggaba,Je,Ji_int) 	!Fourth step.

	do i=1,n 	!Accumulate increments with proper weights.
	  yout(i)=y(i)+h6*(dydx(i)+dyt(i)+2.*dym(i))
	enddo 

	return
	END

! From Numerical Recipes
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	SUBROUTINE rkdumb(vstart,x1,x2,nstep,derivs,gamim,ggaba,Je,
     &	cntJe,Ji_int,cntJi_int,cntggaba,writeyn)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

	include 'ode.par'

	INTEGER nstep
!Maximum number of functions and
!maximum number of values to be stored.
	REAL x1,x2,vstart(nvar),xx(NSTPMX),y(NMAX,NSTPMX)
	EXTERNAL derivs
!	COMMON /path/ xx,y 	!Storage of results.
! USES rk4
!Starting from initial values vstart(1:nvar) known at x1 use fourth-order Runge-Kutta to
!advance nstep equal increments to x2. The user-supplied subroutine derivs(x,v,dvdx)
!evaluates derivatives. Results are stored in the common block path. Be sure to dimension
!the common block appropriately.
	INTEGER i,k,j,cntJe,cntJi_int,cntggaba,writeyn,j_int
	REAL h,x,dv(NMAX),v(NMAX)
	real dVp(NSTPMX),dVpp,freq,tspike(NSTPMX),dVp_int(NSTPMX)
	real gamim,ggaba,Je,Ji_int,dVpp_int,tspike_int(NSTPMX)
	real freq_int
	integer printme

	do i=1,nvar 	!Load starting values.
	  v(i)=vstart(i)
	  y(i,1)=v(i)
	enddo 

	xx(1)=x1
	x=x1
	h=(x2-x1)/nstep	


! write the functions to the output file

	j = 1
	printme = 1

	do k=1,nstep 	!Take nstep steps.

	  if(v(14) .gt. 0.)then		!v(14) is interneuron V
	    v(11) = 1.
	  endif

	  if(v(1) .gt. 0.)then		!v(1) is pyramidal neuron V
	    v(12) = 1.
	    v(13) = 1.
	  endif

	  call derivs(x,v,dv,gamim,ggaba,Je,Ji_int)
	  call rk4(v,dv,nvar,x,h,v,derivs,gamim,ggaba,Je,Ji_int)

	  if(x+h.eq.x)then
	    print*,'stepsize not significant in rkdumb'
	    stop
	  endif

	  x=x+h
	  xx(k+1)=x 	!Store intermediate steps.

	  do i=1,nvar
	    y(i,k+1)=v(i)
	  enddo 

!write membrane potentials to file
	  if(writeyn.eq.1)then
	    write(1000*cntJi_int+100*cntggaba+cntJe,'(6(f6.2,3x))')
     &	x/1000.,y(1,k),y(14,k)		!time,V_pyramid,V_int
	  endif

!calculate frequency as a function of time and write to file
! 1st and 2nd derivs of V for pyramidal neuron
	  if(k.gt.2)then
	    dVp(k) = (y(1,k) - y(1,k-1))/h
	    dVpp = (y(1,k) - 2.*y(1,k-1) + y(1,k-2))/h**2.

	    dVp_int(k) = (y(14,k) - y(14,k-1))/h
	    dVpp_int = (y(14,k) - 2.*y(14,k-1) + y(14,k-2))/h**2.
	  endif

	  if(k.gt.3)then
	    if(dVp(k).lt.0. .and. dVp(k-1).gt.0. .and. dVpp.lt.0. 	!1st deriv changes from + to -
     &	    .and. y(1,k).gt.0.)then					!and 2nd deriv<0 and V_pyramid>0
	      tspike(j) = x/1000.
	
	      if(j.gt.1)then
		freq = 1./(tspike(j) - tspike(j-1))
		write(10000*cntJi_int+1000*cntggaba+10*cntJe,'(5(f6.2,3x))')
     &	Je,Ji_int,ggaba,tspike(j-1),freq

	      endif		!j.gt.1

	      j = j+1
	    endif		!if a spike happens

! interneuron frequency
	    if(dVp_int(k).lt.0. .and. dVp_int(k-1).gt.0.  	!1st deriv changes from + to -
     &	    .and. dVpp_int.lt.0. .and. y(14,k).gt.0.)then	!and 2nd deriv<0 and V_int>0
	      tspike_int(j_int) = x/1000.
	
	      if(j_int.gt.1)then
		freq_int = 1./(tspike_int(j_int) - tspike_int(j_int-1))
		write(100000*cntJi_int+10000*cntggaba+100*cntJe,'(5(f6.2,3x))')
     &	Je,Ji_int,ggaba,tspike_int(j_int-1),freq_int

	      endif	!j_int.gt.1

	      j_int = j_int+1
	    endif	!if a spike happens

	  endif		!k.gt.3	      

	enddo 


	return
	END



Loading data, please wait...