Mathematical Foundations for Oscillators¶
This document provides the mathematical foundations needed to analyze biological oscillators: phase space concepts, linear stability analysis, and eigenvalue interpretation.
Why One-Dimensional Systems Cannot Oscillate¶
Consider a single differential equation:
$$\frac{dx}{dt} = f(x)$$
The state of the system is a single value $x$ on a line. At any point, the derivative $dx/dt = f(x)$ determines which direction $x$ moves. If $f(x) > 0$, then $x$ increases; if $f(x) < 0$, then $x$ decreases; if $f(x) = 0$, we're at a fixed point.
To oscillate, $x(t)$ must return periodically to its starting value. But on a line, once $x$ moves away from a point, it cannot return without passing through a fixed point where the motion stops. The only possible behaviors are monotonic approach to a fixed point or divergence to infinity.
Oscillations require at least two dimensions, where trajectories can form closed loops without intersecting fixed points.
Linear Stability Analysis¶
Linear stability analysis determines the local behavior near fixed points by linearizing the nonlinear dynamics. This approach is justified by the Hartman-Grobman theorem: near a hyperbolic fixed point (one where no eigenvalue has zero real part), the nonlinear system is topologically equivalent to its linearization. The qualitative behavior—stability, instability, and oscillatory dynamics—is captured by the linear approximation.
The Procedure¶
For a system $\dot{\mathbf{x}} = \mathbf{f}(\mathbf{x})$ where $\mathbf{x} = (x_1, x_2, \ldots, x_n)$, we follow four steps. First, find fixed points by solving $\mathbf{f}(\mathbf{x}^*) = 0$. Second, linearize around the fixed point by computing the Jacobian matrix $J_{ij} = \frac{\partial f_i}{\partial x_j}\bigg|_{\mathbf{x}^*}$, which gives the linearized dynamics $d\boldsymbol{\delta}/dt \approx J \cdot \boldsymbol{\delta}$ for small perturbations $\boldsymbol{\delta} = \mathbf{x} - \mathbf{x}^*$. Third, compute eigenvalues $\lambda$ by solving $\det(J - \lambda I) = 0$. Fourth, interpret the eigenvalues $\lambda = \sigma + i\omega$, where $\sigma = \text{Re}(\lambda)$ gives the growth or decay rate, and $\omega = \text{Im}(\lambda)$ gives the oscillation frequency with period $T = 2\pi/|\omega|$ when $\omega \neq 0$.
What This Analysis Reveals¶
Linear stability analysis identifies Hopf bifurcations, points where a fixed point changes stability as parameters vary and a limit cycle is born. When a system has complex eigenvalues ($\omega \neq 0$), trajectories spiral near the fixed point rather than approaching along straight lines. When the real part is positive ($\sigma > 0$), the fixed point is unstable and trajectories spiral outward. In nonlinear systems, saturation effects prevent unbounded growth and create a stable limit cycle at finite amplitude. This is the typical route to oscillations in biological feedback systems.
Two-Dimensional Systems¶
Computing Eigenvalues¶
For a $2 \times 2$ matrix $J = \begin{pmatrix} a & b \\ c & d \end{pmatrix}$, the determinant is $\det(J) = ad - bc$. The trace and determinant satisfy $\text{trace}(J) = a + d = \lambda_1 + \lambda_2$ and $\det(J) = \lambda_1 \cdot \lambda_2$. The characteristic equation $\lambda^2 - \text{trace}(J) \cdot \lambda + \det(J) = 0$ has solutions:
$$\lambda = \frac{\text{trace}(J) \pm \sqrt{\text{trace}(J)^2 - 4\det(J)}}{2}$$
The discriminant $\Delta = \text{trace}(J)^2 - 4\det(J)$ determines whether eigenvalues are real ($\Delta > 0$) or complex ($\Delta < 0$). Complex eigenvalues produce oscillatory behavior. Stability requires $\text{trace}(J) < 0$ and $\det(J) > 0$.
Classification of Fixed Points¶
The sign of the trace, determinant, and discriminant completely characterize the local dynamics:
| Trace | Det | $\Delta$ | Type | Behavior |
|---|---|---|---|---|
| $< 0$ | $> 0$ | $> 0$ | Stable node | Monotonic decay |
| $< 0$ | $> 0$ | $< 0$ | Stable spiral | Damped oscillations |
| $> 0$ | $> 0$ | $< 0$ | Unstable spiral | Growing oscillations |
| $> 0$ | $> 0$ | $> 0$ | Unstable node | Monotonic growth |
| any | $< 0$ | $> 0$ | Saddle | Unstable |
Example: Linear Oscillator¶
Consider $\frac{dx}{dt} = -y$ and $\frac{dy}{dt} = x$. The Jacobian $J = \begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix}$ has characteristic equation $\lambda^2 + 1 = 0$, giving $\lambda = \pm i$. The general solution is:
$$x(t) = A\cos(t) + B\sin(t), \quad y(t) = A\sin(t) - B\cos(t)$$
Trajectories are circles in phase space with period $T = 2\pi$. Purely imaginary eigenvalues produce neutral oscillations in linear systems. In nonlinear biological systems, we typically need $\sigma > 0$ with nonlinear saturation to create stable limit cycles.
Three-Dimensional Systems¶
Computing Determinants for 3×3 Matrices¶
For a $3 \times 3$ matrix $J = \begin{pmatrix} a & b & c \\ d & e & f \\ g & h & i \end{pmatrix}$, the determinant can be computed by cofactor expansion along the first row:
$$\det(J) = a(ei - fh) - b(di - fg) + c(dh - eg)$$
An alternative for $3 \times 3$ matrices only is the Rule of Sarrus: write the first two columns to the right, sum products of down-right diagonals, subtract products of up-right diagonals:
$$\det(J) = (aei + bfg + cdh) - (ceg + afh + bdi)$$
The characteristic equation for a $3 \times 3$ matrix is cubic and has either three real roots or one real root with two complex conjugates. The complex pair can produce oscillatory behavior if its real part is positive.
Dynamics of Complex Eigenvalues¶
Solving the Linearized System¶
For the linearized system $d\boldsymbol{\delta}/dt = J\boldsymbol{\delta}$, the formal solution involves the matrix exponential $\boldsymbol{\delta}(t) = e^{Jt} \boldsymbol{\delta}(0)$. When $J$ has $n$ linearly independent eigenvectors (true for matrices with distinct eigenvalues and many other cases), we can diagonalize as $J = V \Lambda V^{-1}$ where $V$ contains eigenvectors and $\Lambda$ is diagonal with eigenvalues. Then $e^{Jt} = V e^{\Lambda t} V^{-1}$, giving:
$$\boldsymbol{\delta}(t) = c_1 e^{\lambda_1 t} \mathbf{v}_1 + c_2 e^{\lambda_2 t} \mathbf{v}_2 + \cdots + c_n e^{\lambda_n t} \mathbf{v}_n$$
The coefficients $c_i$ are determined from initial conditions. Complex eigenvalues come in conjugate pairs with conjugate eigenvectors, ensuring real-valued solutions. (Defective matrices lacking a full set of independent eigenvectors require Jordan normal form, but typical biological system Jacobians are diagonalizable.)
From Eigenvalues to Oscillations¶
For complex eigenvalues $\lambda = \sigma \pm i\omega$, Euler's formula gives $e^{(\sigma + i\omega)t} = e^{\sigma t}[\cos(\omega t) + i\sin(\omega t)]$. After taking appropriate linear combinations to ensure real solutions, a single component behaves as:
$$x(t) - x^* = Ce^{\sigma t}\cos(\omega t + \phi)$$
Here's the key insight: $e^{\sigma t}$ controls amplitude growth or decay, while $\cos(\omega t + \phi)$ creates oscillation.
Three Regimes¶
The sign of $\sigma = \text{Re}(\lambda)$ determines long-term behavior.
When $\sigma < 0$ (damped oscillations), amplitude decays exponentially as $e^{\sigma t}$. Trajectories spiral inward to the fixed point with decay timescale $\tau = 1/|\sigma|$. After about $5\tau$, oscillations have essentially vanished.
When $\sigma > 0$ (growing oscillations), amplitude grows exponentially in the linear approximation. The fixed point is unstable and trajectories spiral outward. But here's what's interesting: in real nonlinear systems, as amplitude becomes large, nonlinear terms become important and saturate the growth, creating a stable limit cycle at finite amplitude. This is how biological oscillators work.
When $\sigma = 0$ (neutral oscillations), amplitude remains constant and trajectories form closed orbits. This is structurally unstable—any parameter change makes $\sigma \neq 0$.
Visualizing Phase Space Trajectories¶
Figure 1 shows phase space trajectories for a two-dimensional system with complex eigenvalues $\lambda = \sigma \pm i\omega$.
Panel (a) shows a linear system with damped oscillations: $\frac{dx}{dt} = -0.15x - y, \quad \frac{dy}{dt} = x - 0.15y$
The Jacobian is $J = \begin{pmatrix} -0.15 & -1 \\ 1 & -0.15 \end{pmatrix}$ with eigenvalues $\lambda = -0.15 \pm i$, giving $\sigma = -0.15$ and $\omega = 1$.
Panel (b) shows a nonlinear system that produces a limit cycle: $\frac{dx}{dt} = x(1 - x^2 - y^2) - y, \quad \frac{dy}{dt} = y(1 - x^2 - y^2) + x$
Near the origin, this linearizes to $J = \begin{pmatrix} 1 & -1 \\ 1 & 1 \end{pmatrix}$ with eigenvalues $\lambda = 1 \pm i$ (unstable spiral). The nonlinear terms $(1 - x^2 - y^2)$ create a stable limit cycle at radius $r = 1$.
"""
Figure 1: Phase space trajectories showing (a) damped oscillations and (b) limit cycle
Using ODE integration with solve_ivp
"""
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
# Panel (a): Damped oscillations - Linear system
# dx/dt = -0.15x - y, dy/dt = x - 0.15y
ax = axes[0]
def damped_system(t, state):
x, y = state
dxdt = -0.15*x - y
dydt = x - 0.15*y
return [dxdt, dydt]
# Solve ODE
t_span = (0, 40)
t_eval = np.linspace(0, 40, 1000)
initial_state = [1.5, 0]
sol = solve_ivp(damped_system, t_span, initial_state, t_eval=t_eval, method='RK45')
x, y = sol.y
ax.plot(x, y, 'b-', linewidth=2, label='Trajectory')
ax.plot(0, 0, 'ko', markersize=8, label='Fixed point')
# Add arrow to show direction
idx = 600
ax.arrow(x[idx], y[idx], x[idx+1]-x[idx], y[idx+1]-y[idx],
head_width=0.1, head_length=0.08, fc='blue', ec='blue')
ax.set_xlim(-2, 2)
ax.set_ylim(-2, 2)
ax.set_xlabel('$x$', fontsize=14)
ax.set_ylabel('$y$', fontsize=14)
ax.set_title(r'(a) $\sigma < 0$: Damped oscillations', fontsize=12)
ax.text(0, -2.3, 'Stable fixed point', ha='center', fontsize=10)
ax.grid(True, alpha=0.3)
ax.set_aspect('equal')
ax.legend(loc='upper right', fontsize=9)
# Panel (b): Limit cycle - Nonlinear system
# dx/dt = x(1 - x^2 - y^2) - y, dy/dt = y(1 - x^2 - y^2) + x
ax = axes[1]
def limit_cycle_system(t, state):
x, y = state
r_squared = x**2 + y**2
dxdt = x*(1 - r_squared) - y
dydt = y*(1 - r_squared) + x
return [dxdt, dydt]
# Draw limit cycle (analytical: r=1)
theta = np.linspace(0, 2*np.pi, 100)
r_lc = 1.0
x_lc = r_lc * np.cos(theta)
y_lc = r_lc * np.sin(theta)
ax.plot(x_lc, y_lc, 'r-', linewidth=3, label='Limit cycle')
# Add arrow on limit cycle
ax.arrow(r_lc, 0, -0.05, 0.15, head_width=0.1, head_length=0.08,
fc='red', ec='red')
# Solve ODE from inside the limit cycle
t_span = (0, 15)
t_eval = np.linspace(0, 15, 1000)
initial_state = [0.3, 0]
sol = solve_ivp(limit_cycle_system, t_span, initial_state, t_eval=t_eval, method='RK45')
x, y = sol.y
ax.plot(x, y, 'b-', linewidth=2, alpha=0.7, label='Trajectory')
# Add arrow to show direction
idx = 700
ax.arrow(x[idx], y[idx], x[idx+1]-x[idx], y[idx+1]-y[idx],
head_width=0.1, head_length=0.08, fc='blue', ec='blue')
ax.plot(0, 0, 'ko', markersize=8, label='Fixed point')
ax.set_xlim(-2, 2)
ax.set_ylim(-2, 2)
ax.set_xlabel('$x$', fontsize=14)
ax.set_ylabel('$y$', fontsize=14)
ax.set_title(r'(b) $\sigma > 0$: Limit cycle emerges', fontsize=12)
ax.text(0, -2.3, 'Unstable fixed point', ha='center', fontsize=10)
ax.grid(True, alpha=0.3)
ax.set_aspect('equal')
ax.legend(loc='upper right', fontsize=9)
plt.tight_layout()
plt.show()
In Figure 1(a), each loop of the spiral has smaller radius than the previous one. The spacing between loops depends on $|\omega/\sigma|$: large ratios produce tightly wound spirals with many oscillations before reaching the fixed point. In Figure 1(b), trajectories (blue) spiral outward from the unstable fixed point until nonlinear saturation halts the growth, creating a stable limit cycle (red)—a closed periodic orbit that nearby trajectories approach.
Figure 2 shows the evolution through Hopf bifurcation as the parameter $\sigma$ varies. The system is: $\frac{dx}{dt} = \sigma x - y - x(x^2 + y^2), \quad \frac{dy}{dt} = x + \sigma y - y(x^2 + y^2)$
This is the normal form of a supercritical Hopf bifurcation. Near the origin, the Jacobian is $J = \begin{pmatrix} \sigma & -1 \\ 1 & \sigma \end{pmatrix}$ with eigenvalues $\lambda = \sigma \pm i$. As $\sigma$ increases through zero, the system transitions from stable fixed point (damped oscillations) to unstable fixed point with a stable limit cycle. The limit cycle has radius $r = \sqrt{\sigma}$ for $\sigma > 0$.
"""
Figure 2: Evolution through Hopf bifurcation as sigma increases
Using ODE integration to show system behavior across the bifurcation
System: dx/dt = sigma*x - y - x(x^2 + y^2), dy/dt = x + sigma*y - y(x^2 + y^2)
"""
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
fig, axes = plt.subplots(1, 4, figsize=(16, 4))
configs = [
(-0.3, r'$\sigma < 0$'),
(-0.05, r'$\sigma \approx 0$'),
(0.05, r'$\sigma > 0$'),
(0.2, r'$\sigma \gg 0$')
]
def hopf_system(t, state, sigma):
"""Hopf bifurcation normal form"""
x, y = state
r_squared = x**2 + y**2
dxdt = sigma*x - y - x*r_squared
dydt = x + sigma*y - y*r_squared
return [dxdt, dydt]
for idx, (sigma, label) in enumerate(configs):
ax = axes[idx]
# Integrate ODE for trajectory
t_span = (0, 30)
t_eval = np.linspace(0, 30, 1000)
initial_state = [1.2, 0]
sol = solve_ivp(lambda t, y: hopf_system(t, y, sigma),
t_span, initial_state, t_eval=t_eval,
method='RK45', rtol=1e-8, atol=1e-10)
x, y = sol.y
ax.plot(x, y, 'b-', linewidth=2, alpha=0.8)
# For positive sigma, draw the theoretical limit cycle (r = sqrt(sigma))
if sigma > 0:
theta = np.linspace(0, 2*np.pi, 100)
r_lc = np.sqrt(sigma)
x_lc = r_lc * np.cos(theta)
y_lc = r_lc * np.sin(theta)
ax.plot(x_lc, y_lc, 'r--', linewidth=2.5, label=f'Limit cycle\n$r=\\sqrt{{\\sigma}}$', alpha=0.7)
ax.legend(loc='upper right', fontsize=8)
# Mark fixed point
ax.plot(0, 0, 'ko', markersize=6)
ax.set_xlim(-1.5, 1.5)
ax.set_ylim(-1.5, 1.5)
ax.set_xlabel('$x$', fontsize=11)
if idx == 0:
ax.set_ylabel('$y$', fontsize=11)
ax.set_title(label, fontsize=11)
ax.grid(True, alpha=0.3)
ax.set_aspect('equal')
ax.tick_params(labelsize=8)
plt.tight_layout()
plt.show()
Oscillation Quality¶
The ratio $|\omega/(2\sigma)|$ measures oscillation quality. When this ratio is much larger than 1, the system completes many oscillations during one exponential timescale, producing clear oscillatory behavior. When the ratio is near 1, only a few oscillations occur before amplitude changes substantially.
Higher Dimensions¶
In three dimensions with one real eigenvalue and complex conjugates, motion combines contraction or expansion along one direction with spiraling in the orthogonal plane, creating a 3D spiral that can form a limit cycle winding through phase space.
The Poincaré-Bendixson Theorem¶
For two-dimensional systems, this powerful theorem states: if a trajectory is bounded and doesn't approach a fixed point, it must approach a limit cycle. Combined with an unstable fixed point having $\sigma > 0$ and $\omega \neq 0$, plus bounded dynamics from nonlinearity, this guarantees oscillations emerge. In three or more dimensions, the picture is richer and chaos becomes possible, but for biological systems with appropriate feedback structure and dissipation, limit cycles typically emerge via Hopf bifurcation.
Summary¶
Linear stability analysis provides a systematic way to understand how oscillations arise in dynamical systems.
One-dimensional systems cannot oscillate because trajectories on a line cannot form closed loops without passing through fixed points. At least two dimensions are required.
Complex eigenvalues $\lambda = \sigma \pm i\omega$ indicate spiraling dynamics near the fixed point. The imaginary part $\omega$ determines oscillation frequency. The real part $\sigma$ determines stability: negative means trajectories spiral inward (damped oscillations), positive means trajectories spiral outward.
As parameters vary and $\sigma$ crosses zero from negative to positive, a Hopf bifurcation typically occurs and a limit cycle is born. This is the most common route to oscillations in biological systems. Nonlinearity (from Hill functions, resource limits, etc.) prevents unbounded growth and creates stable limit cycles at finite amplitude.
License: © 2025 Matthias Függer and Thomas Nowak. Licensed under CC BY-NC-SA 4.0.