Not sure if this is the right place to ask this but I've been attempting to implement the Pinsky-Rinzel model in the Brian2 simulator. I'm basing my implementation off the original article here: http://www.caam.rice.edu/~cox/neuro/pinskyrinzel.pdf and I've borrowed some of the formalism from the XPP version of the model here: http://senselab.med.yale.edu/ModelDB/showmodel.cshtml?model=35358&file=/b04feb12/booth_bose.ode (Though they are shifted 60mV from where I'm working)
I've got it to replicate the somatic firing (as you can see if you set gc to zero), but I can't figure out how to get the complex dendritc and somatic interplay I'm supposed to see as seen in figure 2 and figure 3 of the model and my q and Ca values don't replicate the model at all. I think there is something wrong in how I'm handling the calcium currents, or the units I'm using there.
Does anyone have any ideas of where I'm going wrong?
My code:
# -*- coding: utf-8 -*-
"""
Created on Sat Oct 24 11:21:30 2015
@author: Conor
"""
from __future__ import division
import numpy as np
from brian2 import *
defaultclock.dt = 0.05*ms
gbarl=0.1*msiemens
gbarna=30.0*msiemens
gbarkDr=15.0*msiemens
gbarCa=10.0*msiemens
gbarkc=15.0*msiemens
gkahp=0.8*msiemens
gNMDA=0*msiemens
gampa=0*msiemens
Vna=120.0*mV
Vca=140.0*mV
Vk=-15.0*mV
Vl=0*mV
Vsyn=60*mV
Id=0*uamp
Is=.15*uamp
gc=2.1*msiemens
p=0.5
Cm=3*uF
eqs='''
alpham=0.32*(13.1-Vs/mV)/(exp((13.1-Vs/mV)/4)-1)/ms : Hz
betam=0.28*(Vs/mV-40.1)/(exp((Vs/mV-40.1)/5)-1)/ms : Hz
alphan=0.016*(35.1-Vs/mV)/(exp((35.1-Vs/mV)/5)-1)/ms : Hz
betan=0.25*exp(0.5-0.025*Vs/mV)/ms : Hz
alphah=0.128*exp((17-Vs/mV)/18)/ms : Hz
betah=4/(1+exp((40-Vs/mV)/5))/ms : Hz
alphas=1.6/(1+exp(-0.072*(Vd/mV-65)))/ms : Hz
betas=.02*(Vd/mV-51.1)/(exp((Vd/mV-51.1)/5)-1)/ms : Hz
alphac=(exp((Vd/mV-10)/11)-exp((Vd/mV-6.5)/27))/(18.975)* (Vd<=50*mV)/ms+2*exp((6.5-Vd/mV)/27)*(Vd>50*mV)/ms : Hz
betac=((2*exp((6.5-Vd/mV)/27))/ms-alphac)*(Vd<(50*mV)) +0/ms : Hz
alphaq=clip(.00002*Ca,0,0.01)/ms :Hz
betaq=0.001/ms : Hz
Xca=clip(Ca/250,0,1) : 1
dh/dt= (alphah*(1-h)-betah*h) : 1
dn/dt = (alphan*(1-n)-betan*n) : 1
ds/dt = alphas*(1-s)-betas*s : 1
dc/dt = alphac*(1-c)-betac*c : 1
dq/dt= (alphaq*(1-q)-betaq*q) : 1
dCa/dt=-0.13*Ica/uamp/ms-0.075*Ca/ms : 1
Minfs=alpham/(alpham+betam) : 1
ILeakVs=gbarl*(Vs-Vl) : amp
ILeakVd=gbarl*(Vd-Vl) : amp
INa=gbarna*(Minfs**2)*h*(Vs-Vna) : amp
Ik_dr=gbarkDr*n*(Vs-Vk) : amp
Ica=gbarCa*s*s*(Vd-Vca) : amp
IK_c=gbarkc*Xca*c*(Vd-Vk) : amp
IK_AHP=gkahp*q*(Vd-Vk) : amp
dVs/dt=(-ILeakVs-INa-Ik_dr+gc/p*(Vd-Vs)+Is/p)/Cm : volt
dVd/dt=(-ILeakVd-Ica-IK_AHP-IK_c+gc/(1-p)*(Vs-Vd)+INMDA/(1-p))/Cm : volt
INMDA=gNMDA*Si*(1+0.28*exp(-0.062*(Vs/mV-60)))**(-1)*(Vd-Vsyn) : amp
Hxs=0<=((Vs/mV-10)*1) : 1
Hxw=0<=((Vs/mV-20)*1) : 1
dSi/dt=(Hxs-Si/150)/ms : 1
'''
neurons = NeuronGroup(1, eqs)
neurons.Vs=-4.6*mV
neurons.Vd=-4.7*mV
neurons.h=.999
neurons.n=.0001
neurons.s=.009
neurons.c=.007
neurons.q=.01
neurons.Ca=0
neurons.Si=.3
M = StateMonitor(neurons, True ,record=True,dt=.1*ms)
run(500*ms, report='text')
plt.cla()
plt.clf()
plot(M.t/ms, M[0].Vs/mV)
plot(M.t/ms, M[0].Vd/mV)
show()
