Hi @foneill. That’s an interesting paper, thanks for sharing. I tried changing a few parameters here and there but couldn’t get any closer to those of the original paper. It looks like reproducing izhikevich neuron can be tricky sometimes, so I modified your code with the model mentioned in this post. Note that the authors in the paper you mentioned used different parameters (delays, dt, and simulation time), so I changed that too.
import matplotlib.pyplot as plt
import numpy as np
from brian2 import NeuronGroup, Synapses, SpikeMonitor, StateMonitor
from brian2 import ms, mV
from brian2 import defaultclock, run
import networkx as nx
def initSmallWorldConnections(k, num_excitatory, num_inhibitory):
excitatory_connectivity_matrix = np.zeros((num_excitatory, num_excitatory))
neuron_indices = np.array(range(0, num_excitatory))
# Split the (excitatory) neurons into modules
neighbourhoods = np.array_split(neuron_indices, k)
# Create a dictionary to map neuron index to module index
neuron_to_module = {}
# Iterate through modules and establish 10% connectivity (within modules)
for module_index, n in enumerate(neighbourhoods):
for i in n:
j = np.random.choice(n, 10, replace=False)
excitatory_connectivity_matrix[i, j] = 1
# Populate the dictionary with neuron index and module index
neuron_to_module.update({i: module_index})
# Introduce the pool of inhibitory neurons
connectivity_matrix = np.pad(excitatory_connectivity_matrix, ((0, num_inhibitory), (0, num_inhibitory)), mode='constant')
# Create focal connections for every four excitatory neurons
current_inhib_idx = num_excitatory
for n in neighbourhoods:
sub_neighbourhoods = np.array_split(n, 25)
for sn in sub_neighbourhoods:
if current_inhib_idx > (num_inhibitory + num_excitatory) - 1:
current_inhib_idx = num_excitatory
for i in sn:
connectivity_matrix[i, current_inhib_idx] = 1
current_inhib_idx += 1
# Create the diffusion connections from inhibitory to excitatory
for i in range(num_excitatory, num_excitatory + num_inhibitory):
for j in range(num_excitatory + num_inhibitory):
connectivity_matrix[i, j] = 1
# Remove reflexive connections
np.fill_diagonal(connectivity_matrix, 0)
return connectivity_matrix, neuron_to_module
# Function to rewire connections according to a probability p.
# If rewired, excitatory neuron disconnects from its target and connects to other excitatory neuron from distant module.
def rewireConnections(p, connectivity_matrix, k, num_excitatory):
tmp = np.copy(connectivity_matrix)
neuron_indices = np.array(range(0,num_excitatory))
# split the (excitatory) neurons into modules
neighbourhoods = np.array_split(neuron_indices, k)
for n_idx, n in enumerate(neighbourhoods):
#filter to only other modules
distant_neurons = neighbourhoods.copy()
distant_neurons.pop(n_idx)
distant_neurons = np.array(distant_neurons)
dis_neuron_flatten = distant_neurons.flatten()
for neuron in n:
excitatory_targets = np.where(tmp[neuron, :num_excitatory] == 1)[0]
rewired_targets = [x for x in excitatory_targets if np.random.rand()<p]
tmp[neuron, rewired_targets] = 0
modified_connectivity_matrix, new_connection = initRandomDistantConnections(neuron, dis_neuron_flatten, tmp, len(rewired_targets))
return modified_connectivity_matrix
# used above
def initRandomDistantConnections(neuron_idx, distant_neurons, connectivity_matrix, new_conns):
new_connection = np.random.choice(distant_neurons, new_conns, replace=False)
#print(f"Neuron {neuron_idx} has made a connection with Neuron {new_connection}\n")
connectivity_matrix[neuron_idx, new_connection] = 1
return connectivity_matrix, new_connection
tfinal = 50000 * ms
Ne = 800
Ni = 200
connections, neuron_to_module = initSmallWorldConnections(8, Ne, Ni)
print("Building connections...")
modified_connections = rewireConnections(0.9, connections, 8, Ne)
print("Print Applying rewiring...")
weights = np.hstack(
[
0.5*np.random.uniform(size=(Ne + Ni, Ne)),
np.random.uniform(size=(Ne + Ni, Ni)),
]
).T
modified_connections *= weights
defaultclock.dt = 0.2 * ms
eqs = """
dv/dt = (0.04*v**2+5.*v+140- u + g_exc - g_inh + noise * randn())/ms : 1 (unless refractory)
du/dt = a*(b*v-u)/ms : 1 (unless refractory)
g_exc :1
g_inh :1
noise :1
a :1
b :1
c :1
d :1
"""
N = NeuronGroup(Ne + Ni, eqs, threshold="v>=30", refractory=0.1*ms,
reset="v=c; u+=d; g_exc=0; g_inh=0", method="euler")
N.v = -65
N_exc = N[:Ne]
N_inh = N[Ne:]
spikemon = SpikeMonitor(N)
statemon = StateMonitor(N, ['v', 'u', 'g_exc', 'g_inh'], record=0, when='after_thresholds')
re = np.random.random(Ne)
ri = np.random.random(Ni)
N_exc.a = 0.02
N_exc.b = 0.2
N_exc.c = -65 + 15 * re**2
N_exc.d = 8 - 6 * re**2
N_exc.noise = 13.0
N_inh.a = 0.02 + 0.08 * ri
N_inh.b = 0.25 - 0.05 * ri
N_inh.c = -65
N_inh.d = 2
#N_inh.noise = 10.0
N_exc.u = "b*v"
N_inh.u = "b*v"
Se = Synapses(N_exc, N, model= "w:1", on_pre= 'g_exc_post += w')
Si = Synapses(N_inh, N, model= "w:1", on_pre= 'g_inh_post += w')
print("Connecting synapses...")
Se.connect(i=np.where(modified_connections[:Ne, :])[0],
j=np.where(modified_connections[:Ne, :])[1])
Si.connect(i=np.where(modified_connections[Ne:, :])[0],
j=np.where(modified_connections[Ne:, :])[1])
Se.delay = 'clip(20*rand(), 1, 20)*ms'
Si.delay = 1*ms
flatten_connections = modified_connections[:Ne, :].flatten()[modified_connections[:Ne, :].flatten() != 0]
Se.w[:] = flatten_connections
flatten_connections = modified_connections[Ne:, :].flatten()[modified_connections[Ne:, :].flatten() != 0]
Si.w[:] = flatten_connections
run(tfinal)
# Extract spike times and neuron indices for excitatory neurons only
spike_times = spikemon.t / ms # Convert to milliseconds
spike_indices = spikemon.i # Only include excitatory neurons
# Create the raster plot
plt.figure(figsize=(40, 24))
plt.scatter(spike_times, spike_indices, s=2, color='k')
plt.xlabel('Time (ms)')
plt.ylabel('Neuron index')
plt.title('Raster plot of excitatory neuron spikes')
plt.ylim(0, Ne+Ni)
plt.show()
plt.imshow(modified_connections)
plt.show()
In the code above I modified your functions a little bit because I had the feeling there were not enough connections in each module. The results are still not what I expected, but at least I could get something like a WTA activity. I hope that helps somehow