Hello, I am learning Brian2 and tried to integrate two coupled Van der Pol oscillators but I think I don’t understand how to write model equations. I tried to plot x and y and I got a blank plot. Could you please help me?
I used the simplest model that starts with
x’ = y;
y’ = u(1-x^2)y - x;
u = 0.5;
Minimal code to reproduce problem
from brian2 import *
import numpy as np
import matplotlib.pyplot as plt
start_scope()
tau = 10*ms
eqs = '''
dx/dt = y/tau : 1
dy/dt = (0.5*(1-(x**2))*y-x)/tau : 1
'''
Neu = NeuronGroup(2, model=eqs,
threshold = 'x>1',
method='euler')
Mx = StateMonitor(Neu, 'x', record=True)
My = StateMonitor(Neu, 'y', record=True)
run(100*ms, report='text')
plot(Mx[0].x / mV, My[0].y / mV)
What you have aready tried
I tried to remove tau parts but I couldn’t manage it.
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 . 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.
Thank you very much! It helps a lot and I got my plot.
I also tried to connect two oscillators.
start_scope()
tau = 10*ms
eqs = '''
dx/dt = y/tau : 1
dy/dt = (0.5*(1-(x**2))*y-x)/tau : 1
'''
Neu = NeuronGroup(2, model=eqs, method='euler', threshold='True')
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)
S = Synapses(Neu,Neu, 'w : 1', on_pre='x += w')
S.connect(i = 0, j = 1, p = 1)
S.w = '0.01'
run(1000*ms, report='text')
plot(M.x.T, M.y.T, linewidth=3)
When I used this ‘w : 1’, on_pre=‘x += w’ code does it mean that it increases x with w for each step or does it connect two oscillators and each calculation step the second oscillator effected by x_2 += w*x_1? If it is the former one could you please help me to build a connected system with a w connectivity weight?
Hi. Brian’s main focus is on spike-based interactions, so the syntax for continuous interactions is a bit less convenient. There are three basic techniques to establish continuous interactions between neurons:
Use threshold=True to trigger an interaction (normally a “spike”) at every time step. This is the approach you used here
Use “linked variables” to get a new name for an existing variable, but referring to different neurons.
The last (linked variable) approach is the most efficient solution but is also less flexible. Most importantly you can only use it to couple single neurons, whereas the other two solutions extend to an arbitrary number of couplings per neuron (summing up the coupling terms), so I will not go into details about it here.
I am not 100% sure I understand your example with the third oscillator. But if you have two oscillators in your NeuronGroup (as in your example code above), and you want to add w times the difference of the two x variables as a coupling term, you can do it like this when using the threshold='True' approach (here I connect them in both directions, not sure whether this is what you want):
. Just to make the link to the mathematic formulation and e.g. the Wikipedia article more clear: By setting 
