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-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 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 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\}.


# 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 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"]
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=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:

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)