| Hosted by CoCalc | Download
Kernel: Python 3 (system-wide)

Neuronal dynamics in small networks: Simple synaptic coupling

Marco Arieli Herrera-Valdez1^1, Fernando Pérez Díaz2^2, Roberto García-Medina1^1, José Alberto Perez Benitez2^2

1^1 Facultad de Ciencias, Universidad Nacional Autónoma de México

2^2 ESIME-SEPI, Instituto Politécnico Nacional

Consider a general minimal model for transmembrane potential for a single cell where Cmtv=(1w)Fm(v)NNaT(v)ψNaT(v)wNKdψKd(v)NNaKψNaK(v)NAmpaψAmpa(v)NGabaAψGabaA(v)ILFP(t)tw=w(Fw(v)w)Rw(v)\begin{darray}{rcl} C_m \partial_t v &=& - (1-w) F_{m}(v) N_{NaT} (v) \psi_{NaT}(v) - w N_{Kd} \psi_{Kd}(v)- N_{NaK} \psi_{NaK}(v) - N_{Ampa} \psi_{Ampa}(v) - N_{GabaA} \psi_{GabaA}(v) - I_{LFP}(t) \\ \partial_t w &=& w \left( F_{w}(v) - w \right) R_w(v) \end{darray} Auxiliary functions ψx(v)=axδx[exp((sx1)vRxδxvvT)exp(sxvRxδxvvT)]x{NaT,Kd,NaK,Ampa,GabaA}Fx(v)=exp(gxvvxvT)1+exp(gxvvxvT)x{m,w}Rw(v)=rw[exp((sw1)gwvvwvT)+exp(swgwvvwvT)]\begin{darray}{rcl} \psi_x (v) &=& a_x \delta_x \left[ \exp \left( (s_x -1) \frac{v_{R_x}-\delta_x v}{v_T} \right) -\exp \left( s_x \frac{v_{R_x} -\delta_x v}{v_T} \right) \right] \quad x \in \left\{ NaT, Kd, NaK, Ampa, GabaA \right\} \\ F_x (v) &=& \frac{\exp \left( g_x \frac{v-v_x}{v_T} \right) }{1+\exp \left( g_x \frac{v-v_x}{v_T} \right) } \quad x \in \left\{ m,w \right\} \\ R_w (v) &=& r_w \left[ \exp \left( (s_w -1) g_w \frac{v-v_w}{v_T} \right) +\exp \left( s_w g_w \frac{v-v_w}{v_T} \right) \right] \end{darray}

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 vxv_x, gxg_x, sxs_x, rxr_x represent, respectively, the half-activation potential, gating charge, de-activation bias, and basal rate for activation of the gate represented by the variable x{m,w}\quad x \in \left\{ m,w \right\}.

Simple network with plasticity

tv=JNaK(v)JKd(v,u)JNaT(v,u)JCaN(v,u,c)JS(v,n,c)tu=u(F(vi,vu,gu)u)R(vi,vu,gu,su)tc=ccτcJCaN(v,c)tn=n(nnτn)n r(c,cm,k)\begin{darray}{rcl} \partial_t v &=& -J_{NaK}(v)- J_{Kd}(v,u)- J_{NaT}(v,u)- J_{CaN}(v,u,c)- J_{S}(v,n,c) \\ \partial_t u &=& u \left(F(v_i,v_u,g_u)- u\right) R(v_i,v_u,g_u,s_u) \\ \partial_t c &=& \frac{c_{\infty}-c}{\tau_c} - J_{CaN}(v,c) \\ \partial_t n &=& n \left( \frac{n_{\infty}-n}{\tau_{n}}\right) - n~r(c,c_m,k) \end{darray}

with auxiliary functions ParseError: KaTeX parse error: Undefined control sequence: \label at position 222: …,u,E \right\}, \̲l̲a̲b̲e̲l̲{uSS} \\ R(v,v_…

\paragraph{Synapses between PINs} Let the connection matrix be w=(w11w1RwR1wRR)\begin{equation} w = \left( \begin{array}{ccc} w_{11}& \cdots & w_{1R} \\ \vdots& \ddots & \vdots \\ w_{R1}& \cdots & w_{RR} \end{array} \right) \end{equation} where ParseError: KaTeX parse error: Undefined control sequence: \lrSet at position 12: w_{ij} \in \̲l̲r̲S̲e̲t̲{0,1} represents the connection weight from neuron jj to neuron ii. For simplification purposes, the activation of the presynaptic terminals in the iith neuron is assumed to depend on the transmembrane potential of the membrane potential of the iith neuron.

The excitatory postsynaptic current for the iith neuron is given by JEi(v)=aEsinh(vivE2vT)j=1Rwijr(cj,cmj,k)nj\begin{darray}{rcl} J_{E i}(v) = a_{E} \sinh \left( \frac{v_i - v_{E}}{2v_T} \right) \sum_{j=1}^R w_{ij} r(c_j, c_{mj}, k) n_j \end{darray} where aE(v)=w b=(j=1Rw1jr(cj,cmj,k)njj=1RwRjr(cj,cmj,k)nj)\begin{darray}{rcl} a_{E}(v) &=& w~b \\ &=& \left( \begin{array}{c} \sum_{j=1}^R w_{1j}r(c_j, c_{mj}, k) n_j\\ \vdots \\ \sum_{j=1}^R w_{Rj} r(c_j, c_{mj}, k) n_j\\ \end{array} \right) \end{darray}

# 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/(CmvT)y= v/(C_m v_T) and y=v/(CmvT)y_* = v_*/(C_m v_T). Then ty=tvCmvT\begin{darray}{rcl} \partial_t y = \frac{\partial_t v}{C_m v_T} \end{darray} As a consequence, the amplitudes in the original system can also be normalized: NxaxAx=NxaxCmvT {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, yy and ww. Then y can be transformed back into vv by setting v=yvTv = 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"] f1=gr.figure(figsize=(15,5)) ax_tv=f1.add_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='')

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=2N=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: CAmpa=[0100],CGabaA=[0010]\begin{equation} 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] \end{equation}

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]) 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)
--------------------------------------------------------------------------- NameError Traceback (most recent call last) <ipython-input-2-731127457a27> in <module>() ----> 1 pars["nCells"]=2 2 v0=-sc.ones(pars["nCells"])*30.0/pars["v_T"] 3 y0= sc.ones(pars["nCells"]) 4 pars["ic"]= sc.array([v0,y0]) 5 pars["timeMax"] = 40.0 NameError: name 'pars' is not defined