Coupled Van der Pol oscillators

Support
2

Hi Cadgas. There’s actually nothing wrong with your code, it works! The only reason you are seeing nothing in your plot is that you are starting your oscillators at x=0, y=0, where they will not do anything interesting and stay at (0, 0). Here’s how to adapt your code to initialize them to other values and plot the results:

Neu = NeuronGroup(2, model=eqs, method='euler')
Neu.x = [0, 0.5]
Neu.y = [1, 0]
# ^ One oscilator starts at x=0, y=1, the other one at x=0.5, y=0
M = StateMonitor(Neu, ['x', 'y'], record=True)

run(100*ms, report='text')

plot(M.x.T, M.y.T)

A few minor comments:

  • you do not need to define a threshold in your NeuronGroup. This is not a model of a spiking neuron so no need to define when a “spike” is emitted.
  • you can use a single StateMonitor for both variables: M = StateMonitor(Neu, ['x', 'y'], record=True)
  • When you plot, do not divide things by mV. This is useful for the membrane potential of neuron models expressed as a voltage, but your model is only using dimensionless variables.

Also a remark on tau, which is always a source of confusion when looking at mathematical formulations like the one in the linked article which do not include time constants like you’d do for the equations of a physical system. When you read these differential equations, you can read them as being relative to some abstract “time unit”. To translate them to Brian, you have to divide the RHS of all differential equations by this “time unit”. This is what you’ve done here by introducing tau and therefore things work :+1:. Just to make the link to the mathematic formulation and e.g. the Wikipedia article more clear: By setting tau = 10*ms, you are defining your “unit of time” to be 10ms. This means when you simulate for 100ms, you simulate for “10 units of time”. If you want to plot something like https://en.wikipedia.org/wiki/Van_der_Pol_oscillator#/media/File:Vanderpol_mu=5.svg, you have to divide the time steps provided by Brian by 10ms to get the same time units on the x axis. Finally, note that Brian by default uses a simulation time step of 0.1ms (a reasonable time step for many neuron models), so with your definitions that means a simulation time step of 0.01 time units (often written as “h=0.01” in mathematical formulations). In your example, this seems to work well, but you might have to adapt it if you encounter problems.