In [1]:

    

### From Burak and Fiete, Plos Comp. Biol. 2009
### Author: Harsha Gurnani


    

In [2]:

    

from numpy import *
from scipy import *


    

In [ ]:

    




    

In [3]:

    

### Setting up model

nsedge = 128
NN = nsedge**2        #no. of neurons

posx = [ int(nn/nsedge)-nsedge/2 for nn in range(NN) ]
posy = [ int(mod(nn,nsedge))-nsedge/2 for nn in range(NN) ]

pref_dirn = [0 for ii in range(NN)]

dirns = [[0,1], [1,0],[0,-1], [-1,0]]
for ii in range(NN):
    shift=int(mod(mod(ii,8),4))
    ind= mod(shift+ int(mod(ii,4)),4)
    pref_dirn[ii] = ind
    
tau=10/1000
dt =0.5/1000


    

In [ ]:

    

### Weight matrix
def get_baselineW(a, gamma, beta, dist2):
    return a*exp(-gamma*dist2) - exp(-beta*dist2)

lambda_net=13
a = 1.0
beta =3/(lambda_net**2)
gamma = 1.05*beta
l = 0.5

vec=zeros(2)
Wt = zeros([NN,NN])
for ii in range(NN):
    for jj in range(NN):
        x1,y1=posx[ii],posy[ii]
        x2,y2=posx[jj],posy[jj]
        pf=dirns[pref_dirn[jj]]
        vec[0] = x1 - x2 - l*pf[0]
        vec[1] = y1 - y2 -l*pf[1]
        dist2 = vec[0]**2 + vec[1]**2
        W = get_baselineW(a,gamma,beta,dist2)
        Wt[ii,jj] = W


    

In [ ]:

    

### Velocity-modulated inputs

B = zeros(NN)
def get_vel_input(NN,A,alpha,edir,vel):
    bv = zeros(NN)
    for ii in range(NN):
        pf = dirns[edir[ii]]
        comp = pf[0]*vel[0] + pf[1]*vel[1]
        bv[ii] = A[ii]*(1+alpha*comp)
    return bv


    

In [ ]:

    

A = [1 for ii in range(NN)]
alpha=0.10315


    

In [ ]:

    

def DiffEqnForRates( s, t, Wt, tau , vel) :
    '''
    eqns for ds/dt
    s: Firing rates of all neurons
    Wt : synaptic weight matrix
    B : Velocity-modulated input
    '''
    Ns= len(s)
    n = int(sqrt(Ns))
    dsdt = zeros(Ns)
    B = get_vel_input(Ns, A, alpha, pref_dirn, vel[t])
    for ii in range(Ns):
        syn_act = B[ii]
        for jj in range(Ns):
            syn_act += Wt[jj,ii]*s[jj]
        if syn_act < 0:
            syn_act = 0;
        dsdt[ii] = (syn_act - s[ii])/tau
        
    return dsdt


    

In [ ]:

    

### Setting up simulations
dt =0.5/1000
t=10
trange=range(int(t/dt))

vel=[(rand(2)-0.5)*0.8 for ii in trange]
init_t=int(0.25/dt)
for ii in range(init_t):
    vel[ii] = [1,0]*0.8
for ii in range(init_t, 2*init_t):
    vel[ii] = [cos(pi/5), sin(pi/5)]*0.8
for ii in range(2*init_t,3*init_t):
    vel[ii] = [cos(pi/2-pi/5), sin(pi/2-pi/5)]*0.8

print shape(vel)

s0 = rand(NN)*5
rates = scipy.integrate.odeint(DiffEqnForRates, s0, t, Wt, tau, vel)
print shape(rates)


    

In [ ]:

    

### Get position from velocity
### Plot rates versus position
### Plot population response maps


    

In [ ]: