Part IV: Electrophysiology & Signaling | Chapter 1

Electrophysiology & Patch Clamp

Patch clamp configurations, single-channel conductance, Markov gating models, and noise analysis — the quantitative foundations of ion channel biophysics

Why Electrophysiology?

The patch clamp technique, developed by Neher and Sakmann (Nobel Prize 1991), revolutionised neuroscience by enabling direct measurement of currents through single ion channels. With picoampere resolution and sub-millisecond time resolution, patch clamping reveals the stochastic behaviour of individual protein molecules in their native membrane environment.

This chapter derives the electrical circuit models underlying patch clamp configurations, extracts single-channel properties from macroscopic recordings through noise analysis, and develops Markov models to describe channel gating kinetics.

1. Patch Clamp Configurations

The patch clamp technique uses a fire-polished glass micropipette pressed against a cell membrane to form an extremely tight seal. The seal resistance, typically in the gigaohm range (G$\Omega$ seal), electrically isolates the patch of membrane beneath the pipette tip, enabling measurement of currents as small as a few picoamperes.

The Four Configurations

Starting from the initial gigaohm seal (cell-attached configuration), three additional configurations can be achieved:

• Cell-attached: Pipette sealed to intact cell. Records single channels in the patch while preserving intracellular milieu. The patch potential is $V_{\text{patch}} = V_{\text{rest}} - V_{\text{pip}}$.
• Whole-cell: Membrane under the pipette ruptured by suction or voltage pulse. The pipette solution dialyses the cell interior. Records total current from all channels in the cell.
• Inside-out: After cell-attached, retract pipette to excise patch. The cytoplasmic face is now exposed to the bath, allowing manipulation of intracellular messengers (Ca$^{2+}$, ATP, etc.).
• Outside-out: After whole-cell, retract pipette; membrane reseals with extracellular face outward. Ideal for studying ligand-gated channels by rapid agonist application.

Derivation: Equivalent Circuit and Seal Resistance

The cell-attached configuration can be modelled as a circuit with three key resistances: the seal resistance $R_{\text{seal}}$, the patch resistance $R_{\text{patch}}$(dominated by channel conductance), and the rest-of-cell resistance $R_{\text{cell}}$.

The measured current consists of the channel current plus a leak through the seal:

$$I_{\text{measured}} = I_{\text{channel}} + I_{\text{leak}} = I_{\text{channel}} + \frac{V_{\text{pip}}}{R_{\text{seal}}}$$

For the channel current to dominate, we require $R_{\text{seal}} \gg R_{\text{channel}}$. A typical single channel has conductance $\gamma \sim 10$–$100$ pS, giving$R_{\text{channel}} \sim 10$–$100$ G$\Omega$ when open. Thus we need:

$$\boxed{R_{\text{seal}} \geq 10 \text{ G}\Omega \quad (\text{gigaohm seal})}$$

The signal-to-noise ratio for detecting single-channel events scales as:

$$\text{SNR} = \frac{i_{\text{channel}}}{\sqrt{4 k_B T \Delta f / R_{\text{seal}} + S_{\text{amp}} \Delta f}}$$

where $i_{\text{channel}}$ is the single-channel current, $\Delta f$ is the recording bandwidth, and $S_{\text{amp}}$ is the amplifier noise spectral density. The thermal (Johnson) noise from the seal resistance sets the fundamental noise floor. For $R_{\text{seal}} = 10$ G$\Omega$ at 1 kHz bandwidth:$\sigma_I = \sqrt{4 k_B T \Delta f / R_{\text{seal}}} \approx 0.13$ pA, well below typical single-channel currents of 1–10 pA.

Derivation: Access Resistance and Series Resistance Compensation

In whole-cell mode, the pipette tip provides access to the cell interior. The access resistance $R_a$ depends on the pipette geometry. For a pipette with tip opening radius $a$ and cone half-angle $\theta$, filled with solution of resistivity $\rho$:

$$R_a = \frac{\rho}{\pi a \tan\theta}$$

For a typical pipette ($a \approx 0.5$ $\mu$m, $\theta \approx 10°$,$\rho \approx 150$ $\Omega$$\cdot$cm):$R_a \approx 5$–$20$ M$\Omega$.

The access resistance causes a voltage error. The actual membrane potential differs from the command voltage:

$$V_m = V_{\text{cmd}} - I_m R_a$$

For a cell with whole-cell current $I_m = 10$ nA and $R_a = 10$ M$\Omega$, the voltage error is $\Delta V = 100$ mV — catastrophically large! Series resistance compensation applies a positive feedback circuit:

$$V_{\text{cmd,corrected}} = V_{\text{cmd}} + \alpha \cdot I_m \cdot R_a$$

where $\alpha$ is the compensation fraction (typically 70–90%). The effective residual series resistance is:

$$\boxed{R_{a,\text{eff}} = R_a (1 - \alpha)}$$

With 80% compensation, a 10 M$\Omega$ access resistance is reduced to an effective 2 M$\Omega$, bringing the voltage error for 10 nA to an acceptable 20 mV. However, excessive compensation ($\alpha > 0.9$) risks oscillation due to the positive feedback loop. The maximum stable compensation depends on the product$R_a C_m$ (the time constant of the series RC filter formed by access resistance and membrane capacitance).

Whole-Cell Capacitance Transients

When a voltage step $\Delta V$ is applied in whole-cell mode, the current response has a fast capacitive transient followed by the steady-state ionic current. The capacitive current decays exponentially:

$$I_{\text{cap}}(t) = \frac{\Delta V}{R_a} \exp\left(-\frac{t}{R_a C_m}\right)$$

From the transient, we extract the cell's electrical parameters:

• Access resistance: $R_a = \Delta V / I_{\text{peak}}$ from the peak current
• Membrane capacitance: $C_m = \tau / R_a$ where $\tau$ is the decay time constant
• Membrane resistance: $R_m = \Delta V / I_{\text{ss}} - R_a$ from the steady-state current

Typical values for a small neuron: $C_m \approx 20$ pF (surface area$\approx 2000$ $\mu$m$^2$ at 1 $\mu$F/cm$^2$),$R_m \approx 500$ M$\Omega$, giving$\tau_m = R_m C_m \approx 10$ ms.

2. Single-Channel Current-Voltage Relationships

The current through a single open channel depends on the membrane potential, ion concentrations, and channel properties. Analysis of the single-channel I-V curve reveals the conductance, selectivity, and rectification properties.

Derivation: Conductance and Reversal Potential

For a channel that obeys Ohm's law, the single-channel current is:

$$i = \gamma (V_m - E_{\text{rev}})$$

where $\gamma$ is the single-channel conductance (slope of the I-V curve) and$E_{\text{rev}}$ is the reversal potential (voltage at which current reverses direction). The conductance is extracted from the slope:

$$\boxed{\gamma = \frac{di}{dV_m} = \frac{i(V_2) - i(V_1)}{V_2 - V_1}}$$

The reversal potential identifies the ion(s) carried by the channel. For a perfectly selective channel, $E_{\text{rev}}$ equals the Nernst potential of the permeant ion. For a channel permeable to multiple ions, the reversal potential is given by the Goldman-Hodgkin-Katz voltage equation. Experimentally, shifting ion concentrations and observing the change in $E_{\text{rev}}$ confirms the permeant ion.

For a K$^+$-selective channel:

$$E_{\text{rev}} = E_K = \frac{RT}{F}\ln\frac{[\text{K}^+]_o}{[\text{K}^+]_i}$$

At room temperature with $[\text{K}^+]_o = 5$ mM and $[\text{K}^+]_i = 140$ mM,$E_K \approx -86$ mV.

Derivation: Rectification from Asymmetric Solutions

Many channels do not show simple ohmic (linear) I-V relationships. Rectification arises from two sources: (1) asymmetric ion concentrations, and (2) intrinsic channel properties.

GHK current equation rectification: Even for a simple pore with no voltage-dependent gating, asymmetric ion concentrations produce a nonlinear I-V. The Goldman-Hodgkin-Katz current equation gives:

$$i = P_s z^2 \frac{F^2 V_m}{RT} \frac{c_i - c_o \exp(-zFV_m/RT)}{1 - \exp(-zFV_m/RT)}$$

where $P_s$ is the single-channel permeability. This produces rectification: with symmetric solutions ($c_i = c_o$), the I-V is linear, but with asymmetric solutions, outward current is larger when $c_i > c_o$ (outward rectification from the perspective of ion flux).

Intrinsic (block) rectification: Inward rectifier K$^+$ channels (K$_{\text{ir}}$) are blocked by intracellular Mg$^{2+}$ and polyamines at depolarised potentials. The effective conductance is:

$$\gamma_{\text{eff}}(V) = \gamma_0 \cdot \frac{1}{1 + [\text{Mg}^{2+}]_i / K_d(V)}$$

where the voltage-dependent dissociation constant is:

$$K_d(V) = K_d(0) \exp\left(\frac{z_B \delta F V}{RT}\right)$$

Here $z_B$ is the charge of the blocker, $\delta$ is the fractional electrical distance of the blocking site from the intracellular side (typically$\delta \approx 0.4$–$0.6$). This produces strong inward rectification: channels conduct inward current readily but pass very little outward current.

3. Channel Gating Kinetics — Markov Models

Ion channels fluctuate stochastically between conducting (open) and non-conducting (closed, inactivated) states. Markov models describe these transitions as memoryless processes with voltage- and ligand-dependent rate constants.

Derivation: Two-State Model (C ⇌ O)

The simplest gating scheme has a single closed (C) and open (O) state:

$$\text{C} \underset{\beta}{\overset{\alpha}{\rightleftharpoons}} \text{O}$$

where $\alpha$ is the opening rate and $\beta$ is the closing rate. The master equation for the probability of being open is:

$$\frac{dP_O}{dt} = \alpha(1 - P_O) - \beta P_O$$

At steady state ($dP_O/dt = 0$):

$$\boxed{P_O^{\infty} = \frac{\alpha}{\alpha + \beta}}$$

The time course following a step change in voltage (which instantly changes $\alpha$and $\beta$) is:

$$P_O(t) = P_O^{\infty} + (P_O(0) - P_O^{\infty})\exp(-t/\tau)$$

where the relaxation time constant is:

$$\boxed{\tau = \frac{1}{\alpha + \beta}}$$

Dwell-time distributions: In a single-channel recording, the channel alternates between open and closed dwell times. For a two-state model, both distributions are single exponentials:

$$f_O(t) = \beta \exp(-\beta t), \quad \langle t_O \rangle = 1/\beta$$

$$f_C(t) = \alpha \exp(-\alpha t), \quad \langle t_C \rangle = 1/\alpha$$

This follows from the memoryless property of the Markov process: the probability of remaining in a state for time $t$ is $\exp(-kt)$ where $k$is the total exit rate from that state. By fitting exponentials to measured dwell-time histograms, we extract the rate constants.

Derivation: Three-State Model (C ⇌ O ⇌ I)

Many voltage-gated channels have an inactivated state that the channel enters from the open state but from which it cannot directly conduct:

$$\text{C} \underset{\beta}{\overset{\alpha}{\rightleftharpoons}} \text{O} \underset{\delta}{\overset{\gamma_I}{\rightleftharpoons}} \text{I}$$

The master equations become a system of ODEs:

$$\frac{dP_C}{dt} = -\alpha P_C + \beta P_O$$

$$\frac{dP_O}{dt} = \alpha P_C - (\beta + \gamma_I) P_O + \delta P_I$$

$$\frac{dP_I}{dt} = \gamma_I P_O - \delta P_I$$

with constraint $P_C + P_O + P_I = 1$. In matrix form:

$$\frac{d\mathbf{P}}{dt} = \mathbf{Q} \cdot \mathbf{P}, \quad \mathbf{Q} = \begin{pmatrix} -\alpha & \beta & 0 \\ \alpha & -(\beta + \gamma_I) & \delta \\ 0 & \gamma_I & -\delta \end{pmatrix}$$

The eigenvalues of $\mathbf{Q}$ give the relaxation time constants. One eigenvalue is always zero (conservation of probability), and the two nonzero eigenvalues$\lambda_1, \lambda_2$ give the fast and slow components of the relaxation:

$$P_O(t) = P_O^{\infty} + A_1 e^{\lambda_1 t} + A_2 e^{\lambda_2 t}$$

Dwell-time consequences: The open-time distribution remains a single exponential (since there is only one open state) with rate constant$\beta + \gamma_I$:

$$\boxed{f_O(t) = (\beta + \gamma_I) \exp(-(\beta + \gamma_I) t)}$$

However, the closed-time distribution is now biexponential because the channel can be in either C or I, each with different exit rates. This is a key diagnostic: if closed-time histograms require multiple exponentials, there are multiple closed/inactivated states.

Maximum Likelihood Estimation of Rate Constants

Given a sequence of $N$ open times $\{t_1^O, t_2^O, \ldots\}$ and closed times $\{t_1^C, t_2^C, \ldots\}$, the log-likelihood for the two-state model is:

$$\ln \mathcal{L}(\alpha, \beta) = \sum_j \ln[\beta \exp(-\beta t_j^O)] + \sum_k \ln[\alpha \exp(-\alpha t_k^C)]$$

$$= N_O \ln \beta - \beta \sum_j t_j^O + N_C \ln \alpha - \alpha \sum_k t_k^C$$

Maximising with respect to $\beta$ and $\alpha$:

$$\frac{\partial \ln \mathcal{L}}{\partial \beta} = \frac{N_O}{\beta} - \sum_j t_j^O = 0 \implies \boxed{\hat{\beta} = \frac{N_O}{\sum_j t_j^O} = \frac{1}{\langle t_O \rangle}}$$

$$\hat{\alpha} = \frac{N_C}{\sum_k t_k^C} = \frac{1}{\langle t_C \rangle}$$

The MLE estimates are simply the reciprocals of the mean dwell times — an intuitive result for exponential distributions. For multi-state models, the likelihood involves matrix exponentials and requires numerical optimisation, but the principle remains the same.

4. Noise Analysis

When a membrane patch contains many channels, individual openings and closings cannot be resolved. However, the stochastic gating produces fluctuations in the macroscopic current. Noise analysis extracts single-channel properties from these fluctuations without needing to resolve individual events.

Derivation: Single-Channel Current from Macroscopic Variance

Consider $N$ identical, independent two-state channels, each with open probability$P_O$ and single-channel current $i$. The number of open channels follows a binomial distribution. The macroscopic current is:

$$I = N i P_O$$

The variance of the macroscopic current is:

$$\sigma^2 = N i^2 P_O (1 - P_O)$$

Dividing variance by mean to eliminate $N$:

$$\frac{\sigma^2}{I} = i(1 - P_O)$$

Rearranging:

$$\boxed{i = \frac{\sigma^2}{I(1 - P_O)}}$$

If $P_O$ is small (as for many channels at low agonist concentrations), then$1 - P_O \approx 1$ and:

$$i \approx \frac{\sigma^2}{I} \quad (\text{low } P_O \text{ approximation})$$

We can also find $N$ by combining mean and variance:

$$N = \frac{I^2}{\sigma^2 \cdot \frac{P_O}{1-P_O} + I \cdot i} = \frac{I}{i P_O}$$

The mean-variance parabola provides a powerful method. By varying $P_O$ (e.g., by changing agonist concentration or voltage) and plotting $\sigma^2$ vs $I$:$$\sigma^2 = iI - \frac{I^2}{N}$$This is a parabola that opens downward. The initial slope gives $i$ and the x-intercept gives $Ni$ (the maximum current when all channels are open).

Campbell's Theorem

Campbell's theorem provides a frequency-domain approach to noise analysis. For a random process consisting of independent, identical events (channel openings) occurring at rate$\nu$, each producing a current pulse $h(t)$:

$$\text{Mean: } \langle I \rangle = \nu \int_0^{\infty} h(t)\, dt$$

$$\text{Variance: } \sigma^2 = \nu \int_0^{\infty} h^2(t)\, dt$$

For rectangular pulses of amplitude $i$ and mean duration $\tau_O = 1/\beta$:

$$\langle I \rangle = \nu i \tau_O, \quad \sigma^2 = \nu i^2 \tau_O$$

The power spectral density of the current fluctuations is:

$$S(f) = \nu |H(f)|^2$$

where $H(f)$ is the Fourier transform of $h(t)$. For a two-state channel, $S(f)$ is a Lorentzian with corner frequency $f_c = (\alpha + \beta)/(2\pi)$:$$S(f) = \frac{S_0}{1 + (f/f_c)^2}$$Fitting the power spectrum yields both the relaxation rate ($\alpha + \beta$) and, combined with the variance, the single-channel current.

5. Applications

Drug Screening & Pharmacology

Patch clamp is the gold standard for characterising drug effects on ion channels. Automated patch clamp platforms (e.g., Nanion, Sophion) enable high-throughput screening of thousands of compounds per day. Key measurements include:

• IC$_{50}$ determination: Concentration-response curves for channel block, fitted to the Hill equation $f_{\text{block}} = 1/(1 + (\text{IC}_{50}/[\text{drug}])^{n_H})$
• Use-dependent block: Drugs that preferentially block open or inactivated states show frequency-dependent effects
• hERG safety screening: All new drugs must be tested for block of the hERG K$^+$ channel, which can cause lethal cardiac arrhythmias (Long QT syndrome)

Channelopathies

Mutations in ion channel genes cause a wide range of diseases (channelopathies). Patch clamp characterisation of mutant channels reveals the biophysical mechanism:

• Epilepsy: Gain-of-function Na$^+$ channel mutations (SCN1A) or loss-of-function K$^+$ channel mutations increase neuronal excitability
• Cystic fibrosis: CFTR (a Cl$^-$ channel) mutations reduce channel expression or conductance
• Long QT syndrome: Mutations in cardiac K$^+$ channels (KCNQ1, hERG) delay repolarisation
• Myotonia: Cl$^-$ channel (ClC-1) mutations reduce muscle membrane conductance, causing hyperexcitability

Cardiac Electrophysiology

The cardiac action potential involves orchestrated activity of multiple ion channels. Patch clamp studies of individual cardiac channel types have revealed the molecular basis of each phase:

• Phase 0 (upstroke): Na$^+$ channel (Na$_v$1.5) opening, rising at $\sim 300$ V/s
• Phase 1 (early repolarisation): Transient outward K$^+$ current (I$_{\text{to}}$)
• Phase 2 (plateau): Balance of L-type Ca$^{2+}$ current and delayed rectifier K$^+$ current
• Phase 3 (repolarisation): I$_{Kr}$ (hERG) and I$_{Ks}$ (KCNQ1) K$^+$ currents dominate
• Phase 4 (resting): I$_{K1}$ (inward rectifier) maintains resting potential at $\sim -85$ mV

6. Python Simulation: Single-Channel Recording Analysis

The following simulation generates a stochastic single-channel recording (random telegraph noise), computes open and closed dwell-time histograms, fits exponential distributions, and performs noise analysis on macroscopic currents from multiple channels.

Single-Channel Recording Simulation & Analysis

Python

script.py255 lines

import numpy as np

# =========================================================
# Single-Channel Recording Simulation & Analysis
# =========================================================

print("=" * 60)
print("SINGLE-CHANNEL RECORDING SIMULATION")
print("=" * 60)

np.random.seed(42)

# Channel parameters
alpha = 500.0    # opening rate (s^-1), C -> O
beta = 200.0     # closing rate (s^-1), O -> C
gamma_sc = 25.0  # single-channel conductance (pS)
V_m = -60.0      # membrane potential (mV)
E_rev = -85.0    # reversal potential (mV, K+ channel)

# Single-channel current amplitude
i_sc = gamma_sc * (V_m - E_rev) * 1e-3  # convert pS*mV to pA
print(f"\nChannel parameters:")
print(f"  Opening rate alpha = {alpha} s^-1")
print(f"  Closing rate beta  = {beta} s^-1")
print(f"  Conductance gamma  = {gamma_sc} pS")
print(f"  V_m = {V_m} mV, E_rev = {E_rev} mV")
print(f"  Single-channel current i = gamma*(V-E_rev) = {i_sc:.3f} pA")

# Theoretical values
Po_theory = alpha / (alpha + beta)
tau_theory = 1.0 / (alpha + beta)
mean_open = 1.0 / beta
mean_closed = 1.0 / alpha
print(f"\nTheoretical predictions:")
print(f"  P_open = alpha/(alpha+beta) = {Po_theory:.4f}")
print(f"  Relaxation tau = 1/(alpha+beta) = {tau_theory*1000:.3f} ms")
print(f"  Mean open time = 1/beta = {mean_open*1000:.3f} ms")
print(f"  Mean closed time = 1/alpha = {mean_closed*1000:.3f} ms")

# --- Simulate single-channel recording ---
print("\n" + "=" * 60)
print("GENERATING SINGLE-CHANNEL RECORDING")
print("=" * 60)

dt = 0.00005  # 50 us time step (20 kHz sampling)
T_total = 2.0  # 2 seconds total recording
N = int(T_total / dt)

# Generate state trajectory using Gillespie-like approach
states = np.zeros(N, dtype=int)  # 0 = closed, 1 = open
state = 0  # start closed

for k in range(1, N):
    if state == 0:  # closed
        p_open = alpha * dt
        if np.random.rand() < p_open:
            state = 1
    else:  # open
        p_close = beta * dt
        if np.random.rand() < p_close:
            state = 0
    states[k] = state

# Generate current trace with noise
noise_std = 0.15  # pA (amplifier noise)
current = states * i_sc + np.random.normal(0, noise_std, N)

# Compute open probability
Po_observed = np.mean(states)
print(f"\nSimulated recording ({T_total} s at {1/dt/1000:.0f} kHz):")
print(f"  Observed P_open = {Po_observed:.4f}")
print(f"  Theory P_open   = {Po_theory:.4f}")

# --- Extract dwell times ---
print("\n" + "=" * 60)
print("DWELL TIME ANALYSIS")
print("=" * 60)

# Find transitions
transitions = np.where(np.diff(states) != 0)[0]
dwell_times = np.diff(transitions) * dt * 1000  # in ms

# Separate open and closed dwell times
open_dwells = []
closed_dwells = []
for k in range(len(transitions) - 1):
    duration = (transitions[k+1] - transitions[k]) * dt * 1000  # ms
    if states[transitions[k] + 1] == 1:  # open state
        open_dwells.append(duration)
    else:
        closed_dwells.append(duration)

open_dwells = np.array(open_dwells)
closed_dwells = np.array(closed_dwells)

print(f"\nNumber of openings detected: {len(open_dwells)}")
print(f"Number of closings detected: {len(closed_dwells)}")

# Fit exponentials (MLE = 1/mean for exponential)
if len(open_dwells) > 5:
    mean_open_obs = np.mean(open_dwells)
    beta_fit = 1.0 / mean_open_obs  # ms^-1
    print(f"\nOpen time distribution:")
    print(f"  Mean open time = {mean_open_obs:.4f} ms")
    print(f"  Fitted beta = 1/<t_open> = {beta_fit:.1f} ms^-1 = {beta_fit*1000:.1f} s^-1")
    print(f"  True beta = {beta} s^-1")
    print(f"  Error: {abs(beta_fit*1000 - beta)/beta * 100:.1f}%")

if len(closed_dwells) > 5:
    mean_closed_obs = np.mean(closed_dwells)
    alpha_fit = 1.0 / mean_closed_obs  # ms^-1
    print(f"\nClosed time distribution:")
    print(f"  Mean closed time = {mean_closed_obs:.4f} ms")
    print(f"  Fitted alpha = 1/<t_closed> = {alpha_fit:.1f} ms^-1 = {alpha_fit*1000:.1f} s^-1")
    print(f"  True alpha = {alpha} s^-1")
    print(f"  Error: {abs(alpha_fit*1000 - alpha)/alpha * 100:.1f}%")

# Dwell time histogram analysis
print("\nOpen dwell time histogram (ms bins):")
if len(open_dwells) > 10:
    max_t = min(np.percentile(open_dwells, 99), 30)
    bins = np.linspace(0, max_t, 15)
    counts_o, edges_o = np.histogram(open_dwells, bins=bins)
    bin_centers_o = 0.5 * (edges_o[:-1] + edges_o[1:])
    print(f"  {'Bin center (ms)':>16} {'Count':>8} {'Expected':>10}")
    total_o = len(open_dwells)
    bin_w = bins[1] - bins[0]
    for bc, c in zip(bin_centers_o, counts_o):
        expected = total_o * (np.exp(-beta * bc / 1000 * 1) - np.exp(-beta * (bc + bin_w/2) / 1000 * 1))
        # simple expected from exponential
        exp_val = total_o * beta/1000 * np.exp(-beta/1000 * bc) * bin_w
        print(f"  {bc:>16.2f} {c:>8d} {exp_val:>10.1f}")

# --- Noise Analysis (multiple channels) ---
print("\n" + "=" * 60)
print("NOISE ANALYSIS: MACROSCOPIC CURRENT")
print("=" * 60)

N_channels = 100
print(f"\nSimulating {N_channels} independent channels...")

# Shorter simulation for speed
dt2 = 0.0001  # 100 us
T2 = 0.5  # 500 ms
N2 = int(T2 / dt2)

# Vary P_o by changing voltage (which changes alpha)
# Simulate at several P_o values
Po_values = [0.1, 0.2, 0.3, 0.5, 0.7, 0.9]
mean_currents = []
var_currents = []

print(f"\n{'P_o':>6} {'I_mean (pA)':>12} {'Var (pA^2)':>12} {'i_est (pA)':>12} {'N_est':>8}")
print("-" * 55)

for Po_target in Po_values:
    # Adjust alpha to get target Po
    alpha_adj = Po_target * beta / (1 - Po_target)

# Simulate N_channels
    macro_current = np.zeros(N2)
    for ch in range(N_channels):
        ch_state = 0
        for t_idx in range(N2):
            if ch_state == 0:
                if np.random.rand() < alpha_adj * dt2:
                    ch_state = 1
            else:
                if np.random.rand() < beta * dt2:
                    ch_state = 0
            macro_current[t_idx] += ch_state * i_sc

I_mean = np.mean(macro_current)
    I_var = np.var(macro_current)
    mean_currents.append(I_mean)
    var_currents.append(I_var)

# Estimate i from variance/mean
    if I_mean > 0.01:
        i_est = I_var / (I_mean * (1 - I_mean / (N_channels * i_sc)))
        N_est = I_mean / (i_est * Po_target) if i_est > 0 else 0
    else:
        i_est = 0
        N_est = 0

print(f"{Po_target:>6.1f} {I_mean:>12.3f} {I_var:>12.4f} {i_est:>12.4f} {N_est:>8.0f}")

# Fit mean-variance parabola: var = i*I - I^2/N
print("\nMean-Variance Parabola Fit:")
print("  sigma^2 = i*I - I^2/N")

I_arr = np.array(mean_currents)
V_arr = np.array(var_currents)

# Fit quadratic: var = a*I + b*I^2
# Using least squares: [I, I^2] * [a; b] = var
A_mat = np.column_stack([I_arr, I_arr**2])
coeffs = np.linalg.lstsq(A_mat, V_arr, rcond=None)[0]
i_fit = coeffs[0]
N_fit = -1.0 / coeffs[1]

print(f"  Fitted i = {i_fit:.4f} pA (true: {i_sc:.3f} pA)")
print(f"  Fitted N = {N_fit:.0f} channels (true: {N_channels})")
print(f"  Error in i: {abs(i_fit - i_sc) / i_sc * 100:.1f}%")
print(f"  Error in N: {abs(N_fit - N_channels) / N_channels * 100:.1f}%")

# --- Power Spectrum Analysis ---
print("\n" + "=" * 60)
print("POWER SPECTRUM ANALYSIS")
print("=" * 60)

# Single channel current for spectral analysis
channel_trace = states[:20000].astype(float) * i_sc
fs = 1.0 / dt  # sampling frequency

# Compute power spectrum using numpy FFT
n_fft = len(channel_trace)
freqs = np.fft.rfftfreq(n_fft, d=dt)
fft_vals = np.fft.rfft(channel_trace - np.mean(channel_trace))
psd = 2.0 * np.abs(fft_vals)**2 / (n_fft * fs)  # one-sided PSD

# Bin the spectrum for cleaner display
n_bins_spec = 20
log_freq_bins = np.logspace(np.log10(freqs[1]), np.log10(freqs[-1]), n_bins_spec + 1)
binned_freq = []
binned_psd = []
for b in range(n_bins_spec):
    mask = (freqs >= log_freq_bins[b]) & (freqs < log_freq_bins[b+1])
    if np.sum(mask) > 0:
        binned_freq.append(np.mean(freqs[mask]))
        binned_psd.append(np.mean(psd[mask]))

# Theoretical Lorentzian
f_c_theory = (alpha + beta) / (2 * np.pi)
S0_theory = 4 * i_sc**2 * Po_theory * (1 - Po_theory) * tau_theory

print(f"\nCorner frequency f_c = (alpha+beta)/(2*pi) = {f_c_theory:.1f} Hz")
print(f"DC spectral density S(0) = {S0_theory:.6f} pA^2/Hz")

print(f"\n{'Freq (Hz)':>12} {'PSD (pA^2/Hz)':>16} {'Lorentzian':>14}")
print("-" * 45)
for f_val, p_val in zip(binned_freq[:12], binned_psd[:12]):
    lor = S0_theory / (1 + (f_val / f_c_theory)**2)
    print(f"{f_val:>12.1f} {p_val:>16.8f} {lor:>14.8f}")

print("\n" + "=" * 60)
print("SUMMARY")
print("=" * 60)
print(f"Single-channel current: {i_sc:.3f} pA")
print(f"Conductance: {gamma_sc} pS")
print(f"Open probability: {Po_theory:.3f} (theory), {Po_observed:.4f} (observed)")
print(f"Rate constants recovered from dwell times agree")
print(f"Noise analysis recovers i and N from macroscopic fluctuations")
print(f"Power spectrum is Lorentzian with f_c = {f_c_theory:.0f} Hz")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

7. Three-State Markov Model Simulation

This simulation implements the C ⇌ O ⇌ I three-state Markov model, generates dwell-time histograms, and demonstrates how inactivation produces biexponential closed-time distributions.

Three-State Markov Model: C ↔ O ↔ I

Python

script.py228 lines

import numpy as np

# =========================================================
# Three-State Markov Model: C <-> O <-> I
# =========================================================

print("=" * 60)
print("THREE-STATE MARKOV MODEL SIMULATION")
print("=" * 60)

np.random.seed(123)

# Rate constants (s^-1)
alpha = 800.0    # C -> O
beta_r = 300.0   # O -> C
gamma_i = 150.0  # O -> I (inactivation)
delta = 50.0     # I -> O (recovery)

print(f"\nRate constants:")
print(f"  C -> O: alpha    = {alpha} s^-1")
print(f"  O -> C: beta     = {beta_r} s^-1")
print(f"  O -> I: gamma_I  = {gamma_i} s^-1")
print(f"  I -> O: delta    = {delta} s^-1")

# Build Q matrix
Q = np.array([
    [-alpha, beta_r, 0],
    [alpha, -(beta_r + gamma_i), delta],
    [0, gamma_i, -delta]
])

print(f"\nQ matrix:")
for row in Q:
    print(f"  [{row[0]:>8.1f} {row[1]:>8.1f} {row[2]:>8.1f}]")

# Steady-state: solve Q^T * pi = 0, sum(pi) = 1
# Use null space approach
# Add constraint row
A_eq = np.vstack([Q.T, np.ones(3)])
b_eq = np.zeros(4)
b_eq[3] = 1.0
# Solve overdetermined system
pi_ss = np.linalg.lstsq(A_eq, b_eq, rcond=None)[0]

print(f"\nSteady-state probabilities:")
print(f"  P_C = {pi_ss[0]:.4f}")
print(f"  P_O = {pi_ss[1]:.4f}")
print(f"  P_I = {pi_ss[2]:.4f}")
print(f"  Sum = {np.sum(pi_ss):.4f}")

# Eigenvalues for relaxation times
eigenvalues = np.linalg.eigvals(Q)
eigenvalues_sorted = np.sort(np.real(eigenvalues))
print(f"\nEigenvalues of Q:")
for ev in eigenvalues_sorted:
    if abs(ev) > 1e-10:
        print(f"  lambda = {ev:.2f} s^-1, tau = {-1/ev*1000:.3f} ms")
    else:
        print(f"  lambda = {ev:.6f} (steady state)")

# --- Simulate single-channel trajectory ---
print("\n" + "=" * 60)
print("SIMULATING SINGLE-CHANNEL TRAJECTORY")
print("=" * 60)

dt = 0.00005  # 50 us
T_total = 5.0  # 5 seconds
N = int(T_total / dt)

state = 0  # start in C
states = np.zeros(N, dtype=int)

for k in range(N):
    states[k] = state
    r = np.random.rand()
    if state == 0:  # Closed
        if r < alpha * dt:
            state = 1  # -> Open
    elif state == 1:  # Open
        if r < beta_r * dt:
            state = 0  # -> Closed
        elif r < (beta_r + gamma_i) * dt:
            state = 2  # -> Inactivated
    else:  # Inactivated
        if r < delta * dt:
            state = 1  # -> Open

# Observed occupancies
P_C_obs = np.mean(states == 0)
P_O_obs = np.mean(states == 1)
P_I_obs = np.mean(states == 2)

print(f"\nObserved vs theoretical occupancies:")
print(f"  {'State':>8} {'Observed':>10} {'Theory':>10}")
print(f"  {'C':>8} {P_C_obs:>10.4f} {pi_ss[0]:>10.4f}")
print(f"  {'O':>8} {P_O_obs:>10.4f} {pi_ss[1]:>10.4f}")
print(f"  {'I':>8} {P_I_obs:>10.4f} {pi_ss[2]:>10.4f}")

# --- Extract dwell times by state ---
transitions = []
current_state = states[0]
start_idx = 0

for k in range(1, N):
    if states[k] != current_state:
        transitions.append((current_state, start_idx, k))
        current_state = states[k]
        start_idx = k

open_dwells = []
closed_dwells = []
inact_dwells = []

for st, s_idx, e_idx in transitions:
    dur = (e_idx - s_idx) * dt * 1000  # ms
    if st == 0:
        closed_dwells.append(dur)
    elif st == 1:
        open_dwells.append(dur)
    else:
        inact_dwells.append(dur)

open_dwells = np.array(open_dwells)
closed_dwells = np.array(closed_dwells)
inact_dwells = np.array(inact_dwells)

print(f"\nDwell time statistics:")
print(f"  Open dwells: n={len(open_dwells)}")
if len(open_dwells) > 5:
    mean_to = np.mean(open_dwells)
    exit_rate_O = beta_r + gamma_i
    print(f"    Mean = {mean_to:.4f} ms")
    print(f"    Expected = 1/(beta+gamma_I) = {1.0/exit_rate_O*1000:.4f} ms")
    print(f"    Fitted exit rate = {1.0/mean_to*1000:.1f} s^-1 (true: {exit_rate_O:.1f} s^-1)")

print(f"  Closed dwells: n={len(closed_dwells)}")
if len(closed_dwells) > 5:
    mean_tc = np.mean(closed_dwells)
    print(f"    Mean = {mean_tc:.4f} ms")
    print(f"    Expected = 1/alpha = {1.0/alpha*1000:.4f} ms")

print(f"  Inactivated dwells: n={len(inact_dwells)}")
if len(inact_dwells) > 5:
    mean_ti = np.mean(inact_dwells)
    print(f"    Mean = {mean_ti:.4f} ms")
    print(f"    Expected = 1/delta = {1.0/delta*1000:.4f} ms")

# Combined "non-conducting" dwell time distribution
# The channel appears closed when in either C or I state
# This produces a biexponential distribution
print("\n" + "=" * 60)
print("BIEXPONENTIAL CLOSED-TIME DISTRIBUTION")
print("=" * 60)

# Combine C and I intervals that are adjacent (no intervening O)
# Actually, from the open channel perspective, each closing event
# ends up in either C or I, then returns to O
# The apparent closed time = time from O->C or O->I until next opening

# Rebuild apparent closed times (gap between openings)
open_mask = (states == 1)
transitions_to_open = np.where(np.diff(open_mask.astype(int)) == 1)[0]
transitions_from_open = np.where(np.diff(open_mask.astype(int)) == -1)[0]

if len(transitions_from_open) > 0 and len(transitions_to_open) > 1:
    # Gap = time from channel closing to next opening
    gaps = []
    for k in range(len(transitions_from_open)):
        close_time = transitions_from_open[k]
        # Find next opening after this closing
        next_opens = transitions_to_open[transitions_to_open > close_time]
        if len(next_opens) > 0:
            gap = (next_opens[0] - close_time) * dt * 1000  # ms
            gaps.append(gap)

gaps = np.array(gaps)
    if len(gaps) > 20:
        print(f"\nApparent gap (non-conducting) times: n={len(gaps)}")
        print(f"  Mean gap = {np.mean(gaps):.4f} ms")

# Histogram
        max_gap = min(np.percentile(gaps, 98), 50)
        bins_g = np.linspace(0, max_gap, 12)
        counts_g, edges_g = np.histogram(gaps, bins=bins_g)
        bin_c_g = 0.5 * (edges_g[:-1] + edges_g[1:])

print(f"\n  Gap time histogram:")
        print(f"  {'Bin (ms)':>10} {'Count':>8}")
        for bc, c in zip(bin_c_g, counts_g):
            bar = '#' * min(int(c / max(counts_g) * 30), 30)
            print(f"  {bc:>10.2f} {c:>8d} {bar}")

print(f"\n  Note: This distribution is BIEXPONENTIAL because")
        print(f"  the channel can be in C (exit rate alpha={alpha} s^-1)")
        print(f"  or I (exit rate delta={delta} s^-1)")
        print(f"  Short gaps ~ C state (tau = {1/alpha*1000:.2f} ms)")
        print(f"  Long gaps ~ I state (tau = {1/delta*1000:.2f} ms)")

# --- Macroscopic current relaxation ---
print("\n" + "=" * 60)
print("MACROSCOPIC CURRENT RELAXATION")
print("=" * 60)

# Solve ODE: dP/dt = Q*P for step from all-closed
t_relax = np.linspace(0, 20, 200)  # ms
P0 = np.array([1.0, 0.0, 0.0])  # start all closed

# Matrix exponential via eigendecomposition
eigvals, eigvecs = np.linalg.eig(Q)
eigvals_ms = eigvals / 1000  # convert to ms^-1

P_t = np.zeros((len(t_relax), 3))
for k, t_val in enumerate(t_relax):
    exp_diag = np.diag(np.exp(eigvals_ms * t_val))
    P_t[k] = np.real(eigvecs @ exp_diag @ np.linalg.inv(eigvecs) @ P0)

print(f"\nMacroscopic P_O(t) after step activation:")
print(f"  {'t (ms)':>8} {'P_C':>8} {'P_O':>8} {'P_I':>8}")
for idx in [0, 5, 10, 20, 40, 60, 80, 100, 150, 199]:
    t_val = t_relax[idx]
    print(f"  {t_val:>8.2f} {P_t[idx,0]:>8.4f} {P_t[idx,1]:>8.4f} {P_t[idx,2]:>8.4f}")

print(f"\n  P_O peaks at t ~ {t_relax[np.argmax(P_t[:,1])]:.2f} ms")
print(f"  Peak P_O = {np.max(P_t[:,1]):.4f}")
print(f"  Steady-state P_O = {pi_ss[1]:.4f}")
print(f"\n  The transient overshoot of P_O demonstrates")
print(f"  macroscopic inactivation: channels open then inactivate.")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Chapter Summary

• Patch clamp configurations (cell-attached, whole-cell, inside-out, outside-out) each offer unique experimental advantages. The gigaohm seal is essential for single-channel resolution.
• Access resistance causes voltage errors in whole-cell mode; series resistance compensation with $R_{a,\text{eff}} = R_a(1-\alpha)$ reduces but cannot eliminate this error.
• Single-channel I-V curves reveal conductance $\gamma$ and reversal potential $E_{\text{rev}}$; rectification arises from asymmetric solutions (GHK) or voltage-dependent block.
• Markov models describe gating kinetics; rate constants are extracted from dwell-time histograms via MLE ($\hat{\beta} = 1/\langle t_O \rangle$).
• Noise analysis extracts single-channel current from macroscopic variance: $i = \sigma^2 / I(1-P_O)$. The mean-variance parabola gives both $i$ and $N$.
• Campbell's theorem links mean and variance to the shot-noise properties of channel gating; the Lorentzian power spectrum has corner frequency $f_c = (\alpha+\beta)/(2\pi)$.

← X-ray Crystallography Neural Biophysics →