Description of problem
Hello All.
I am trying to implement a model from this paper.
A brief description of the network:
The model has 1000 neurons, 800 excitatory and 200 inhibitory.
Excitatory neurons are broken up into 8 modules, with 100 neurons each. Each model has a 10% edge density.
Inhibitory neurons exist in an inhibitory pool consisting of 200 neurons.
Focal connections exist from excitatory to inhibitory; every four neurons from the same excitatory module connects to the same inhibitory neuron.
Diffusion connections exist from inhibitory to excitatory; each inhibitory neuron connects to every other neuron in the network.
Finally, before initialization each excitatory neuron is subject to a rewiring process according to a probability p.
If the excitatory neuron is rewired: it detaches from it’s target and connects to a neuron in a distant module.
Different values of p dictate what regime the network undertakes. Low p values imply a competitive network, high p values imply a cooperative, synchronized network alike to a PING oscillator.
On page 5, Fig 2 of the paper you may observe a roster plot for different values of p. High values of p result in highly synchronized activity.
However, upon implementing this I am getting strange results. No matter the strength of p, the same level of activity and nothing like the raster plots seen in Fig 2.
I took the sample code of the Izhikevich model 2003 from Brian2’s website and applied the connectivity pattern explained above.
I rewired at p=0.9 as it is the most identifible raster plot pattern seen in Fig 2.
As you can see, my result is just random.
Minimal code to reproduce problem
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:
possible_targets = [j for j in n if i != j]
if possible_targets:
j = np.random.choice(possible_targets)
if np.random.rand() < 0.1:
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):
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:
#check if there is an inter-module relationship
distant_connections = np.nonzero(connectivity_matrix[neuron, dis_neuron_flatten] == 1)[0]
if np.random.rand() < p:
excitatory_targets = np.where(connectivity_matrix[neuron, :num_excitatory] == 1)[0]
#print(f"{neuron } has target {excitatory_targets}")
if len(distant_connections) != 0:
excitatory_targets = np.where(connectivity_matrix[neuron, :num_excitatory] == 1)[0]
random_connection = np.random.choice(distant_connections)
connectivity_matrix[neuron, excitatory_targets[0]] = 0
modified_connectivity_matrix, new_connection = initRandomDistantConnections(neuron, dis_neuron_flatten, connectivity_matrix)
return modified_connectivity_matrix
# used above
def initRandomDistantConnections(neuron_idx, distant_neurons, connectivity_matrix):
new_connection = np.random.choice(distant_neurons)
#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 = 1000 * 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...")
re = np.random.uniform(size=Ne)
ri = np.random.uniform(size=Ni)
weights = np.hstack(
[
0.5 * np.random.uniform(size=(Ne + Ni, Ne)),
-np.random.uniform(size=(Ne + Ni, Ni)),
]
).T
defaultclock.dt = 1 * ms
eqs = """dv/dt = (0.04*v**2 + 5*v + 140 - u + I + I_noise )/ms : 1
du/dt = (a*(b*v - u))/ms : 1
I : 1
I_noise : 1
a : 1
b : 1
c : 1
d : 1
"""
N = NeuronGroup(Ne + Ni, eqs, threshold="v>=30", reset="v=c; u+=d", method="euler")
N.v = -65
N_exc = N[:Ne]
N_inh = N[Ne:]
spikemon = SpikeMonitor(N)
statemon = StateMonitor(N, 'v', record=0, when='after_thresholds')
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_inh.a = 0.02 + 0.08 * ri
N_inh.b = 0.25 - 0.05 * ri
N_inh.c = -65
N_inh.d = 2
N_exc.u = "b*v"
N_inh.u = "b*v"
S = Synapses(
N,
N,
"w : 1",
on_pre={"up": "I += w", "down": "I -= w"},
delay={"up": 0 * ms, "down": 1 * ms},
)
print("Connecting synapses...")
S.connect(i=np.where(modified_connections)[0],
j=np.where(modified_connections)[1])
flatten_connections = modified_connections.flatten()[modified_connections.flatten() != 0]
S.w[:] = flatten_connections
N_exc.run_regularly("I_noise = 5*randn()", dt=1 * ms)
N_inh.run_regularly("I_noise = 2*randn()", dt=1 * ms)
run(tfinal)
# Extract spike times and neuron indices for excitatory neurons only
exc_spike_times = spikemon.t[spikemon.i < Ne] / ms # Convert to milliseconds
exc_spike_indices = spikemon.i[spikemon.i < Ne] # Only include excitatory neurons
# Create the raster plot
plt.figure(figsize=(40, 24))
plt.scatter(exc_spike_times, exc_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)
plt.show()
What you have aready tried
I have checked the firing rates for anything unusual. I took the mean firing rates of each module and aggregated. They seem to nearly always be the same for each time step.
I have tried verifying the connectivity pattern is correct multiple times now. I can confirm it is correct.
Expected output (if relevant)
I can’t actually include a second image as a new user. But the raster plot should be similar the third raster plot of Fig 2 in the paper attached above.
Actual output (if relevant)
What it actually looks like:
