Wow thank you both for your input! I have exactly the same setup as @adam (well, the variable names are different, of course), but knowing that someone more experienced has used this already puts me a bit at ease.

@mstimberg I understand your attached code, thank you for sharing! There are a few pieces here and there that I have not seen before (i.e. the `(summed)`

keyword in the equation definitions), but overall it looks much nicer (and neater) than my own implementation.

As for the [quote=āmstimberg, post:5, topic:496ā] `network_operation`

[/quote] ātrickā, I found out about it just today. My version is very hacky, but it works so far and I can leave it here for future reference (or for laughs).

```
from brian2 import *
# Kuramoto oscillators
kuramoto_eqs_stim = '''
dTheta/dt = ((omega + (kN * PIF) - I_stim*X*sin(Theta)) * second**-1) : 1
PIF = .5 * (sin(ThetaPreInput - Theta)) : 1
Vm = sin(Theta)*mV : volt
ThetaPreInput : 1
omega : 1
kN : 1
I_stim : amp
X = pulse_train(t) : amp**-1
'''
# parameters
duration = 5*second
defaultclock.dt = 1*ms
# Inputs setup
dt_stim = 1*ms
I0 = 10*amp
tv = linspace(0, duration/second, int(duration/(dt_stim))+1)
xstim = 1. * logical_and(tv>3, tv<3.1)
pulse_train = TimedArray(xstim*amp**-1, dt=dt_stim)
# Oscillators
seed(42)
N = 50
f0 = 4 # center freq [Hz]
sigma = 0.5 # normal std
net_kur = Network()
G = NeuronGroup(N, kuramoto_eqs_stim, threshold='True', method='euler', name='Kuramoto_N_%d' %N)
G.Theta = '2*pi*rand()' # uniform U~[0,2Ļ]
G.omega = '2*pi*(f0+sigma*randn())' # normal N~(f0,Ļ)
G.kN = 10
G.I_stim = I0
# hacks start here
G.namespace['order_param'] = zeros(int(duration/defaultclock.dt), dtype='complex')
G.namespace['cnt'] = 0
# synapses
S = Synapses(G, G, on_pre = '''ThetaPreInput_post = Theta_pre''', method='euler')
S.connect(condition='i!=j')
# monitors
M = StateMonitor(G, ['Theta'], record=True)
@network_operation(dt=1*ms)
def update_active():
G.namespace['order_param'][G.namespace['cnt']] = 1/G.N * sum(exp(1j*G.Theta)) # order parameter r(t)
G.namespace['cnt'] += 1
# add the above to the network
net_kur.add(G)
net_kur.add(S)
net_kur.add(M)
net_kur.add(update_active)
# run the simulation
net_kur.run(duration, report='text', report_period=10*second, profile=True)
print("Simulation done")
# plot the results
fig, axs = subplots(2,1)
axs[0].plot(M.t/second, mean(sin(M.Theta), axis=0))
axs[0].plot(M.t/second, sin(imag(G.namespace['order_param'])), '--')
axs[1].plot(M.t/second, abs(G.namespace['order_param']), '--', label='N=%d'%N)
# labels
axs[0].set_ylabel("Ensemble Theta Rhythm")
axs[1].set_ylabel("Kuramoto Order Parameter")
axs[1].set_xlabel("Time [s]")
axs[0].set_ylim([-1,1])
axs[1].set_ylim([0,1])
axs[0].axvline(x=3, ymin=-1, ymax=1, c="red", linewidth=2, zorder=0, clip_on=False)
axs[1].axvline(x=3, ymin=-1, ymax=1, c="red", linewidth=2, zorder=0, clip_on=True)
# make things pretty
legend()
grid()
# show
show()
```

I went around the imaginary number restrictions by defining the vector of interest as a numpy array that holds complex numbers as its data type. The vector is kept in the groupās namespace, which if Iām not mistaken is leveraged to the `run()`

workspace at simulation time. Itās a very crude way to do things, but it produces an interesting result, shown below.

Thank you both again for your time and detailed explanations! Your advice was extremely useful, and I learned a few new tricks!