Spotlight #35 – Neuroscience & Neural Circuits: Hodgkin-Huxley, Synaptic Plasticity and Whole-Brain Models

Neuroscience is unusual among the natural sciences in that it spans more than ten orders of magnitude in spatial scale: from synaptic vesicles (50 nm) to white matter tracts (metres in length), from the millisecond timescale of an action potential to the years of memory formation. Computational approaches are essential for integrating phenomena across these scales. The Hodgkin-Huxley model, the Ising model of neuroscience, has been the template for everything from single compartment models to multi-million neuron network simulations.

1. The Hodgkin-Huxley Model

Alan Hodgkin and Andrew Huxley performed their landmark experiments on the squid giant axon between 1947 and 1952, using voltage clamp to isolate the sodium and potassium conductances. Their four-equation model, published in the Journal of Physiology in 1952 and rewarded with the 1963 Nobel Prize, quantitatively reproduced the action potential for the first time.

Membrane current balance:
  C_m dV/dt = I_ext − I_Na − I_K − I_L

Ion currents:
  I_Na = g_Na m³h (V − E_Na)    E_Na = +50 mV
  I_K  = g_K  n&sup4;  (V − E_K )    E_K  = −77 mV
  I_L  = g_L       (V − E_L )    E_L  = −54 mV

Maximal conductances (squid axon):
  g_Na = 120 mS/cm²  g_K = 36 mS/cm²  g_L = 0.3 mS/cm²

Gating variable ODEs (x ∈ {m, h, n}):
  dx/dt = α_x(V)(1−x) − β_x(V) x
  x_∞ = α_x / (α_x + β_x),  τ_x = 1/(α_x + β_x)

Example rate functions (at 6.3°C):
  α_m = 0.1(V+40) / (1 − exp(−(V+40)/10))
  β_m = 4 exp(−(V+65)/18)
  α_h = 0.07 exp(−(V+65)/20)
  β_h = 1 / (1 + exp(−(V+35)/10))
          

The model generates action potentials via a positive feedback loop: depolarisation opens Na⁺ channels (fast, activation gate m), which depolarise the membrane further, opening more channels. This regenerative process is terminated by Na⁺ channel inactivation (gate h) and the delayed opening of K⁺ channels (gate n), which repolarise the membrane. The refractory period following each spike prevents backward propagation and limits firing rates.

2. Synaptic Transmission

When an action potential reaches a presynaptic terminal, voltage-gated Ca²⁺ channels open and calcium influx triggers the fusion of neurotransmitter-filled vesicles with the membrane (exocytosis). Neurotransmitter diffuses across the 20 nm synaptic cleft and binds to postsynaptic receptors, generating an excitatory or inhibitory postsynaptic potential (EPSP or IPSP).

Quantal model (del Castillo & Katz 1954):
  EPSP amplitude A = n · p · q
  n = number of release sites
  p = release probability per site
  q = quantal amplitude (single vesicle response)

Coefficient of variation: CV = √[(1−p)/(np)]
  (high CV → low p, few sites; diagnostic of presynaptic locus)

Short-term synaptic depression (STD):
  dx/dt = (1 − x)/τ_D − u_SE · x · δ(t − t_sp)
  u_SE = base utilisation fraction (typically 0.5)
  x decreases with each spike: fatigue at repetitive stimulation

Short-term facilitation (STF):
  du/dt = −u/τ_F + U(1−u) δ(t − t_sp)
  u increases with each spike: Ca²+ accumulation in terminal

NMDA receptor: dual Mg²+/ligand gate
  Activates only when Vm > −50 mV AND glutamate bound
  Ca²+ permeable: coincidence detector for plasticity
          

3. Hebbian Learning and Long-Term Potentiation

Donald Hebb proposed in 1949 that synapses strengthen when the pre- and postsynaptic neurons fire together: “neurons that fire together, wire together.” This principle was given a cellular basis when Bliss and Lømo discovered long-term potentiation (LTP) in the rabbit hippocampus in 1973.

Hebb rule:
  Δw_ij = η x_i x_j   (x = firing rate)
  Problem: unbounded weight growth (no forgetting)

BCM (Bienenstock-Cooper-Munro 1982) theory:
  Δw_ij = φ(x_j, θ_M) x_i
  φ < 0 if x_j < θ_M (depression);  φ > 0 if x_j > θ_M (potentiation)
  θ_M is a sliding threshold: θ_M ~ ⟨x_j²⟩
  Explains ocular dominance plasticity, orientation tuning

Spike-Timing-Dependent Plasticity (STDP, Bi & Poo 1998):
  Δw = A+ exp(−Δt/τ+)   if Δt = t_post − t_pre > 0  (LTP)
  Δw = −A− exp(Δt/τ−)   if Δt < 0  (LTD)
  A+ ≈ 0.005,  τ+ ≈ 20 ms,  A− ≈ 0.005,  τ− ≈ 20 ms

NMDA receptor as coincidence detector:
  Requires simultaneous glutamate AND postsynaptic depolarisation
  → natural implementation of Hebb’s temporal correlation rule
          

LTP is considered the cellular basis of learning and memory, particularly explicit memory mediated by the hippocampus. NMDA receptor-dependent LTP at CA1 synapses requires postsynaptic Ca²⁺ influx, which activates CaMKII — a kinase that phosphorylates AMPA receptors and recruits additional AMPA receptors to the synapse, thus increasing synaptic strength. Long-term depression (LTD), induced by low-frequency stimulation, reverses this process.

4. Neural Oscillations

Brain rhythms span multiple frequency bands: delta (0.5–4 Hz, deep sleep), theta (4–8 Hz, hippocampal navigation and working memory), alpha (8–12 Hz, cortical idling), beta (13–30 Hz, sensorimotor), and gamma (30–100 Hz, feature binding and attention). These rhythms emerge from the synchronised activity of neuronal populations and shape information processing.

Wilson-Cowan model (1972) — mean-field E-I network:
  τ_E dE/dt = −E + (1 − r_E E) S_E(w_EE E − w_EI I + I_E)
  τ_I dI/dt = −I + (1 − r_I I) S_I(w_IE E − w_II I + I_I)
  S(x) = 1/(1 + exp(−(x−θ)/σ))  (sigmoidal gain function)

Kuramoto model of N coupled oscillators:
  dθ_i/dt = ω_i + (K/N) ∑_j sin(θ_j − θ_i)
  Order parameter: r = (1/N)|∑_j exp(iθ_j)|
  Phase transition at K_c = 2/(πg(Ω)) where g = freq distribution width

Cross-frequency coupling:
  Theta/gamma coupling: gamma amplitude modulated by theta phase
  Modulation index MI = |⟨A_γ e^{iφ_θ}⟩|
  Hippocampal theta (δ6 Hz) organises gamma (δ80 Hz) sequences
          

5. Connectomics

Connectomics is the systematic mapping of all synaptic connections in a nervous system. Only two complete connectomes exist: C. elegans (302 neurons, ~7000 synapses, White et al. 1986) and the adult Drosophila brain (139,255 neurons, ~54.5 million synapses, FlyWire consortium 2023). Human whole-brain connectomics at synaptic resolution remains a grand challenge.

Structural connectivity matrix W_ij:
  W_ij = number of fibres (DTI tractography) or synapse count
  Adjacency A_ij = 1 if W_ij > threshold, else 0

Graph measures:
  Degree k_i = ∑_j A_ij
  Path length L = ⟨d_ij⟩  (d = shortest path)
  Clustering coefficient C_i = (triangles through i) / (k_i(k_i−1)/2)
  Small-world index σ = (C/C_rand) / (L/L_rand)

C. elegans connectome properties:
  Density 0.072  (vs random: same density)
  Clustering C = 0.308  (vs random C_rand = 0.054)
  Path length L = 2.65  (vs random L_rand = 2.25)
  → Small-world network (σ ≈ 5.6)

Rich club organisation:
  Hubs (high-degree nodes) are more densely interconnected
  than expected: φ(k) = E_>k / [k_>(k_>−1)/2]
  Human connectome rich club: precuneus, superior frontal, cingulate
          

The Human Connectome Project (HCP, 2010–2017) produced high-resolution diffusion MRI tractography and resting-state fMRI for 1200 subjects, enabling statistical analysis of structural-functional relationships. The Allen Brain Atlas provides spatially resolved gene expression data that can be integrated with connectivity to understand how molecular profiles relate to network topology.

6. Whole-Brain Modelling

The Virtual Brain (TVB) platform (Sanz Leon et al. 2013) simulates resting-state brain dynamics by placing neural mass models (Wilson-Cowan, Jansen-Rit, Hindmarsh-Rose) at each parcel of a parcellated connectome. The structural connectivity matrix sets the coupling strengths between regions, and simulated BOLD signals can be compared directly to empirical fMRI.

Neural mass model (per region i):
  Each region described by mean excitatory E_i and inhibitory I_i rate
  Coupling between regions: input to i from j weighted by W_ij

Combined dynamics:
  τ dE_i/dt = −E_i + S(∑_j w_ij E_j + w_EE E_i − w_EI I_i + I_ext)
  τ dI_i/dt = −I_i + S(w_IE E_i − w_II I_i)

Haemodynamic model (Balloon-Windkessel, Buxton-Friston):
  df/dt = s − (f−1)/τ_s − (f−f_0)/τ_f
  dv/dt = (f − v^(1/α)) / τ_0
  dq/dt = (f·E(f)/E_0 − q·v^(1/α−1)) / τ_0
  BOLD = V_0 [k_1(1−q) + k_2(1−q/v) + k_3(1−v)]

Model fitting: optimise global coupling G, local E/I balance
  to reproduce empirical functional connectivity matrix FC = corr(BOLD_i, BOLD_j)
          

Try These Simulations