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 suponiendo que las sinápsis se activan con respecto a funciones temporales

"""
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
"""
"""
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
mInf = expSigmoid_(p['gain_Act_NaT']*(y-p['y_Half_Act_NaT']))
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_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 $\partial_t y = \frac{\partial_t v}{C_m v_T}$ 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"]
print("max dv/dt=%g"%(dvdt.max()))
#print("orbit: ", orbit)
mInf = expSigmoid_(p['gain_Act_NaT']*(y-p['y_Half_Act_NaT']))
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
iNaKa = p["N_NaKa"]*jNaKa
upInds = sc.where(dvdt>0)[0]
dnInds = sc.where(dvdt<0)[0]
rcNaT = iNaT.sum()/iTot.sum()
effNaT = iNaT[upInds].sum()/iTot[upInds].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,
"iTot":iTot, #"vNull":vNull,
#"jAMPA":jAMPA, "jGabaA":jGabaA
}
if graphs>0:
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):
gr.ioff()
# Generar los planos para las graficas

# Asignar a cada plano una grafica
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,
'v_Ka':-90.0, 'v_Na':60.0, 'v_ATP':-430.0, 'v_Ampa':0.0, 'v_GabaA':-70.0,
}
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= 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.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: $C_{Ampa} =\left[ \begin{array}{cc} 0&1\\0&0 \end{array} \right], \quad C_{GabaA} =\left[ \begin{array}{cc} 0&0\\1&0 \end{array} \right]$

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
"""
"""
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"]]
mInf = expSigmoid_(p['gain_Act_NaT']*(y-p['y_Half_Act_NaT']))
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_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"])