Hello, everybody!
I’m trying to implement at Brian2 neuron model proposed by Komendatov and Kononenko (equations with parameters at this article Effects of Synaptic Time Constant on Firing Activities of a minimal Central Pattern Generator | IEEE Conference Publication | IEEE Xplore ), but i have some trouble with physical units.
This is my current version of code:
from brian2 import *
from brian2.units.constants import faraday_constant as F
import numpy as np
start_scope()
defaultclock.dt = 0.01*ms
# Parameters of neuron
area = 20000*umetre**2
Cm = 0.02*ufarad*cm**-2 * area
gNAV = 0.12*siemens*cm**-2 * area
gNA = 0.032*siemens*cm**-2 * area
gK = 0.25*siemens*cm**-2 * area
gB = 0.115*siemens*cm**-2 * area
vNa = 40*mV
vK = -70*mV
vB = -58*mV
gNaTTX = 400*siemens*cm**-2 * area
gKTEA = 10*siemens*cm**-2 * area
gCa = 1*siemens*cm**-2 * area
vCa =150 *mV
gcaca = 0.01*siemens*cm**-2 * area
Kbetta = 15000 * 1/mM
betta = 0.00004 * mM
ro = 0.02
Ks = 50
# The model
neuron_eqs = Equations('''
dv/dt = (-InaV - Ina - Ik - Ib - InaTTX - IkTEA -Ica - Icaca + I)/Cm : volt
InaV = gNAV * (1/(1+exp(-0.2*1/mV*(v +45*mV))))*(v - vNa) : amp
Ina = gNA * (v - vNa) : amp
Ik = gK*(v - vK) : amp
Ib = gB*mB*hB*(v - vB) :amp
dmB/dt = (1/ (1 + exp(0.4*(v + 34*mV)/mV)) - mB)/0.05*hertz :1
dhB/dt = (1/(1 + exp(-0.55*(v + 43*mV)/mV)) - hB)/1.5*hertz :1
InaTTX = gNaTTX*m*m*m*h*(v- vNa) : amp
IkTEA = gKTEA*n*n*n*n*(v - vK) : amp
dm/dt = (1/(1+exp(-0.4*(v+31*mV)/mV))-m)/0.0005*hertz :1
dh/dt = (1/(1+exp(0.25*(v+45*mV)/mV))-h)/0.01*hertz :1
dn/dt = (1/(1+exp(-0.18*(v+25*mV)/mV))-n)/0.015*hertz :1
Ica = gCa*mCa*mCa*(v-vCa) : amp
dmCa/dt = (1/(1+exp(-0.2*v/mV))-mCa)/0.01*hertz :1
Icaca = (gcaca*(v-vCa))/((1+exp(-0.06*1/mV*(v+45*mV)))*(1+exp(Kbetta*(Ca*mM - betta)))) :amp
dCa/dt = ro*(-Ica/(2*F*V*mol/ms)-Ks*Ca)/ms :1
V = (4*pi*0.1*1/(mm*mm*mm)*mm*mm*mm)/3 :1
I : amp
''')
N = NeuronGroup(1, model=neuron_eqs, threshold='v>=-40*mV', refractory=3*ms, method='euler')
# Initialization
N.v = 0*mV
N.I = 1*amp
M = StateMonitor(N, 'v', record=True)
duration = 50*ms
run(duration, report = 'text')
plot(M.t/ms, M.v[0])
Iam be very pleasent if somebody helping with it.
I’am trying to reduce units from model as in original paper Deterministic Chaos in Mathematical Model of Pacemaker Activity in Bursting Neurons of Snail,Helix Pomatia - ScienceDirect
But calculation is very long (150 seconds with rk4 method and dt = 0.0001) . How to accelerate calculation?
What’s wrong at this code ? Could you @adam kindly help with this trouble?
from brian2 import *
import numpy as np
start_scope()
defaultclock.dt = 0.1*ms
# Parameters of neuron
Cm=0.02
R=0.1
F=96485
gK=0.25
gNaTTX=400.0
gKTEA=10.0
VNa=40.0
VK=-70.0
VB=-58.0
VCa=150.0
ks=50.0
rho=0.002
kbeta=15000
beta=0.00004
gNa=0.0231
gNaV=0.11
gB=0.1650
gCa=1.5
gCaCa=0.02
# The model
neuron_eqs = Equations('''
dV/dt=-(INaTTX+IKTEA+IK+INa+INaV+IB+ICa+ICaCa+Iapp)*hertz/Cm :1
dCa/dt=rho*(-ICa/(2*F*vol)-ks*Ca)/ms :1
dmB/dt=(1/(1+exp(0.4*(V+34)))-mB)/ms/0.05 :1
dhB/dt=(1/(1+exp(-0.55*(V+43)))-hB)/ms/1.5 :1
dm/dt=(1/(1+exp(-0.4*(V+31)))-m)/ms/0.0005 :1
dh/dt=(1/(1+exp(0.25*(V+45)))-h)/ms/0.01 :1
dn/dt=(1/(1+exp(-0.18*(V+25)))-n)/ms/0.015 :1
dmCa/dt=(1/(1+exp(-0.2*V))-mCa)/ms/0.01 :1
INaV=gNaV*(1/(1+exp(-0.2*(V+45))))*(V-VNa) :1
IK=gK*(V-VK) :1
INa=gNa*(V-VNa) :1
IB=gB*mB*hB*(V-VB) :1
INaTTX=gNaTTX*m*m*m*h*(V-VNa) :1
IKTEA=gKTEA *n*n*n*n*(V-VK) :1
ICa=gCa*mCa*mCa*(V-VCa) :1
ICaCa=gCaCa*(1/(1+exp(-0.06*(V+45))))*(1/(1+exp(kbeta*(Ca-beta))))*(V-VCa) :1
vol=4/3*pi*R*R*R :1
Iapp :1
''')
N = NeuronGroup(1, model=neuron_eqs, threshold='V>=-20', refractory=3*ms, method='gsl_rk4')
# Initialization
N.V = -42
N.Ca = 6e-05
N.mB = 0.95
N.hB = 0.77
N.m = 0.14
N.n = 0.048
N.mCa = 0.0002
@network_operation(dt=defaultclock.dt)
def update_input(t):
if t>=25*second and t<25*second:
N.Iapp = -0.5
else:
N.Iapp = 0
M = StateMonitor(N, 'V', record=True)
duration = 50*second
run(duration, report = 'text')
plot(M.t/second, M.V[0])
adam
13 April 2022 15:25
4
Hi Serge,
At the moment I don’t have time to take a detailed look at your code. But one quick thing that I spotted which may help speed up the simulation is replacing your update_input network_operation with a TimedArray :
https://brian2.readthedocs.io/en/stable/resources/tutorials/3-intro-to-brian-simulations.html
Something like
Iapp_ta = TimedArray([0,-0.5,0], dt=25*second)
(then use Iapp_ta in the equations for neurons somewhere)
if t>=25*second and t<25*second:
which I suspect will always be False. Did you want the current on for a single dt timestep at t=25sec ? or on for 25sec < t < 50sec?
For the trouble with the physical units, you’ll have to be more specific about what it is you’re having trouble with. It looks like in the most recent version of your code you’ve made all quantities “unitless,” what error’s did you have when you tried to add units?
This code is my attempt to adapt the original article code for XXPOUT from the link ModelDB: Show Model
#################################### Control parameters ############################################
# Example 1. Chaotic activity
#p gB=0.1372
#p gNa=0.0231, gNaV=0.11, gCa=1.5, gCaCa=0.02, Iapp=0.0
# Example 2. Spiking (beating) activity
#p gB=0.11
# Example 3. Period-two spiking activity
# gB=0.124
# Example 4. Period-four spiking activity
# gB=0.130
# Example 5. Bursting activity
# p gB=0.1650
###################################################################################################
# Example 6. Mode transition from chaotic activity into bursting one evoked by short depolarization)
p gNa=0.0231,gNaV=0.11,gB=0.1372,gCa=1.5,gCaCa=0.02,Iapp=-0.5
# Example 7. Chaotic bursting mode
#p gNa=0.02, gNaV=0.13, gB=0.18, gCa=1.0, gCaCa=0.01, Iapp=0.0
############################# Fixed parameter #######################################################
number Cm=0.02,R=0.1,F=96485
p gK=0.25,gNaTTX=400.0,gKTEA=10.0
p VNa=40.0,VK=-70.0,VB=-58.0,VCa=150.0
p ks=50.0,rho=0.002,kbeta=15000,beta=0.00004
########################### Equations
vol=4/3*pi*R*R*R
Iappx=if((t>=70.0)&(t<=72.0))then(Iapp)else(0.0)
########################### Currents ####################################
INaV=gNaV*(1/(1+exp(-0.2*(V+45))))*(V-VNa)
IK=gK*(V-Vk)
INa=gNa*(V-VNa)
IB=gB*mB*hB*(V-VB)
INaTTX=gNaTTX*m*m*m*h*(V-VNa)
IKTEA=gKTEA *n*n*n*n*(V-VK)
ICa=gCa*mCa*mCa*(V-VCa)
ICaCa=gCaCa*(1/(1+exp(-0.06*(V+45))))*(1/(1+exp(kbeta*(Ca-beta))))*(V-VCa)
########################## Differential equations ####################
V'=-(INaTTX+IKTEA+IK+INa+INaV+IB+ICa+ICaCa+Iappx)/Cm
Ca'=rho*(-ICa/(2*F*vol)-ks*Ca)
mB'=(1/(1+exp(0.4*(V+34)))-mB)/0.05
hB'=(1/(1+exp(-0.55*(V+43)))-hB)/1.5
m'=(1/(1+exp(-0.4*(V+31)))-m)/0.0005
h'=(1/(1+exp(0.25*(V+45)))-h)/0.01
n'=(1/(1+exp(-0.18*(V+25)))-n)/0.015
mCa'=(1/(1+exp(-0.2*V))-mCa)/0.01
####################### Initial conditions ###########################
# initial conditions: Examples 1-6.
V(0)=-42
Ca(0)=6e-05
mB(0)=0.95
hB(0)=0.77
m(0)=0.14
n(0)=0.048
mCa(0)=0.0002
# initial conditions: Example 7, chaotic bursting
#V(0)=-55.56913
#Ca(0)=3.593358e-05
#mB(0)=0.0
#hB(0)=0.0
#m(0)=0.0
#n(0)=0.0
#mCa(0)=0.0
@ MAXSTOR=10000000
@ TOTAL=150.0
@ DT=0.0001
@ XLO=0.0, XHI=150.0, YLO=-65.0, YHI=55.0
done
I’ve tried your solution, but so far haven’t been able to get the change in potential in the form of bursts, as expected in the original article.
Could @mstimberg , if possible, see the implementation of this model in the simulator?
I found a bug and now everything works, thanks for your attention @adam @mstimberg