The Van der Pol Oscillator#
Modeling relaxation oscillations and nonlinear damping in neural systems.
Learning Objectives#
By the end of this tutorial, you will be able to:
Understand the Van der Pol equations and the role of parameter \(\mu\)
Distinguish between harmonic oscillations (\(\mu \approx 0\)) and relaxation oscillations (\(\mu >> 1\))
Simulate Van der Pol oscillators across different parameter regimes
Explain why relaxation oscillations are relevant for neural dynamics
Background / Theory#
The Van der Pol Oscillator#
The Van der Pol oscillator was originally developed by Balthasar van der Pol in 1926 to model electrical circuits with vacuum tubes. It has since become a paradigmatic model for relaxation oscillations - oscillations characterized by slow buildup and rapid discharge phases.
Original Form (Second-Order)#
\[
\frac{d^2x}{dt^2} - \mu(1 - x^2)\frac{dx}{dt} + x = 0
\]
Two-Dimensional Form (Lienard Transformation)#
BrainMass uses the Lienard form:
\[\begin{split}
\begin{aligned}
\dot{x} &= \mu\left(x - \frac{x^3}{3} - y\right) + I_x(t) \\
\dot{y} &= \frac{x}{\mu} + I_y(t)
\end{aligned}
\end{split}\]
The Key Parameter: \(\mu\)#
\(\mu\) Value |
Oscillation Type |
Characteristics |
|---|---|---|
\(\mu \to 0\) |
Nearly harmonic |
Sinusoidal waveform |
\(\mu \approx 1\) |
Transition |
Mixed character |
\(\mu >> 1\) |
Relaxation |
Sharp transitions, flat plateaus |
Neuroscience Relevance#
The Van der Pol oscillator captures features of neural dynamics:
Slow-fast dynamics: Rapid action potentials followed by slower recovery
Threshold behavior: Excitation must exceed threshold for transition
Refractory period: Post-spike period of reduced excitability
It is closely related to the FitzHugh-Nagumo model (Tutorial 014).
import brainmass
import brainstate
import brainunit as u
import matplotlib.pyplot as plt
import numpy as np
# Set simulation time step
brainstate.environ.set(dt=0.1 * u.ms)
# Create Van der Pol oscillator
node = brainmass.VanDerPolStep(
1, # single node
mu=1.0, # moderate nonlinearity
)
# Initialize states
node.init_all_states()
# Set initial condition away from origin
node.x.value = np.array([0.5])
node.y.value = np.array([0.0])
print(f"mu = {node.mu.value()}")
mu = 1.0
# Define simulation step
def step_run(i):
node.update()
return node.x.value, node.y.value
# Run simulation (2000 ms for slower dynamics)
indices = np.arange(20000)
x_trace, y_trace = brainstate.transform.for_loop(step_run, indices)
# Define simulation step
def step_run(i):
node.update()
return node.x.value, node.y.value
# Run simulation (2000 ms for slower dynamics)
indices = np.arange(20000)
x_trace, y_trace = brainstate.transform.for_loop(step_run, indices)
t_ms = indices * brainstate.environ.get_dt()
fig, axes = plt.subplots(1, 2, figsize=(12, 4))
# Time series
axes[0].plot(t_ms, x_trace[:, 0], 'b-', label='x(t)', linewidth=1)
axes[0].plot(t_ms, y_trace[:, 0], 'r-', label='y(t)', linewidth=1, alpha=0.7)
axes[0].set_xlabel('Time (ms)')
axes[0].set_ylabel('Activity')
axes[0].set_title(f'Van der Pol Oscillator (mu={node.mu.value()})')
axes[0].legend()
axes[0].set_xlim([0, 1000])
# Phase portrait
axes[1].plot(x_trace[:, 0], y_trace[:, 0], 'k-', linewidth=0.5, alpha=0.7)
# Mark nullclines
x_null = np.linspace(-3, 3, 100)
y_null = x_null - x_null ** 3 / 3 # x-nullcline: y = x - x^3/3
axes[1].plot(x_null, y_null, 'g--', label='x-nullcline', linewidth=2)
axes[1].axvline(0, color='orange', linestyle='--', label='y-nullcline', linewidth=2)
axes[1].set_xlabel('x')
axes[1].set_ylabel('y')
axes[1].set_title('Phase Portrait with Nullclines')
axes[1].legend()
axes[1].set_xlim([-3, 3])
axes[1].set_ylim([-3, 3])
plt.tight_layout()
plt.show()
t_ms = indices * brainstate.environ.get_dt()
fig, axes = plt.subplots(1, 2, figsize=(12, 4))
# Time series
axes[0].plot(t_ms, x_trace[:, 0], 'b-', label='x(t)', linewidth=1)
axes[0].plot(t_ms, y_trace[:, 0], 'r-', label='y(t)', linewidth=1, alpha=0.7)
axes[0].set_xlabel('Time (ms)')
axes[0].set_ylabel('Activity')
axes[0].set_title(f'Van der Pol Oscillator (mu={node.mu.value()})')
axes[0].legend()
axes[0].set_xlim([0, 1000])
# Phase portrait
axes[1].plot(x_trace[:, 0], y_trace[:, 0], 'k-', linewidth=0.5, alpha=0.7)
# Mark nullclines
x_null = np.linspace(-3, 3, 100)
y_null = x_null - x_null ** 3 / 3 # x-nullcline: y = x - x^3/3
axes[1].plot(x_null, y_null, 'g--', label='x-nullcline', linewidth=2)
axes[1].axvline(0, color='orange', linestyle='--', label='y-nullcline', linewidth=2)
axes[1].set_xlabel('x')
axes[1].set_ylabel('y')
axes[1].set_title('Phase Portrait with Nullclines')
axes[1].legend()
axes[1].set_xlim([-3, 3])
axes[1].set_ylim([-3, 3])
plt.tight_layout()
plt.show()
@brainstate.transform.jit(static_argnums=1, static_argnames='n_step')
def simulate(mu_val, n_step: int = 30000):
model = brainmass.VanDerPolStep(1, mu=mu_val)
model.init_all_states()
model.x.value = np.array([0.5])
model.y.value = np.array([0.0])
# Simulate
def step(i):
model.update()
return model.x.value, model.y.value
return brainstate.transform.for_loop(step, np.arange(n_step))
# Different mu values
mu_values = [0.1, 0.5, 1.0, 3.0, 5.0]
fig, axes = plt.subplots(2, len(mu_values), figsize=(15, 6))
for idx, mu_val in enumerate(mu_values):
# Create oscillator
x, y = simulate(mu_val, 30000)
# Time series (top row)
t = np.arange(30000) * brainstate.environ.get_dt()
axes[0, idx].plot(t[-10000:], x[-10000:, 0], 'b-', linewidth=0.8)
axes[0, idx].set_title(f'$\mu$ = {mu_val}')
axes[0, idx].set_xlabel('Time (ms)')
if idx == 0:
axes[0, idx].set_ylabel('x(t)')
# Phase portrait (bottom row)
axes[1, idx].plot(x[-10000:, 0], y[-10000:, 0], 'k-', linewidth=0.3, alpha=0.7)
# Nullcline
x_null = np.linspace(-3, 3, 100)
axes[1, idx].plot(x_null, x_null - x_null ** 3 / 3, 'g--', linewidth=1.5, alpha=0.7)
axes[1, idx].set_xlabel('x')
if idx == 0:
axes[1, idx].set_ylabel('y')
axes[1, idx].set_xlim([-3, 3])
axes[1, idx].set_ylim([-3, 3])
plt.suptitle('Transition from Harmonic to Relaxation Oscillations', y=1.02, fontsize=14)
plt.tight_layout()
plt.show()
<>:13: SyntaxWarning: invalid escape sequence '\m'
<>:13: SyntaxWarning: invalid escape sequence '\m'
C:\Users\adadu\AppData\Local\Temp\ipykernel_32900\1202397000.py:13: SyntaxWarning: invalid escape sequence '\m'
axes[0, idx].set_title(f'$\mu$ = {mu_val}')
# Different mu values
mu_values = [0.1, 0.5, 1.0, 3.0, 5.0]
fig, axes = plt.subplots(2, len(mu_values), figsize=(15, 6))
for idx, mu_val in enumerate(mu_values):
# Create oscillator
model = brainmass.VanDerPolStep(1, mu=mu_val)
brainstate.nn.init_all_states(model)
model.x.value = np.array([0.5])
model.y.value = np.array([0.0])
# Simulate
def step(i):
with brainstate.environ.context(i=i, t=i * brainstate.environ.get_dt()):
model.update()
return model.x.value, model.y.value
x, y = brainstate.transform.for_loop(step, np.arange(30000))
t = np.arange(30000) * brainstate.environ.get_dt()
# Time series (top row)
axes[0, idx].plot(t[-10000:], x[-10000:, 0], 'b-', linewidth=0.8)
axes[0, idx].set_title(f'$\mu$ = {mu_val}')
axes[0, idx].set_xlabel('Time (ms)')
if idx == 0:
axes[0, idx].set_ylabel('x(t)')
# Phase portrait (bottom row)
axes[1, idx].plot(x[-10000:, 0], y[-10000:, 0], 'k-', linewidth=0.3, alpha=0.7)
# Nullcline
x_null = np.linspace(-3, 3, 100)
axes[1, idx].plot(x_null, x_null - x_null ** 3 / 3, 'g--', linewidth=1.5, alpha=0.7)
axes[1, idx].set_xlabel('x')
if idx == 0:
axes[1, idx].set_ylabel('y')
axes[1, idx].set_xlim([-3, 3])
axes[1, idx].set_ylim([-3, 3])
plt.suptitle('Transition from Harmonic to Relaxation Oscillations', y=1.02, fontsize=14)
plt.tight_layout()
plt.show()
<>:26: SyntaxWarning: invalid escape sequence '\m'
<>:26: SyntaxWarning: invalid escape sequence '\m'
C:\Users\adadu\AppData\Local\Temp\ipykernel_32900\4055819095.py:26: SyntaxWarning: invalid escape sequence '\m'
axes[0, idx].set_title(f'$\mu$ = {mu_val}')
x, y = simulate(5.0, 50000)
t = np.arange(50000) * brainstate.environ.get_dt()
# Focus on steady state
x_ss, y_ss, t_ss = x[-20000:, 0], y[-20000:, 0], t[-20000:]
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
# Time series with annotations
axes[0].plot(t_ss, x_ss, 'b-', linewidth=1)
axes[0].set_xlabel('Time (ms)')
axes[0].set_ylabel('x(t)')
axes[0].set_title('Relaxation Oscillation ($\mu=5$)')
# Annotate phases
axes[0].annotate('Fast jump', xy=(t_ss[2000], 1.5), fontsize=10, color='red')
axes[0].annotate('Slow phase', xy=(t_ss[8000], 1.8), fontsize=10, color='green')
# Phase portrait showing slow-fast structure
axes[1].plot(x_ss, y_ss, 'b-', linewidth=1, alpha=0.7)
# Nullcline (cubic)
x_null = np.linspace(-2.5, 2.5, 200)
y_null = x_null - x_null ** 3 / 3
axes[1].plot(x_null, y_null, 'r-', linewidth=3, label='Slow manifold (x-nullcline)')
# Mark regions
axes[1].fill_between(x_null[x_null > 1], y_null[x_null > 1], 5, alpha=0.2, color='green', label='Slow (right branch)')
axes[1].fill_between(x_null[x_null < -1], y_null[x_null < -1], -5, alpha=0.2, color='green')
axes[1].set_xlabel('x')
axes[1].set_ylabel('y')
axes[1].set_title('Slow-Fast Dynamics in Phase Space')
axes[1].legend(loc='upper right')
axes[1].set_xlim([-3, 3])
axes[1].set_ylim([-2, 2])
plt.tight_layout()
plt.show()
<>:13: SyntaxWarning: invalid escape sequence '\m'
<>:13: SyntaxWarning: invalid escape sequence '\m'
C:\Users\adadu\AppData\Local\Temp\ipykernel_32900\2480255506.py:13: SyntaxWarning: invalid escape sequence '\m'
axes[0].set_title('Relaxation Oscillation ($\mu=5$)')
C:\Users\adadu\AppData\Local\Temp\ipykernel_32900\2480255506.py:13: SyntaxWarning: invalid escape sequence '\m'
axes[0].set_title('Relaxation Oscillation ($\mu=5$)')
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
File ~\miniconda3\envs\brainx\Lib\site-packages\matplotlib\artist.py:1402, in _get_tightbbox_for_layout_only(obj, *args, **kwargs)
1401 try:
-> 1402 return obj.get_tightbbox(*args, **{**kwargs, "for_layout_only": True})
1403 except TypeError:
File ~\miniconda3\envs\brainx\Lib\site-packages\matplotlib\axes\_base.py:4587, in _AxesBase.get_tightbbox(self, renderer, call_axes_locator, bbox_extra_artists, for_layout_only)
4586 for a in bbox_artists:
-> 4587 bbox = a.get_tightbbox(renderer)
4588 if (bbox is not None
4589 and 0 < bbox.width < np.inf
4590 and 0 < bbox.height < np.inf):
File ~\miniconda3\envs\brainx\Lib\site-packages\matplotlib\text.py:2030, in Annotation.get_tightbbox(self, renderer)
2028 def get_tightbbox(self, renderer=None):
2029 # docstring inherited
-> 2030 if not self._check_xy(renderer):
2031 return Bbox.null()
File ~\miniconda3\envs\brainx\Lib\site-packages\matplotlib\text.py:1590, in _AnnotationBase._check_xy(self, renderer)
1588 if b or (b is None and self.xycoords == "data"):
1589 # check if self.xy is inside the Axes.
-> 1590 xy_pixel = self._get_position_xy(renderer)
1591 return self.axes.contains_point(xy_pixel)
File ~\miniconda3\envs\brainx\Lib\site-packages\matplotlib\text.py:1581, in _AnnotationBase._get_position_xy(self, renderer)
1580 """Return the pixel position of the annotated point."""
-> 1581 return self._get_xy(renderer, self.xy, self.xycoords)
File ~\miniconda3\envs\brainx\Lib\site-packages\matplotlib\text.py:1475, in _AnnotationBase._get_xy(self, renderer, xy, coords)
1474 if xcoord == 'data':
-> 1475 x = float(self.convert_xunits(x))
1476 if ycoord == 'data':
File ~\miniconda3\envs\brainx\Lib\site-packages\matplotlib\artist.py:278, in Artist.convert_xunits(self, x)
277 return x
--> 278 return ax.xaxis.convert_units(x)
File ~\miniconda3\envs\brainx\Lib\site-packages\matplotlib\axis.py:1802, in Axis.convert_units(self, x)
1800 def convert_units(self, x):
1801 # If x is natively supported by Matplotlib, doesn't need converting
-> 1802 if munits._is_natively_supported(x):
1803 return x
File ~\miniconda3\envs\brainx\Lib\site-packages\matplotlib\units.py:64, in _is_natively_supported(x)
62 if np.iterable(x):
63 # Assume lists are homogeneous as other functions in unit system.
---> 64 for thisx in x:
65 if thisx is ma.masked:
File ~\miniconda3\envs\brainx\Lib\site-packages\saiunit\_base_quantity.py:1195, in Quantity.__iter__(self)
1194 if self.ndim == 0:
-> 1195 raise TypeError("iteration over a 0-d Quantity is not allowed")
1196 for i in range(self.shape[0]):
TypeError: iteration over a 0-d Quantity is not allowed
During handling of the above exception, another exception occurred:
TypeError Traceback (most recent call last)
Cell In[10], line 38
35 axes[1].set_xlim([-3, 3])
36 axes[1].set_ylim([-2, 2])
---> 38 plt.tight_layout()
39 plt.show()
File ~\miniconda3\envs\brainx\Lib\site-packages\matplotlib\pyplot.py:2843, in tight_layout(pad, h_pad, w_pad, rect)
2835 @_copy_docstring_and_deprecators(Figure.tight_layout)
2836 def tight_layout(
2837 *,
(...) 2841 rect: tuple[float, float, float, float] | None = None,
2842 ) -> None:
-> 2843 gcf().tight_layout(pad=pad, h_pad=h_pad, w_pad=w_pad, rect=rect)
File ~\miniconda3\envs\brainx\Lib\site-packages\matplotlib\figure.py:3640, in Figure.tight_layout(self, pad, h_pad, w_pad, rect)
3638 previous_engine = self.get_layout_engine()
3639 self.set_layout_engine(engine)
-> 3640 engine.execute(self)
3641 if previous_engine is not None and not isinstance(
3642 previous_engine, (TightLayoutEngine, PlaceHolderLayoutEngine)
3643 ):
3644 _api.warn_external('The figure layout has changed to tight')
File ~\miniconda3\envs\brainx\Lib\site-packages\matplotlib\layout_engine.py:188, in TightLayoutEngine.execute(self, fig)
186 renderer = fig._get_renderer()
187 with getattr(renderer, "_draw_disabled", nullcontext)():
--> 188 kwargs = get_tight_layout_figure(
189 fig, fig.axes, get_subplotspec_list(fig.axes), renderer,
190 pad=info['pad'], h_pad=info['h_pad'], w_pad=info['w_pad'],
191 rect=info['rect'])
192 if kwargs:
193 fig.subplots_adjust(**kwargs)
File ~\miniconda3\envs\brainx\Lib\site-packages\matplotlib\_tight_layout.py:266, in get_tight_layout_figure(fig, axes_list, subplotspec_list, renderer, pad, h_pad, w_pad, rect)
261 return {}
262 span_pairs.append((
263 slice(ss.rowspan.start * div_row, ss.rowspan.stop * div_row),
264 slice(ss.colspan.start * div_col, ss.colspan.stop * div_col)))
--> 266 kwargs = _auto_adjust_subplotpars(fig, renderer,
267 shape=(max_nrows, max_ncols),
268 span_pairs=span_pairs,
269 subplot_list=subplot_list,
270 ax_bbox_list=ax_bbox_list,
271 pad=pad, h_pad=h_pad, w_pad=w_pad)
273 # kwargs can be none if tight_layout fails...
274 if rect is not None and kwargs is not None:
275 # if rect is given, the whole subplots area (including
276 # labels) will fit into the rect instead of the
(...) 280 # auto_adjust_subplotpars twice, where the second run
281 # with adjusted rect parameters.
File ~\miniconda3\envs\brainx\Lib\site-packages\matplotlib\_tight_layout.py:82, in _auto_adjust_subplotpars(fig, renderer, shape, span_pairs, subplot_list, ax_bbox_list, pad, h_pad, w_pad, rect)
80 for ax in subplots:
81 if ax.get_visible():
---> 82 bb += [martist._get_tightbbox_for_layout_only(ax, renderer)]
84 tight_bbox_raw = Bbox.union(bb)
85 tight_bbox = fig.transFigure.inverted().transform_bbox(tight_bbox_raw)
File ~\miniconda3\envs\brainx\Lib\site-packages\matplotlib\artist.py:1404, in _get_tightbbox_for_layout_only(obj, *args, **kwargs)
1402 return obj.get_tightbbox(*args, **{**kwargs, "for_layout_only": True})
1403 except TypeError:
-> 1404 return obj.get_tightbbox(*args, **kwargs)
File ~\miniconda3\envs\brainx\Lib\site-packages\matplotlib\axes\_base.py:4587, in _AxesBase.get_tightbbox(self, renderer, call_axes_locator, bbox_extra_artists, for_layout_only)
4584 bbox_artists = self.get_default_bbox_extra_artists()
4586 for a in bbox_artists:
-> 4587 bbox = a.get_tightbbox(renderer)
4588 if (bbox is not None
4589 and 0 < bbox.width < np.inf
4590 and 0 < bbox.height < np.inf):
4591 bb.append(bbox)
File ~\miniconda3\envs\brainx\Lib\site-packages\matplotlib\text.py:2030, in Annotation.get_tightbbox(self, renderer)
2028 def get_tightbbox(self, renderer=None):
2029 # docstring inherited
-> 2030 if not self._check_xy(renderer):
2031 return Bbox.null()
2032 return super().get_tightbbox(renderer)
File ~\miniconda3\envs\brainx\Lib\site-packages\matplotlib\text.py:1590, in _AnnotationBase._check_xy(self, renderer)
1587 b = self.get_annotation_clip()
1588 if b or (b is None and self.xycoords == "data"):
1589 # check if self.xy is inside the Axes.
-> 1590 xy_pixel = self._get_position_xy(renderer)
1591 return self.axes.contains_point(xy_pixel)
1592 return True
File ~\miniconda3\envs\brainx\Lib\site-packages\matplotlib\text.py:1581, in _AnnotationBase._get_position_xy(self, renderer)
1579 def _get_position_xy(self, renderer):
1580 """Return the pixel position of the annotated point."""
-> 1581 return self._get_xy(renderer, self.xy, self.xycoords)
File ~\miniconda3\envs\brainx\Lib\site-packages\matplotlib\text.py:1475, in _AnnotationBase._get_xy(self, renderer, xy, coords)
1473 xcoord, ycoord = coords if isinstance(coords, tuple) else (coords, coords)
1474 if xcoord == 'data':
-> 1475 x = float(self.convert_xunits(x))
1476 if ycoord == 'data':
1477 y = float(self.convert_yunits(y))
File ~\miniconda3\envs\brainx\Lib\site-packages\matplotlib\artist.py:278, in Artist.convert_xunits(self, x)
276 if ax is None or ax.xaxis is None:
277 return x
--> 278 return ax.xaxis.convert_units(x)
File ~\miniconda3\envs\brainx\Lib\site-packages\matplotlib\axis.py:1802, in Axis.convert_units(self, x)
1800 def convert_units(self, x):
1801 # If x is natively supported by Matplotlib, doesn't need converting
-> 1802 if munits._is_natively_supported(x):
1803 return x
1805 if self._converter is None:
File ~\miniconda3\envs\brainx\Lib\site-packages\matplotlib\units.py:64, in _is_natively_supported(x)
61 # Matplotlib natively supports all number types except Decimal.
62 if np.iterable(x):
63 # Assume lists are homogeneous as other functions in unit system.
---> 64 for thisx in x:
65 if thisx is ma.masked:
66 continue
File ~\miniconda3\envs\brainx\Lib\site-packages\saiunit\_base_quantity.py:1195, in Quantity.__iter__(self)
1186 """Solve the issue of DeviceArray.__iter__.
1187
1188 Details please see JAX issues:
(...) 1191 - https://github.com/google/jax/pull/3821
1192 """
1194 if self.ndim == 0:
-> 1195 raise TypeError("iteration over a 0-d Quantity is not allowed")
1196 for i in range(self.shape[0]):
1197 yield Quantity(self.mantissa[i], unit=self.unit)
TypeError: iteration over a 0-d Quantity is not allowed
Error in callback <function _draw_all_if_interactive at 0x00000205F9EF1B20> (for post_execute), with arguments args (),kwargs {}:
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
File ~\miniconda3\envs\brainx\Lib\site-packages\matplotlib\pyplot.py:278, in _draw_all_if_interactive()
276 def _draw_all_if_interactive() -> None:
277 if matplotlib.is_interactive():
--> 278 draw_all()
File ~\miniconda3\envs\brainx\Lib\site-packages\matplotlib\_pylab_helpers.py:131, in Gcf.draw_all(cls, force)
129 for manager in cls.get_all_fig_managers():
130 if force or manager.canvas.figure.stale:
--> 131 manager.canvas.draw_idle()
File ~\miniconda3\envs\brainx\Lib\site-packages\matplotlib\backend_bases.py:1893, in FigureCanvasBase.draw_idle(self, *args, **kwargs)
1891 if not self._is_idle_drawing:
1892 with self._idle_draw_cntx():
-> 1893 self.draw(*args, **kwargs)
File ~\miniconda3\envs\brainx\Lib\site-packages\matplotlib\backends\backend_agg.py:382, in FigureCanvasAgg.draw(self)
379 # Acquire a lock on the shared font cache.
380 with (self.toolbar._wait_cursor_for_draw_cm() if self.toolbar
381 else nullcontext()):
--> 382 self.figure.draw(self.renderer)
383 # A GUI class may be need to update a window using this draw, so
384 # don't forget to call the superclass.
385 super().draw()
File ~\miniconda3\envs\brainx\Lib\site-packages\matplotlib\artist.py:94, in _finalize_rasterization.<locals>.draw_wrapper(artist, renderer, *args, **kwargs)
92 @wraps(draw)
93 def draw_wrapper(artist, renderer, *args, **kwargs):
---> 94 result = draw(artist, renderer, *args, **kwargs)
95 if renderer._rasterizing:
96 renderer.stop_rasterizing()
File ~\miniconda3\envs\brainx\Lib\site-packages\matplotlib\artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
68 if artist.get_agg_filter() is not None:
69 renderer.start_filter()
---> 71 return draw(artist, renderer)
72 finally:
73 if artist.get_agg_filter() is not None:
File ~\miniconda3\envs\brainx\Lib\site-packages\matplotlib\figure.py:3257, in Figure.draw(self, renderer)
3254 # ValueError can occur when resizing a window.
3256 self.patch.draw(renderer)
-> 3257 mimage._draw_list_compositing_images(
3258 renderer, self, artists, self.suppressComposite)
3260 renderer.close_group('figure')
3261 finally:
File ~\miniconda3\envs\brainx\Lib\site-packages\matplotlib\image.py:134, in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
132 if not_composite or not has_images:
133 for a in artists:
--> 134 a.draw(renderer)
135 else:
136 # Composite any adjacent images together
137 image_group = []
File ~\miniconda3\envs\brainx\Lib\site-packages\matplotlib\artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
68 if artist.get_agg_filter() is not None:
69 renderer.start_filter()
---> 71 return draw(artist, renderer)
72 finally:
73 if artist.get_agg_filter() is not None:
File ~\miniconda3\envs\brainx\Lib\site-packages\matplotlib\axes\_base.py:3226, in _AxesBase.draw(self, renderer)
3223 if artists_rasterized:
3224 _draw_rasterized(self.get_figure(root=True), artists_rasterized, renderer)
-> 3226 mimage._draw_list_compositing_images(
3227 renderer, self, artists, self.get_figure(root=True).suppressComposite)
3229 renderer.close_group('axes')
3230 self.stale = False
File ~\miniconda3\envs\brainx\Lib\site-packages\matplotlib\image.py:134, in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
132 if not_composite or not has_images:
133 for a in artists:
--> 134 a.draw(renderer)
135 else:
136 # Composite any adjacent images together
137 image_group = []
File ~\miniconda3\envs\brainx\Lib\site-packages\matplotlib\artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
68 if artist.get_agg_filter() is not None:
69 renderer.start_filter()
---> 71 return draw(artist, renderer)
72 finally:
73 if artist.get_agg_filter() is not None:
File ~\miniconda3\envs\brainx\Lib\site-packages\matplotlib\text.py:1990, in Annotation.draw(self, renderer)
1988 if renderer is not None:
1989 self._renderer = renderer
-> 1990 if not self.get_visible() or not self._check_xy(renderer):
1991 return
1992 # Update text positions before `Text.draw` would, so that the
1993 # FancyArrowPatch is correctly positioned.
File ~\miniconda3\envs\brainx\Lib\site-packages\matplotlib\text.py:1590, in _AnnotationBase._check_xy(self, renderer)
1587 b = self.get_annotation_clip()
1588 if b or (b is None and self.xycoords == "data"):
1589 # check if self.xy is inside the Axes.
-> 1590 xy_pixel = self._get_position_xy(renderer)
1591 return self.axes.contains_point(xy_pixel)
1592 return True
File ~\miniconda3\envs\brainx\Lib\site-packages\matplotlib\text.py:1581, in _AnnotationBase._get_position_xy(self, renderer)
1579 def _get_position_xy(self, renderer):
1580 """Return the pixel position of the annotated point."""
-> 1581 return self._get_xy(renderer, self.xy, self.xycoords)
File ~\miniconda3\envs\brainx\Lib\site-packages\matplotlib\text.py:1475, in _AnnotationBase._get_xy(self, renderer, xy, coords)
1473 xcoord, ycoord = coords if isinstance(coords, tuple) else (coords, coords)
1474 if xcoord == 'data':
-> 1475 x = float(self.convert_xunits(x))
1476 if ycoord == 'data':
1477 y = float(self.convert_yunits(y))
File ~\miniconda3\envs\brainx\Lib\site-packages\matplotlib\artist.py:278, in Artist.convert_xunits(self, x)
276 if ax is None or ax.xaxis is None:
277 return x
--> 278 return ax.xaxis.convert_units(x)
File ~\miniconda3\envs\brainx\Lib\site-packages\matplotlib\axis.py:1802, in Axis.convert_units(self, x)
1800 def convert_units(self, x):
1801 # If x is natively supported by Matplotlib, doesn't need converting
-> 1802 if munits._is_natively_supported(x):
1803 return x
1805 if self._converter is None:
File ~\miniconda3\envs\brainx\Lib\site-packages\matplotlib\units.py:64, in _is_natively_supported(x)
61 # Matplotlib natively supports all number types except Decimal.
62 if np.iterable(x):
63 # Assume lists are homogeneous as other functions in unit system.
---> 64 for thisx in x:
65 if thisx is ma.masked:
66 continue
File ~\miniconda3\envs\brainx\Lib\site-packages\saiunit\_base_quantity.py:1195, in Quantity.__iter__(self)
1186 """Solve the issue of DeviceArray.__iter__.
1187
1188 Details please see JAX issues:
(...) 1191 - https://github.com/google/jax/pull/3821
1192 """
1194 if self.ndim == 0:
-> 1195 raise TypeError("iteration over a 0-d Quantity is not allowed")
1196 for i in range(self.shape[0]):
1197 yield Quantity(self.mantissa[i], unit=self.unit)
TypeError: iteration over a 0-d Quantity is not allowed
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
File ~\miniconda3\envs\brainx\Lib\site-packages\IPython\core\formatters.py:402, in BaseFormatter.__call__(self, obj)
400 pass
401 else:
--> 402 return printer(obj)
403 # Finally look for special method names
404 method = get_real_method(obj, self.print_method)
File ~\miniconda3\envs\brainx\Lib\site-packages\IPython\core\pylabtools.py:170, in print_figure(fig, fmt, bbox_inches, base64, **kwargs)
167 from matplotlib.backend_bases import FigureCanvasBase
168 FigureCanvasBase(fig)
--> 170 fig.canvas.print_figure(bytes_io, **kw)
171 data = bytes_io.getvalue()
172 if fmt == 'svg':
File ~\miniconda3\envs\brainx\Lib\site-packages\matplotlib\backend_bases.py:2157, in FigureCanvasBase.print_figure(self, filename, dpi, facecolor, edgecolor, orientation, format, bbox_inches, pad_inches, bbox_extra_artists, backend, **kwargs)
2154 # we do this instead of `self.figure.draw_without_rendering`
2155 # so that we can inject the orientation
2156 with getattr(renderer, "_draw_disabled", nullcontext)():
-> 2157 self.figure.draw(renderer)
2158 if bbox_inches:
2159 if bbox_inches == "tight":
File ~\miniconda3\envs\brainx\Lib\site-packages\matplotlib\artist.py:94, in _finalize_rasterization.<locals>.draw_wrapper(artist, renderer, *args, **kwargs)
92 @wraps(draw)
93 def draw_wrapper(artist, renderer, *args, **kwargs):
---> 94 result = draw(artist, renderer, *args, **kwargs)
95 if renderer._rasterizing:
96 renderer.stop_rasterizing()
File ~\miniconda3\envs\brainx\Lib\site-packages\matplotlib\artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
68 if artist.get_agg_filter() is not None:
69 renderer.start_filter()
---> 71 return draw(artist, renderer)
72 finally:
73 if artist.get_agg_filter() is not None:
File ~\miniconda3\envs\brainx\Lib\site-packages\matplotlib\figure.py:3257, in Figure.draw(self, renderer)
3254 # ValueError can occur when resizing a window.
3256 self.patch.draw(renderer)
-> 3257 mimage._draw_list_compositing_images(
3258 renderer, self, artists, self.suppressComposite)
3260 renderer.close_group('figure')
3261 finally:
File ~\miniconda3\envs\brainx\Lib\site-packages\matplotlib\image.py:134, in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
132 if not_composite or not has_images:
133 for a in artists:
--> 134 a.draw(renderer)
135 else:
136 # Composite any adjacent images together
137 image_group = []
File ~\miniconda3\envs\brainx\Lib\site-packages\matplotlib\artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
68 if artist.get_agg_filter() is not None:
69 renderer.start_filter()
---> 71 return draw(artist, renderer)
72 finally:
73 if artist.get_agg_filter() is not None:
File ~\miniconda3\envs\brainx\Lib\site-packages\matplotlib\axes\_base.py:3226, in _AxesBase.draw(self, renderer)
3223 if artists_rasterized:
3224 _draw_rasterized(self.get_figure(root=True), artists_rasterized, renderer)
-> 3226 mimage._draw_list_compositing_images(
3227 renderer, self, artists, self.get_figure(root=True).suppressComposite)
3229 renderer.close_group('axes')
3230 self.stale = False
File ~\miniconda3\envs\brainx\Lib\site-packages\matplotlib\image.py:134, in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
132 if not_composite or not has_images:
133 for a in artists:
--> 134 a.draw(renderer)
135 else:
136 # Composite any adjacent images together
137 image_group = []
File ~\miniconda3\envs\brainx\Lib\site-packages\matplotlib\artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
68 if artist.get_agg_filter() is not None:
69 renderer.start_filter()
---> 71 return draw(artist, renderer)
72 finally:
73 if artist.get_agg_filter() is not None:
File ~\miniconda3\envs\brainx\Lib\site-packages\matplotlib\text.py:1990, in Annotation.draw(self, renderer)
1988 if renderer is not None:
1989 self._renderer = renderer
-> 1990 if not self.get_visible() or not self._check_xy(renderer):
1991 return
1992 # Update text positions before `Text.draw` would, so that the
1993 # FancyArrowPatch is correctly positioned.
File ~\miniconda3\envs\brainx\Lib\site-packages\matplotlib\text.py:1590, in _AnnotationBase._check_xy(self, renderer)
1587 b = self.get_annotation_clip()
1588 if b or (b is None and self.xycoords == "data"):
1589 # check if self.xy is inside the Axes.
-> 1590 xy_pixel = self._get_position_xy(renderer)
1591 return self.axes.contains_point(xy_pixel)
1592 return True
File ~\miniconda3\envs\brainx\Lib\site-packages\matplotlib\text.py:1581, in _AnnotationBase._get_position_xy(self, renderer)
1579 def _get_position_xy(self, renderer):
1580 """Return the pixel position of the annotated point."""
-> 1581 return self._get_xy(renderer, self.xy, self.xycoords)
File ~\miniconda3\envs\brainx\Lib\site-packages\matplotlib\text.py:1475, in _AnnotationBase._get_xy(self, renderer, xy, coords)
1473 xcoord, ycoord = coords if isinstance(coords, tuple) else (coords, coords)
1474 if xcoord == 'data':
-> 1475 x = float(self.convert_xunits(x))
1476 if ycoord == 'data':
1477 y = float(self.convert_yunits(y))
File ~\miniconda3\envs\brainx\Lib\site-packages\matplotlib\artist.py:278, in Artist.convert_xunits(self, x)
276 if ax is None or ax.xaxis is None:
277 return x
--> 278 return ax.xaxis.convert_units(x)
File ~\miniconda3\envs\brainx\Lib\site-packages\matplotlib\axis.py:1802, in Axis.convert_units(self, x)
1800 def convert_units(self, x):
1801 # If x is natively supported by Matplotlib, doesn't need converting
-> 1802 if munits._is_natively_supported(x):
1803 return x
1805 if self._converter is None:
File ~\miniconda3\envs\brainx\Lib\site-packages\matplotlib\units.py:64, in _is_natively_supported(x)
61 # Matplotlib natively supports all number types except Decimal.
62 if np.iterable(x):
63 # Assume lists are homogeneous as other functions in unit system.
---> 64 for thisx in x:
65 if thisx is ma.masked:
66 continue
File ~\miniconda3\envs\brainx\Lib\site-packages\saiunit\_base_quantity.py:1195, in Quantity.__iter__(self)
1186 """Solve the issue of DeviceArray.__iter__.
1187
1188 Details please see JAX issues:
(...) 1191 - https://github.com/google/jax/pull/3821
1192 """
1194 if self.ndim == 0:
-> 1195 raise TypeError("iteration over a 0-d Quantity is not allowed")
1196 for i in range(self.shape[0]):
1197 yield Quantity(self.mantissa[i], unit=self.unit)
TypeError: iteration over a 0-d Quantity is not allowed
<Figure size 1200x500 with 2 Axes>
# High mu relaxation oscillator
model = brainmass.VanDerPolStep(1, mu=5.0)
brainstate.nn.init_all_states(model)
model.x.value = np.array([0.5])
model.y.value = np.array([0.0])
def step(i):
with brainstate.environ.context(i=i, t=i * brainstate.environ.get_dt()):
model.update()
return model.x.value, model.y.value
x, y = brainstate.transform.for_loop(step, np.arange(50000))
t = np.arange(50000) * brainstate.environ.get_dt()
# Focus on steady state
x_ss, y_ss, t_ss = x[-20000:, 0], y[-20000:, 0], t[-20000:]
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
# Time series with annotations
axes[0].plot(t_ss, x_ss, 'b-', linewidth=1)
axes[0].set_xlabel('Time (ms)')
axes[0].set_ylabel('x(t)')
axes[0].set_title('Relaxation Oscillation ($\mu=5$)')
# Annotate phases
axes[0].annotate('Fast jump', xy=(t_ss[2000], 1.5), fontsize=10, color='red')
axes[0].annotate('Slow phase', xy=(t_ss[8000], 1.8), fontsize=10, color='green')
# Phase portrait showing slow-fast structure
axes[1].plot(x_ss, y_ss, 'b-', linewidth=1, alpha=0.7)
# Nullcline (cubic)
x_null = np.linspace(-2.5, 2.5, 200)
y_null = x_null - x_null ** 3 / 3
axes[1].plot(x_null, y_null, 'r-', linewidth=3, label='Slow manifold (x-nullcline)')
# Mark regions
axes[1].fill_between(x_null[x_null > 1], y_null[x_null > 1], 5, alpha=0.2, color='green', label='Slow (right branch)')
axes[1].fill_between(x_null[x_null < -1], y_null[x_null < -1], -5, alpha=0.2, color='green')
axes[1].set_xlabel('x')
axes[1].set_ylabel('y')
axes[1].set_title('Slow-Fast Dynamics in Phase Space')
axes[1].legend(loc='upper right')
axes[1].set_xlim([-3, 3])
axes[1].set_ylim([-2, 2])
plt.tight_layout()
plt.show()
# Measure periods for different mu
mu_range = np.array([0.5, 1.0, 2.0, 3.0, 4.0, 5.0])
periods = []
dt = brainstate.environ.get_dt() / u.ms
for mu_val in mu_range:
x, y = simulate(mu_val, 100000)
# Find period from zero crossings
x_ss = x[-50000:, 0]
crossings = np.where(np.diff(np.sign(x_ss)) > 0)[0]
if len(crossings) > 2:
period = np.mean(np.diff(crossings)) * dt
periods.append(period)
else:
periods.append(np.nan)
periods = np.array(periods)
# Plot period vs mu
plt.figure(figsize=(8, 5))
plt.plot(mu_range, periods, 'bo-', markersize=8, linewidth=2, label='Measured')
# Theoretical approximation for large mu: T ~ (3 - 2*log(2)) * mu ~ 1.614 * mu
mu_theory = np.linspace(0.5, 5, 50)
T_theory = (3 - 2 * np.log(2)) * mu_theory
plt.plot(mu_theory, T_theory, 'r--', linewidth=2, label=r'Theory: $T \approx 1.614 \mu$ (large $\mu$)')
plt.xlabel('$\mu$', fontsize=12)
plt.ylabel('Period (ms)', fontsize=12)
plt.title('Oscillation Period vs Nonlinearity Parameter')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# Measure periods for different mu
mu_range = np.array([0.5, 1.0, 2.0, 3.0, 4.0, 5.0])
periods = []
for mu_val in mu_range:
model = brainmass.VanDerPolStep(1, mu=mu_val)
brainstate.nn.init_all_states(model)
model.x.value = np.array([0.5])
def step(i):
with brainstate.environ.context(i=i, t=i * brainstate.environ.get_dt()):
model.update()
return model.x.value
x = brainstate.transform.for_loop(step, np.arange(100000))
# Find period from zero crossings
x_ss = x[-50000:, 0]
crossings = np.where(np.diff(np.sign(x_ss)) > 0)[0]
if len(crossings) > 2:
period = np.mean(np.diff(crossings)) * float(brainstate.environ.get_dt() / u.ms)
periods.append(period)
else:
periods.append(np.nan)
periods = np.array(periods)
# Plot period vs mu
plt.figure(figsize=(8, 5))
plt.plot(mu_range, periods, 'bo-', markersize=8, linewidth=2, label='Measured')
# Theoretical approximation for large mu: T ~ (3 - 2*log(2)) * mu ~ 1.614 * mu
mu_theory = np.linspace(0.5, 5, 50)
T_theory = (3 - 2 * np.log(2)) * mu_theory
plt.plot(mu_theory, T_theory, 'r--', linewidth=2, label=r'Theory: $T \approx 1.614 \mu$ (large $\mu$)')
plt.xlabel('$\mu$', fontsize=12)
plt.ylabel('Period (ms)', fontsize=12)
plt.title('Oscillation Period vs Nonlinearity Parameter')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
Summary#
In this tutorial, you learned:
Van der Pol equations: The Lienard form with nonlinear damping
Parameter \(\mu\): Controls transition from harmonic to relaxation oscillations
Relaxation oscillations: Characterized by slow phases and fast jumps
Slow-fast dynamics: The trajectory follows the nullcline during slow phases
Key Insights#
Van der Pol captures threshold dynamics relevant to neural excitability
Large \(\mu\) creates spike-like waveforms similar to action potentials
The model exhibits limit cycle oscillations for all \(\mu > 0\)
References#
van der Pol, B. (1926). On relaxation-oscillations. The London, Edinburgh, and Dublin Philosophical Magazine, 2(11), 978-992.
FitzHugh, R. (1961). Impulses and physiological states in theoretical models of nerve membrane. Biophysical Journal, 1(6), 445-466.
Strogatz, S. H. (2015). Nonlinear Dynamics and Chaos. Westview Press. (Chapter 7)