SharedsimpleNeuronalCoupling.ipynbOpen in CoCalc
Tutorial to build simple networks of synaptically coupled neurons. The membrane potential dynamics are obtained using the thermodynamic transmembrane transport by Herrera-Valdez (PeerJ Preprints, 2017)

# Neuronal dynamics in small networks: Simple synaptic coupling¶

### Marco Arieli Herrera-Valdez$^1$, Fernando Pérez Díaz$^2$, César Flóres López$^1$, José Alberto Perez Benitez$^2$¶

#### $^2$ ESIME-SEPI, Instituto Politécnico Nacional¶

Consider a general minimal model for transmembrane potential for a single cell where Auxiliary functions

The amplitude of the synaptic current is described by a sum of an Ornstein-Uhlenbeck process representing non-specific inputs and a sum of synaptic currents that depend on the activity of the other cells. The parameters $v_x$, $g_x$, $s_x$, $r_x$ represent, respectively, the half-activation potential, gating charge, de-activation bias, and basal rate for activation of the gate represented by the variable $\quad x \in \left\{ m,w \right\}$.

# Import necessary modules
import scipy as sc
import matplotlib.pylab as gr
import time


### Evolution of the system assuming that synapses activate as functions of time¶

"""
Change in membrane potential for one neuron. The function depends on the states of the system, represented by the vector U, time, and a dictionary containing parameters for the cell
"""
def NaTKaDNaKaAMPAGabaA(u,t,p):
"""
Notes:
The voltages have been normalized y = v/vT and y* = v*/vT. Then dy/dt = dv/dt / vT
The current amplitudes for the currents should be already normalized with respect to the membrane capacitance.
A -> a = A/(vT Cm)
"""
y,w = u
yyKaD = p['gain_Act_KaD']*(y-p['y_Half_Act_KaD'])
wInf = expSigmoid_(yyKaD)
wRate = p['rate_Act_KaD'] * expSum_(yyKaD,s=p['bias_Act_KaD'])
mInf = expSigmoid_(p['gain_Act_NaT']*(y-p['y_Half_Act_NaT']))
jKaD= transmembraneFlux_(y,p["eta_KaD"]*p['y_Ka'], eta=p["eta_KaD"], rect=p["rect_KaD"], a=p["A_Single_KaD"])
jNaT= transmembraneFlux_(y,p["eta_NaT"]*p['y_Na'], eta=p["eta_NaT"], rect=p["rect_NaT"], a=p["A_Single_NaT"])
jNaKa= transmembraneFlux_(y,p['y_NaKa'], eta=p["eta_NaKa"], rect=p["rect_NaKa"], a=p["A_Single_NaKa"])
jAMPA= transmembraneFlux_(y,p['y_Ampa'], eta=p["eta_Ampa"], rect=p["rect_Ampa"], a=p["A_Single_Ampa"])
jGabaA= transmembraneFlux_(y,p['y_GabaA'], eta=p["eta_GabaA"], rect=p["rect_GabaA"], a=p["A_Single_GabaA"])
dy_Base = -p["N_NaT"]*mInf*(1-w)*jNaT - p["N_KaD"]*w*jKaD - p["N_NaKa"]* jNaKa
dy_Syn = - p['N_Ampa']*p['p_Ampa'](t) * jAMPA - p['N_GabaA']*p['p_GabaA'](t) * jGabaA - p["A_LFP"](t)
dw = (w**p['exp_Act_KaD']) * wRate * (wInf-w)
return sc.array([dy_Base+dy_Syn,dw])


#### Funciones auxiliares¶

# Function to calculate the transmembrane fluxes
def transmembraneFlux_(y,y_r, eta, rect=0.5, a=1.0, degree="exp"):
aa= (y_r- eta*y)
A= a*eta
if degree=="exp":
i = A *((sc.exp((rect-1) * aa)-(sc.exp((rect* aa)))))
elif degree=="3":
i = A * aa * (1+ ((1-2* rect )/2)* aa + ((3* rect **2 -3* rect +1)/6)* (aa**2))
elif degree=="1":
i =  A * aa
return i

def expSum_(y, s=0.5):
ee = sc.exp(y)
return ee**(s-1)+ ee**(s)

def expSigmoid_(y):
ee = sc.exp(y)
return ee/(1+ee)

def wInverse_(w,vHalf,gCharge,vT):
v= vHalf - (vT/gCharge) * sc.log( (1-w)/w)
return v

def expSub_(y, s=0.5):
ee = sc.exp(y)
return ee**(s-1) - ee**(s)



### Solución numérica de ecuaciones¶

##### Change of variables to simplify the numerics¶

Let $y= v/(C_m v_T)$ and $y_* = v_*/(C_m v_T)$. Then As a consequence, the amplitudes in the original system can also be normalized: ${N_x a_x} \mapsto A_x = \frac{N_x a_x}{C_m v_T}$

def normalizeVolts(parDict):
"""Normalize voltages with respect to the Boltzmann potential kT/q"""
nDict={}
vTCm= parDict['Cm'] * parDict['v_T']
nDict['vTCm']=vTCm
for k,i in parDict.items():
if ((k.find('v_')==0)):
#print(k)
nDict["y_"+k[2:]]=parDict[k]/parDict['v_T']

parDict.update(nDict)
return parDict

# Normalization of potentials and amplitudes for the numerical solution of the system
def normalizeAmps(parDict):
"""
Normalizes amplitudes for transmembrane currents with respect to the membrane capacitance.
"""
nDict={}
vC= parDict['Cm'] * parDict['v_T']
nDict['vTCm']=vC

for k,i in parDict.items():
if k.find('a_')==0:
str0="A_" + k[2:]
#if type(i)==float:
nDict[str0]= i/vC
#print(k,i); #print(str0, nDict[str0])
parDict.update(nDict)
return parDict



#### Numerical method: Fourth order Runge-Kutta¶

# ------------------------------------------
# Runge-Kutta numerics
# ------------------------------------------
def RK4(f):
return lambda y, t, dt: (
lambda dy1: (
lambda dy2: (
lambda dy3: (
lambda dy4: (dy1 + 2*dy2 + 2*dy3 + dy4)/6
)( dt * f( y + dy3, t + dt ) )
)( dt * f( y + dy2/2, t + dt/2 ) )
)( dt * f( y + dy1/2, t + dt/2 ) )
)( dt * f( y, t ) )

def solveRK4(p, f, ic, timeStep, nSteps):
"""
Example:
u=solveRK4(p, f,ic,timeStep, nSteps)
p is a dictionary containing the parameters for the system
"""
U=sc.zeros((nSteps, sc.prod(sc.shape(ic))),"float64")
if len(p.keys()) * len(p.values()) >0:
ff = lambda yy,tt: f(yy,tt,p)
dU = RK4(ff)
U[0]=ic
for n in sc.arange(0,nSteps-1):
changeU =dU( U[n], n*timeStep, timeStep)
U[n+1] = U[n] + changeU

return U.transpose()


def calc2DNrnOrbit(p,graphs=1,fs=16, maxT=20, show_wInv=0, plotCurrents=0, plotTimeCurrents=0, figName=''):
"""
Examples:
q=calc2DNrnOrbit(p,graphs=0)
q=calc2DNrnOrbit(p,graphs=1)
"""
orbit=solveRK4(p=p, f=p["rhs"],ic=p["ic"], timeStep=p["timeStep"], nSteps=p["nSteps"])
y = orbit[0]; v = y*p["v_T"]
w = orbit[1]
dvdt=sc.zeros(len(v)); dvdt[1:]= (v[1:]-v[:-1])/p["timeStep"]
dwdt=sc.zeros(len(v)); dwdt[1:]= (w[1:]-w[:-1])/p["timeStep"]
#wInv= wInverse(w,vHalf=p["v_Half_Act_KaD"],vT=p["v_T"],gCharge=p["gain_Act_KaD"])
print("max dv/dt=%g"%(dvdt.max()))
#print("orbit: ", orbit)
yyKaD = p['gain_Act_KaD']*(y-p['y_Half_Act_KaD'])
wInf = expSigmoid_(yyKaD)
wRate = p['rate_Act_KaD'] * expSum_(yyKaD,s=p['bias_Act_KaD'])
mInf = expSigmoid_(p['gain_Act_NaT']*(y-p['y_Half_Act_NaT']))
jKaD= transmembraneFlux_(y,p['y_Ka'], eta=p["eta_KaD"], rect=p["rect_KaD"], a=p["a_Single_KaD"],degree="exp")
jNaT= transmembraneFlux_(y,-p['y_Na'], eta=p["eta_NaT"], rect=p["rect_NaT"], a=p["a_Single_NaT"],degree="exp")
jNaKa= transmembraneFlux_(y,p['y_NaKa'], eta=p["eta_NaKa"], rect=p["rect_NaKa"], a=p["a_Single_NaKa"],degree="exp")
iNaT = p["N_NaT"]*(1-w)*mInf *jNaT
iKaD = p["N_KaD"]*w* jKaD
iNaKa = p["N_NaKa"]*jNaKa
iTot= (iNaT+ iKaD + iNaKa)
upInds = sc.where(dvdt>0)[0]
dnInds = sc.where(dvdt<0)[0]
rcNaT = iNaT.sum()/iTot.sum()
rcKd = iKaD.sum()/iTot.sum()
effNaT = iNaT[upInds].sum()/iTot[upInds].sum()
effKd = iKaD[dnInds].sum()/iTot[dnInds].sum()
print("Na efficiency = %g"%((effNaT)))
print("K efficiency = %g"%((effKd)))
#jAMPA = p["Cm"]*p['AMPA'](timeSamp) * nExpSub(y- p['y_r_AMPA'], s=p['bias_AMPA'])
#jGabaA = p["Cm"]*p['GabaA'](timeSamp) * nExpSub(y- p['y_r_GabaA'], s=p['bias_GabaA'])
p['timeSamples'] =sc.arange(0,p['timeMax'],p['timeStep'])
q={"timeSamples":p["timeSamples"],
"v":v, "w":w, "dvdt":dvdt, "dwdt":dwdt, #"wInv":wInv,
"mInf":mInf, "wInf":wInf,
"iNaT":iNaT, "iKaD":iKaD, "iNaKa":iNaKa,
"iTot":iTot, #"vNull":vNull,
#"jAMPA":jAMPA, "jGabaA":jGabaA
}
if graphs>0:
maxC = sc.array([iNaT.max(),iKaD.max()]).max()
f1=gr.figure(figsize=(17,5));
q["fig1"]=f1;
if plotCurrents>0:
f2=gr.figure(figsize=(17,5));q["fig3"]=f3;
q["fig2"]=f2
r=2; c=2; ax=list()
for n in range(r*c):
ax.append(f2.add_subplot(r,c,n+1))
gr.ioff()
# Generar los planos para las graficas
axV=f1.add_subplot(1,1,1)

# Asignar a cada plano una grafica
#xx=wInv*p["rate_Act_KaD"]
axV.plot(p['timeSamples'], v, 'k', lw=1, alpha=1.0, label=r'$(t,v)$');

ax[0].plot(dvdt, v, 'k', lw=1, alpha=1.0, label=r'$(\partial_t v,v)$');
ax[0].plot([dvdt.min(),dvdt.max()], [0,0], 'k:');
ax[1].plot(w, v, 'k', lw=1, alpha=1.0, label=r'$(w,v)$');
if show_wInv>0:
#axV.plot(p['timeSamples'], wInv, 'k', lw=3, alpha=0.35, label=r'$(t,w^{-1}(v))$');
ax[1].plot(w, wInv, 'k', lw=3, alpha=0.35, label=r'$(w,w^{-1}(v))$');
ax[0].plot(dvdt, wInv, 'k', lw=3, alpha=0.3, label=r'$(\partial_t v,w^{-1}(v))$');
if plotTimeCurrents>0:
ax.plot(p['timeSamples'], iKaD, 'r', lw=2, alpha=1.0, label=r'$(t,i_{Kd})$');
ax.plot(p['timeSamples'], iNaT, 'g', lw=2, alpha=1.0, label=r'$(t,i_{NaT})$');
ax.plot(p['timeSamples'], iNaKa, 'b', lw=2, alpha=1.0, label=r'$(t,i_{NaK})$');
ax.plot(p['timeSamples'],iTot, 'k', lw=1, alpha=1.0,  label=r'$(t,i_{M})$');
ax.plot([0,p['timeMax']], [0,0], 'k:');
ax.legend(ncol=1,loc="upper right",fontsize=fs)
ax.set_xlim(p["timeSamples"].min(),maxT)
ax.set_ylabel('nA',fontsize=fs); axC.set_xlabel('ms',fontsize=fs)
f3.text(xStart,0.975,'B',fontsize=fs)

ax[2].plot(dvdt, iKaD, 'r', lw=2, alpha=1.0, label=r'$(\partial_t v,i_{Kd})$');
ax[2].plot(dvdt, iNaT, 'g', lw=2, alpha=1.0, label=r'$(\partial_t v,i_{NaT})$');
ax[2].plot(dvdt, iNaKa, 'b', lw=2, alpha=1.0, label=r'$(\partial_t v,i_{NaK})$');
ax[2].plot(dvdt, iTot, 'k', lw=1, alpha=1.0, label=r'$(\partial_t v,i_{M})$');
ax[2].plot([dvdt.min(),dvdt.max()], [0,0], 'k:');
ax[3].plot(w, iKaD, 'r', lw=2, alpha=1.0, label=r'$(\partial_t w,i_{Kd})$');
ax[3].plot(w, iNaT, 'g', lw=2, alpha=1.0, label=r'$(\partial_t w,i_{NaT})$');
ax[3].plot(w, iNaKa, 'b', lw=2, alpha=1.0, label=r'$(\partial_t w,i_{NaK})$');
ax[3].plot(w, iTot, 'k', lw=1, alpha=1.0, label=r'$(\partial_t w,i_{M})$');
ax[3].plot([w.min(),w.max()], [0,0], 'k:');
axV.legend(ncol=1,loc="upper right",fontsize=fs)

ax[0].legend(ncol=1,loc="lower right",fontsize=fs)
ax[1].legend(ncol=1,loc="lower right",fontsize=fs)
ax[2].legend(ncol=1,loc="upper right",fontsize=fs)
ax[3].legend(ncol=1,loc="upper left",fontsize=fs)
ax[2].set_ylim(-maxC,maxC)
ax[3].set_ylim(-maxC,maxC)
axV.set_ylabel('$v$ (mV)',fontsize=fs); axV.set_xlabel('ms',fontsize=fs);
axV.set_xlim(p["timeSamples"].min(),maxT);
ax[0].set_ylabel('$v$ (mV)',fontsize=fs); ax[0].set_xlabel('V/s',fontsize=fs);
ax[1].set_ylabel('$v$ (mV)',fontsize=fs); ax[1].set_xlabel('$w$',fontsize=fs);
ax[2].set_ylabel('pA',fontsize=fs); ax[2].set_xlabel('V/s',fontsize=fs)
ax[3].set_ylabel('pA',fontsize=fs); ax[3].set_xlabel('$\partial_t w$ (1/s)',fontsize=fs)
f1.subplots_adjust(left=0.05, bottom=0.05, right=0.95, top=0.95, wspace=0.2, hspace=0.13)
f2.subplots_adjust(left=0.05, bottom=0.05, right=0.95, top=0.95, wspace=0.2, hspace=0.13)
xStart=0.01; xDiv=0.485; yDiv=0.485
f1.text(xStart,0.975,'A',fontsize=fs)
f2.text(xStart,0.975,'C',fontsize=fs)
f2.text(xStart+xDiv,0.975,'D',fontsize=fs)
f2.text(xStart,yDiv,'E',fontsize=fs)
f2.text(xStart+xDiv,yDiv,'F',fontsize=fs)
gr.ion(); gr.draw();
if len(figName)>0:
n1=figName[:-4]+'_p1.png'
n2=figName[:-4]+'_p2.png'
f1.savefig(n1)
f2.savefig(n2)
print("Saved figure to files %s and %s"%(n1,n2))
if plotCurrents:
n3=figName[:-4]+'_currents.png'
f3.savefig(n3)
print("and the currents to %s and %s"%(n3))

return q


### Parameters¶

# Dictionaries with parameters and functions
pars={'Cm':20.0, 'v_T':26.7,
'gain_Act_NaT':5.0, 'v_Half_Act_NaT': -17.0,
'gain_Act_KaD':3.5, 'v_Half_Act_KaD': 0.0, 'rate_Act_KaD': 2.0, 'bias_Act_KaD': 0.3,'exp_Act_KaD':1.0,
'v_Ka':-90.0, 'v_Na':60.0, 'v_ATP':-430.0, 'v_Ampa':0.0, 'v_GabaA':-70.0,
'a_Single_NaT':1.2, 'a_Single_KaD':1.0,'a_Single_NaKa':1.0, 'a_Single_Ampa':1.0,'a_Single_GabaA':1.0,
'N_NaT':1400.0, 'N_KaD':4500.0, 'N_NaKa':200.0, 'N_Ampa':0.0, 'N_GabaA':0.0,
'eta_NaT':-1.0, 'eta_KaD':1.0, 'eta_NaKa':1.0, 'eta_Ampa':-1.0, 'eta_GabaA':1.0,
'rect_KaD':0.5,'rect_NaT':0.5,'rect_NaKa':0.5,'rect_Ampa':0.5,'rect_GabaA':0.5,
}
pars["v_NaKa"]= 3*pars["v_Na"] - 2*pars["v_Ka"] + pars["v_ATP"]
funcs ={"LFP": lambda t: 0.0}
descriptions ={'Assumptions':"Potassium KD channel activation describes NaT channel inactivation.",
'gain_Act_NaT': "Ganancia en el estado estable de activación de canales de sodio",
'gain_Act_KaD': "Ganancia en el estado estable de activación de canales de potasio",
'v_Half_Act_NaT': "Voltaje de activación media de canales de sodio",
'v_Half_Act_KaD': "Voltaje de activación media de canales de potasio",
'v_ATP': "Voltaje de reversa de la ATPasa de Na-K"}

##### Calculation of the solution¶

Initial conditions and number of steps

pars["nCells"]=1
v0=-30.0/pars["v_T"]
y0= 0.001
pars["ic"]= sc.array([v0,y0])
pars["timeMax"] = 40.0
pars["timeStep"]=1/70.0; pars["nSteps"]=10000;
pars["nSteps"]= sc.int32(pars["timeMax"]/pars["timeStep"])
pars["timeSamples"]=sc.arange(0,pars["timeMax"],pars["timeStep"])
pars['N_Ampa']=10.0;
pars["p_Ampa"]= lambda t: sc.int32((t>20)&(t<22)) + sc.int32((t>30)&(t<32))
pars["p_GabaA"]= lambda t: 0.0
pars["A_LFP"] = lambda t: 0.0
pars["rhs"]=NaTKaDNaKaAMPAGabaA
pars= normalizeAmps(pars)
pars= normalizeVolts(pars)


The solution to the system is composed of two vectors that correspond to the two variables of the system, $y$ and $w$. Then y can be transformed back into $v$ by setting $v = y\cdot v_T$.

orbit=solveRK4(p=pars, f=pars["rhs"],ic=pars["ic"], timeStep=pars["timeStep"], nSteps=pars["nSteps"])
y = orbit[0]; v = y*pars["v_T"]
w = orbit[1]
dvdt=sc.zeros(len(v)); dvdt[1:]= (v[1:]-v[:-1])/pars["timeStep"]
gr.figure(figsize=(15,5))
ax_tv=_subplot(121); ax_dvdt=f1.add_subplot(122)
ax_tv.plot(pars["timeSamples"],v,label=r'$(t,v)$')
ax_dvdt.plot(dvdt,v,label=r'$(\partial_t v,v)$')
ax_tv.legend(); ax_dvdt.legend()
#q=calc2DNrnOrbit(pars,graphs=1,fs=16, maxT=40, show_wInv=0, plotCurrents=0, figName='')

--------------------------------------------------------------------------- NameError Traceback (most recent call last) <ipython-input-9-8e22c940ab5e> in <module>() 4 dvdt=sc.zeros(len(v)); dvdt[1:]= (v[1:]-v[:-1])/pars["timeStep"] 5 gr.figure(figsize=(15,5)) ----> 6 ax_tv=_subplot(121); ax_dvdt=f1.add_subplot(122) 7 ax_tv.plot(pars["timeSamples"],v,label=r'$(t,v)$') 8 ax_dvdt.plot(dvdt,v,label=r'$(\partial_t v,v)$') NameError: name '_subplot' is not defined 
<matplotlib.figure.Figure at 0x7fded3a2c990>

### Network dynamics with simple synaptic coupling¶

The activation of synapses depends only on the voltage of the presynaptic neuron. Assume for simplicity that there are $N=2$ neurons in the network. Construct two matrices for the connections such that the first neuron excites the second, and the second neuron excites the first:

pars["gain_Act_AMPA"] = 5.0; pars["v_Half_Act_AMPA"] = -20.0
pars["gain_Act_GabaA"] = 5.0; pars["v_Half_Act_GabaA"] = -20.0
pars=normalizeVolts(pars)

# Connection matrices
pars["conn_Ampa"]= sc.matrix([[0,1],[0,0]])
print(pars["conn_Ampa"])
print(sc.array([1,2]))
sc.dot(pars["conn_Ampa"],sc.array([1,2]))

"""
Change in membrane potential for one neuron. The function depends on the states of the system, represented by the vector U, time, and a dictionary containing parameters for the cell
"""
def NaTKaDNaKaAMPAGabaA_Net(u,t,p):
"""
Notes:
The voltages have been normalized y = v/vT and y* = v*/vT. Then dy/dt = dv/dt / vT
The current amplitudes for the currents should be already normalized with respect to the membrane capacitance.
A -> a = A/(vT Cm)
"""
y= u[0:pars["nCells"]]
w= u[pars["nCells"]: 2*pars["nCells"]]
yyKaD = p['gain_Act_KaD']*(y-p['y_Half_Act_KaD'])
wInf = expSigmoid_(yyKaD)
wRate = p['rate_Act_KaD'] * expSum_(yyKaD,s=p['bias_Act_KaD'])
mInf = expSigmoid_(p['gain_Act_NaT']*(y-p['y_Half_Act_NaT']))
jKaD= transmembraneFlux_(y,p["eta_KaD"]*p['y_Ka'], eta=p["eta_KaD"], rect=p["rect_KaD"], a=p["A_Single_KaD"])
jNaT= transmembraneFlux_(y,p["eta_NaT"]*p['y_Na'], eta=p["eta_NaT"], rect=p["rect_NaT"], a=p["A_Single_NaT"])
jNaKa= transmembraneFlux_(y,p['y_NaKa'], eta=p["eta_NaKa"], rect=p["rect_NaKa"], a=p["A_Single_NaKa"])
jAMPA= transmembraneFlux_(y,p['y_Ampa'], eta=p["eta_Ampa"], rect=p["rect_Ampa"], a=p["A_Single_Ampa"])
jGabaA= transmembraneFlux_(y,p['y_GabaA'], eta=p["eta_GabaA"], rect=p["rect_GabaA"], a=p["A_Single_GabaA"])
dw = (w**p['exp_Act_KaD']) * wRate * (wInf-w)
pAmpa_ss = sc.dot(pars["conn_Ampa"], expSigmoid_(p['gain_Act_Ampa']*(y-p['y_Half_Act_Ampa'])))
pGabaA_ss = sc.dot(pars["conn_GabaA"], expSigmoid_(p['gain_Act_GabaA']*(y-p['y_Half_Act_GabaA'])))
pGabaA_ss = expSigmoid_(p['gain_Act_GabaA']*(y-p['y_Half_Act_GabaA']))
dy_Base = -p["N_NaT"]*mInf*(1-w)*jNaT - p["N_KaD"]*w*jKaD - p["N_NaKa"]* jNaKa
dy_Syn = - p['N_Ampa']*p['p_Ampa'](t) * jAMPA - p['N_GabaA']*p['p_GabaA'](t) * jGabaA - p["A_LFP"](t)
return sc.array([dy_Base+dy_Syn,dw])


pars["nCells"]=2
v0=-sc.ones(pars["nCells"])*30.0/pars["v_T"]
y0= sc.ones(pars["nCells"])
pars["ic"]= sc.array([v0,y0])
print(pars["ic"])

pars["nCells"]=2
v0=-sc.ones(pars["nCells"])*30.0/pars["v_T"]
y0= sc.ones(pars["nCells"])
pars["ic"]= sc.array([v0,y0])
print(pars["v0"])
pars["timeMax"] = 40.0
pars["timeStep"]=1/70.0; pars["nSteps"]=10000;
pars["nSteps"]= sc.int32(pars["timeMax"]/pars["timeStep"])
pars["rhs"]=NaTKaDNaKaAMPAGabaA
pars= normalizeAmps(pars)
pars= normalizeVolts(pars)