<jupyter_start><jupyter_text>**The Izhikevich Model**(1) dv/dt = 0.04v^2 + 5v + 140 - u + I\(2) du/dt = a(bv-u)\(3) if v = 30mV, then: v <- c; u <- u+d**where:**\v = membrane potential \u = recovery variable \a = time scale of recovery \b = sensitivity of recovery to subthreshold oscillations \c = after-spike reset value of v \d = after-spike reset increment of u<jupyter_code>def thalamic_input(signal):
    out_signal = np.zeros_like(signal)
    i = np.nonzero(signal)[0]
    spikes = np.random.choice(len(i), int(np.floor(0.6 * len(i))), replace=False)
    out_signal[i[spikes]] = 1
    return out_signal

def timeStep(recov_scale, recov_sensitivity, voltage_reset, recov_reset,
             voltage, recov, inp, thresh=30, substeps=2):
    """ worker function for simulation. given parameters and current state variables, compute next ms """
    
    fired = voltage > thresh # array of indices of spikes
    #print('fired', sum(fired))
    voltage_next = voltage # next step of membrane potential v
    recov_next = recov # next step of recovery variable u
    
    ### Action potentials ###
    voltage[fired] = thresh
    voltage_next[fired] = voltage_reset[fired] # reset the voltage of any neuron that fired to c
    recov_next[fired] = recov[fired] + recov_reset[fired] # reset the recovery variable of any neuron that fired to u_0 + d
    
    # sum spontanous thalamic input and weighted inputs from all other neurons
    inp_next = inp + np.sum(connect[:,fired], axis=1) 
    
    ### Step forward ###
    for i in range(substeps):  # for numerical stability, execute at least two substeps per ms
        voltage_next += (1.0/substeps) * (0.04*(voltage_next**2) + (5*voltage_next) + 140 - recov + inp_next)
        recov_next += (1.0/substeps) * recov_scale * (recov_sensitivity*voltage_next-recov_next)
    return voltage_next, recov_next, fired

def simulate(recov_scale, recov_sensitivity, voltage_reset, recov_reset,
             voltage, recov, length, verbose=True):
    global all_inp, img
    """ input:
    - for each of N neurons: parameters recov_scale, recov_sensitivity, voltage_reset, recov_reset,
                             voltage, recov
    - length of simulation in milliseconds
    - NxN synaptic weight matrix
    
    processing:
        - simulates network evolution with spike-timing-dependent plasticity
    
    output:
        - Nxlength matrix of membrane voltages over time
        - Nxlength matrix of recovery variables over time
        - Nxlength matrix of spikes over time
    """
    voltage_out = np.zeros((N,length),dtype=np.double)
    voltage_out[:,0] = voltage  # initial voltages
    recov_out = np.zeros((N,length),dtype=np.double)
    recov_out[:,0] = recov
    firings_out = np.zeros((N,length),dtype=np.double)
    
    t0 = time.clock()
    
    ## simulate
    for t in range(1,length):
        random_inp = thalamic_input(filtered) + np.random.rand(N) * 0.3 # thalamic input
        #random_inp = random_inp > 0.7
        all_inp += random_inp
        all_inp -= all_inp * decay
        #Itrace[t,:] = It
        
        voltage_out[:,t], recov_out[:,t], firings_out[:,t] = timeStep(recov_scale, 
                                                                    recov_sensitivity, 
                                                                    voltage_reset, 
                                                                    recov_reset,
                                                                    voltage_out[:,t-1],
                                                                    recov_out[:,t-1],
                                                                    all_inp)
        
        # report progress
        if verbose and t % 100 == 0:
            print("Simulated " + str(t) + "ms of braintime in " + str(time.clock()-t0) + "s of computer time.") 
    
    t1 = time.clock()
    print("Simulation took " + str((t1-t0)) + "s")
    return voltage_out, recov_out, firings_out
#spatial distances

l = np.array(list(itertools.product(range(height), range(width))))
distances = sp.spatial.distance.cdist(l,l,'chebyshev')
distances = np.tile(distances, (len(filters),len(filters)))
spatial_connect = 1/(distances+1)
#spatial_connect[spatial_connect != 0.5] = 0 #connect only immediate neighbors
inp_strength, connect_strength = 1, 1.5

#angle differences

def angle_diff(x,y):
    x,y = sorted([x,y])
    abs_diff = 180 - abs(abs(x-y) - 180)
    if (x == 0 and y > 90) or (x < 90 and y == 180):
        return (180 - abs_diff)/45
    return abs_diff/45   

vec_len = width*height
angles = [180//4*i for i in range(len(filters))]
angle_diffs = np.zeros((vec_len*len(filters), vec_len*len(filters)))

for i in range(len(angles)):
    for j in range(len(angles)):
        angle_diffs[vec_len*i:vec_len*(i+1),vec_len*j:vec_len*(j+1)] = angle_diff(angles[i], angles[j])
angle_connect = 1/(angle_diffs+1)        
        
#connections
connect = angle_connect * spatial_connect * connect_strength
connect.shape
#image generation
inp_strength = 1
img = np.zeros(shape=(height,width))

rr, cc, _ = line_aa(0, int(width*3/8), height-1, int(width*3/8))  
img[rr, cc] = inp_strength

rr, cc, _ = line_aa(0, int(width/2), height-1, int(width/2))  
img[rr, cc] = inp_strength

rr, cc, _ = line_aa(0, int(width*5/8), height-1, int(width*5/8))  
img[rr, cc] = inp_strength

rr, cc, _ = line_aa(int(height*3/8), 0, int(height*3/8), width-1)  
img[rr, cc] = inp_strength

rr, cc, _ = line_aa(int(height/2), 0, int(height/2), width-1)  
img[rr, cc] = inp_strength

rr, cc, _ = line_aa(int(height*5/8), 0, int(height*5/8), width-1)  
img[rr, cc] = inp_strength

plt.imshow(img)
plt.show()
#image generation
img = np.zeros(shape=(height,width))

for i in range(5):
    #img[line_0(height, width)] = 1
    #img[line_90(height, width)] = 1
    img[line_45(height, width)] = 1
    img[line_135(height, width)] = 1

plt.imshow(img)
plt.show()
#applying convolutional filters
filtered_arr = get_angles(img)
filtered_arr[filtered_arr != 2] = 0
filtered_arr /= 2
filtered = filtered_arr.reshape(-1)

length = 1000    # simulation time in ms
N = height * width * len(filters) # total number of neurons

#firings = np.zeros((length,N))
re = np.array(np.random.rand(N),dtype=np.double)  # uniformly distributed random doubles [0,1)
    
recov_scale = 0.02+0.001*re     # a = time scale of recovery
recov_sensitivity = 0.2+0.001*re     # b = sensitivity of recovery to subthreshold oscillations
voltage_reset = -65+5*(re**2)  # c = after-spike reset value of v
recov_reset = 8-6*(re**2)    # d = after-spike reset increment of u
voltage = -65.0 * np.ones(N,dtype=np.double) # Initial values of v
recov = recov_sensitivity*voltage    # Initial values of u
    
all_inp = np.zeros(shape=(N,))
#Itrace  = np.zeros(shape=(length, N))
decay = 0.1

voltage_out, recov_out, firings_out = simulate(recov_scale,
                                               recov_sensitivity,
                                               voltage_reset,
                                               recov_reset,
                                               voltage,
                                               recov,
                                               length,
                                               verbose=True)
voltage_out1 = voltage_out.copy()
voltage_out1[voltage_out1 >= 30] = 30
#OUTPUT FOR HORIZONTAL AND VERTICAL LINES

h, w = voltage_out.shape #firings[:,1]

fig, ax = plt.subplots(figsize=(8,8))
plt.imshow(voltage_out1.T) 
#plt.imshow(firings_out.T) 
plt.colorbar()

ax.set_aspect(w/h)
plt.show()
s1 = []

for i in range(len(firings_out)):
    s = sum(firings_out[i] > 0)
    if s > 0:
        s1.append(s)
        
print('max', max(s1), 'min', min(s1), 'mean', np.mean(s1), 'median', np.median(s1))
#divide Izhikevich in model.py into Izhikevich and Simulation
from model import Izhikevich
from network import Connectivity
from utils import *
izh = Izhikevich(img, filters)
izh.simulat<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12245153220762976, 'fraction_numerical': 0.023639774859287054, 'mean_word_length': 2.496283340620901, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 26, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '262770', 'license': 'MIT', 'n_tokens_mistral': 2798, 'n_tokens_neox': 2626, 'n_words': 652, 'path': 'prototyping/izhikevich_matlab.ipynb', 'repo_name': 'two-7182/synchrony', 'size': '7995'}