Mathematical ECG and Heart Models

The heart is a nonlinear oscillator. Its rhythm depends on autonomous pacemaker cells that fire action potentials spontaneously, coupled electrically to the rest of the myocardium. The resulting wave of electrical depolarisation — recorded at the skin surface as the electrocardiogram (ECG) — encodes the health of the conduction system in its waveform. Remarkably, the essential dynamics of this system can be captured by the Van der Pol oscillator, a classical model of self-sustained oscillations, extended and coupled to reproduce the PQRST morphology of a realistic ECG trace.

1. Cardiac electrical activity

The heart's mechanical contraction is initiated and coordinated by electrical signals. Unlike skeletal muscle, which requires direct neural stimulation for each contraction, cardiac muscle is autorhythmic: specialised cells in the sinoatrial (SA) node spontaneously generate action potentials at a natural rate of 60–100 beats per minute in a resting adult.

This electrical signal propagates through the atria, inducing atrial contraction and converging at the atrioventricular (AV) node, where it is deliberately delayed (~120 ms) to allow the ventricles time to fill. The signal then travels rapidly through the His-Purkinje system — a specialised network of fibres with high conduction velocity — to reach ventricular myocardium simultaneously from both bottom corners, producing a coordinated squeeze.

Hierarchical pacemakers: The SA node is the primary pacemaker, but subsidiary pacemakers exist at all levels of the conduction system. If the SA node fails, the AV node takes over at 40–60 bpm (junctional rhythm). If the AV node also fails, ventricular cells pace at 20–40 bpm (idioventricular rhythm). This hierarchy provides redundant protection against cardiac arrest.

The electrocardiogram (ECG) records the net electrical dipole of the heart as seen from different surface electrode positions. Each deflection in the ECG trace corresponds to a specific event in the conduction sequence — the P wave (atrial depolarisation), QRS complex (ventricular depolarisation), and T wave (ventricular repolarisation).

2. Cardiac action potential

Cardiac action potentials are significantly different from neuronal action potentials. Ventricular myocytes have a characteristic plateau phase — a prolonged period of maintained depolarisation that lasts 200–400 ms, compared to the ~1 ms neuronal spike. This long duration prevents tetanic contraction (sustained contraction that would prevent the heart from refilling) and coincides with the absolute refractory period, protecting the heart from premature re-excitation.

The five phases of the ventricular action potential are:

Phase	Duration	Current	Voltage change
0 — Upstroke	1–2 ms	Fast Na⁺ inward (I_Na)	−90 → +30 mV
1 — Early repol.	10 ms	Transient K⁺ outward (I_to)	+30 → +10 mV
2 — Plateau	150–200 ms	L-type Ca²⁺ in vs slow K⁺ out	~0 mV
3 — Repolarisation	100 ms	Rapid K⁺ outward (I_Kr, I_Ks)	0 → −90 mV
4 — Rest	—	I_K1 (inward rectifier)	−90 mV (stable)

SA node cells differ: they lack the fast Na⁺ current responsible for the rapid upstroke. Their spontaneous depolarisation (phase 4) is driven by the funny current I_f — a hyperpolarisation-activated mixed cation current — causing a slow pacemaker potential that eventually triggers another action potential. This is inherently oscillatory behaviour.

3. Conduction system: SA, AV, His-Purkinje

The cardiac conduction system is a specialised network that generates and propagates electrical impulses with precise timing. From a dynamical systems perspective, it is a network of coupled oscillators with different intrinsic frequencies, in which the fastest oscillator (SA node, 60–100 bpm) normally enslaves the slower nodes through entrainment.

SA node (sinoatrial): upper right atrium; primary pacemaker; 60–100 bpm intrinsic rate; driven by I_f (HCN channels). Responds to autonomic input: sympathetic stimulation (adrenaline) increases rate; vagal stimulation (acetylcholine) decreases rate.
AV node (atrioventricular): boundary between atria and ventricles; introduces ~120–200 ms delay (visible as PR interval on ECG); intrinsic rate 40–60 bpm if not driven by SA.
Bundle of His: penetrates the fibrous skeleton; only electrical pathway from atria to ventricles; divides into left and right bundle branches.
Purkinje fibres: rapidly conduct to ventricular apex and walls; conduction velocity ~4 m/s (vs. ~0.05 m/s in nodal tissue); intrinsic rate 20–40 bpm.

Conduction velocities: The variation in conduction velocity is functionally critical. The slow AV nodal conduction (0.05 m/s) provides the PR delay. The fast Purkinje system (4 m/s) achieves near-simultaneous ventricular activation. Disruption of any part produces characteristic ECG abnormalities: AV block (prolonged PR or dropped beats), bundle branch block (widened QRS), and various arrhythmias.

4. The Van der Pol oscillator

The Van der Pol oscillator, introduced by Balthasar van der Pol in 1926 to model oscillations in vacuum tube circuits, has become a canonical model for biological oscillators including the heart. The equation is:

d²x/dt² - μ(1 - x²) dx/dt + x = 0 or as a 2D system: dx/dt = y dy/dt = μ(1 - x²)y - x where μ > 0 controls the strength of nonlinear damping.

The key feature is the nonlinear damping term −μ(1−x²)ẋ. For |x| < 1, this acts as negative damping (energy injection), amplifying small oscillations. For |x| > 1, it is positive damping (energy dissipation), limiting large oscillations. The result is a stable limit cycle: all initial conditions converge to the same periodic orbit.

For large μ (strongly nonlinear), the oscillation develops a fast-slow character: the trajectory spends most time near the two slow branches of the cubic nullcline, with rapid jumps between them. This produces a relaxation oscillation — precisely the shape of a cardiac action potential: slow phase 4 diastolic depolarisation, then a rapid spike, then slow plateau, then rapid repolarisation.

History: Van der Pol himself used his oscillator to model the heartbeat as early as 1928, noting that irregular heartbeats occur when a pacemaker is forced at a frequency near a sub-harmonic of its natural frequency — an early description of cardiac arrhythmias using nonlinear dynamics.

5. Dynamical ECG model

McSharry et al. (2003, IEEE Transactions on Biomedical Engineering) developed an influential dynamical model that generates realistic ECG waveforms using a three-dimensional system of ODEs. The model uses a trajectory on a limit cycle in 3D space to track the cardiac state, and adds Gaussian-shaped deflections to produce the PQRST morphology.

The model evolves on the unit circle in the (x, y) plane (representing the underlying pacemaker cycle) and generates the ECG signal as the z-coordinate:

dx/dt = α \cdot x - ω \cdot y dy/dt = α \cdot y + ω \cdot x dz/dt = -Σ_i a_i \cdot Δθ_i \cdot exp(-Δθ_i² / (2b_i²)) - (z - z₀) where α = 1 - \sqrt(x² + y²) (keeps trajectory on unit circle) ω = 2πf (angular frequency of heart rate f) Δθ_i = (θ - θ_i) mod 2π (angular distance to each ECG event) a_i, b_i, θ_i (amplitude, width, angle of each peak)

Five Gaussian terms i ∈ {P, Q, R, S, T} represent the five major ECG deflections. The angular positions θ_i determine the timing of each feature within the cardiac cycle. As the state point θ = atan2(y, x) sweeps around the circle at rate ω, it encounters each Gaussian in sequence, producing the characteristic PQRST waveform in z(t).

Event	θ_i (radians)	a_i (amplitude)	b_i (width)	Physiological event
P	−π/3	1.2	0.25	Atrial depolarisation
Q	−π/12	−5.0	0.10	Septal depolarisation
R	0	30.0	0.10	Ventricular depolarisation peak
S	π/12	−7.5	0.10	Basal ventricular depolarisation
T	π/2	0.75	0.40	Ventricular repolarisation

6. QRS complex and PQRST morphology

The QRS complex is the dominant ECG feature: a sharp, narrow (80–120 ms) deflection representing ventricular depolarisation. Its large amplitude reflects the simultaneous activation of the thick ventricular wall. The individual components are:

Q wave: small initial negative deflection; represents septal depolarisation travelling left-to-right (away from typical left-sided leads).
R wave: large positive peak; main left ventricular depolarisation wave travelling toward left-sided leads.
S wave: terminal negative deflection; basal depolarisation moving superiorly and rightward.

Surrounding the QRS, the ECG shows:

P wave: small, rounded, positive; atrial depolarisation; duration ~100 ms.
PR interval: from P onset to QRS onset; corresponds to AV conduction delay; normal 120–200 ms.
ST segment: isoelectric segment between S and T; all ventricular cells depolarised; elevation indicates ischaemia or infarction.
T wave: broad, positive; ventricular repolarisation; occurs in opposite direction to depolarisation wave due to epicardial cells repolarising first.
QT interval: QRS onset to T wave end; represents total ventricular electrical systole; prolongation (QTc > 440 ms) risks torsades de pointes arrhythmia.

7. Arrhythmia mechanisms

Arrhythmias — abnormal heart rhythms — arise from disturbances in impulse generation, conduction, or both. From a dynamical systems perspective, arrhythmias represent bifurcations, instabilities, or pathological attractors of the cardiac conduction system.

Enhanced automaticity and triggered activity

Enhanced automaticity occurs when non-pacemaker cells develop spontaneous depolarisation due to elevated intracellular Ca²⁺, hypokalaemia, or catecholamine excess. Additional ectopic foci compete with the SA node, producing premature beats. In the Van der Pol framework, this corresponds to increasing μ or adding a drive term to normally quiescent cells.

Triggered activity involves afterdepolarisations — membrane oscillations that can reach threshold and trigger extra spikes. Early afterdepolarisations (EADs) occur during the plateau phase (phase 2 or 3), prolonged by QT-prolonging drugs or hypokalaemia, and can initiate torsades de pointes. Delayed afterdepolarisations (DADs) occur in phase 4, driven by Ca²⁺ overload.

Re-entry circuits

Re-entry is the most common mechanism underlying sustained tachyarrhythmias including atrial flutter, atrial fibrillation, ventricular tachycardia, and ventricular fibrillation. It requires:

An anatomical or functional circuit.
Unidirectional block in one limb.
Slowed conduction allowing time for recovery at the block site.

Mathematically, re-entry is a rotating wave (spiral wave) in the cardiac excitable medium. The spiral tip traces a meander path; when the tip reaches a singularity it can break up into multiple wavelets — this transition represents the onset of fibrillation.

Rotor hypothesis of fibrillation: Modern theory posits that atrial and ventricular fibrillation are maintained by one or a few high-frequency rotors (spiral waves) whose periodic activity is fragmented into the apparently random multiple-wavelet pattern. Catheter ablation targeting rotor cores has been explored as a treatment for persistent atrial fibrillation.

8. JavaScript ECG simulation

The simulation below implements the McSharry dynamical ECG model using a 4th-order Runge-Kutta integrator. The scrolling trace shows the synthetic ECG signal (z-coordinate). Adjust the heart rate and the amplitude of the QRS complex to see how the ECG morphology changes. Toggle the arrhythmia mode to simulate irregular RR intervals.

Heart rate (bpm) 75 bpm QRS amplitude 30 PR interval (ms) 160 ms

The ECG trace is generated by numerically solving three coupled ODEs using Euler integration (sufficient for visualisation). The angular frequency ω = 2π · (HR / 60) rad/s drives the cardiac phase variable. Each time the phase variable crosses the angle of each PQRST Gaussian, that deflection contributes to the z-coordinate.

// McSharry dynamical ECG model
const events = [
  { name:'P', theta:-Math.PI/3,  a:1.2,  b:0.25 },
  { name:'Q', theta:-Math.PI/12, a:-5.0, b:0.10 },
  { name:'R', theta:0,           a:30.0, b:0.10 },
  { name:'S', theta:Math.PI/12,  a:-7.5, b:0.10 },
  { name:'T', theta:Math.PI/2,   a:0.75, b:0.40 },
];

function ecgStep(state, omega, dt) {
  const { x, y, z } = state;
  const theta = Math.atan2(y, x);
  const alpha = 1.0 - Math.sqrt(x * x + y * y);

  let dz = -(z - 0);  // baseline restoration
  for (const ev of events) {
    let dth = theta - ev.theta;
    // wrap to [-π, π]
    while (dth >  Math.PI) dth -= 2 * Math.PI;
    while (dth < -Math.PI) dth += 2 * Math.PI;
    dz -= ev.a * dth * Math.exp(-dth * dth / (2 * ev.b * ev.b));
  }

  return {
    x: x + dt * (alpha * x - omega * y),
    y: y + dt * (alpha * y + omega * x),
    z: z + dt * dz,
  };
}