Hidden Markov Model

A Hidden Markov model (HMM) is a statistical generative model within computation and algorithms that represents sequential or time-series data as the observable output of an underlying, unobservable Markov chain. Formulated mathematically by Leonard E. Baum and colleagues at the Institute for Defense Analyses between 1966 and 1972 and popularized in speech recognition by Lawrence Rabiner (1989) and in computational molecular biology by Anders Krogh, David Haussler, Richard Durbin, and Sean Eddy (1994–1998), a Hidden Markov model couples a discrete latent state trajectory π = (π₁, π₂, ..., π_T) with state-dependent emission distributions that generate the observed sequence X = (x₁, x₂, ..., x_T).

In the natural sciences and bioinformatics, Hidden Markov models provide the probabilistic engine for remote protein domain homology detection (Pfam and HMMER profile HMMs), ab initio eukaryotic exon–intron gene finding (GENSCAN, Augustus), CpG island and chromatin-state genome segmentation (ChromHMM), single-molecule ion-channel patch-clamp current analysis, and coalescent ancestral demographic inference (Pairwise Sequentially Markovian Coalescent, PSMC).

State Space, Transition Matrix, and Emission Distributions

A discrete first-order Hidden Markov model θ = (S, Σ, π₀, A, E) is specified by five components:

1. Latent State Alphabet (S): A finite set of K = |S| hidden states S = {s₁, s₂, ..., s_K} (for example, {α-helix, β-sheet, random coil} in protein secondary structure, or {Match, Insert, Delete} in a profile alignment model).
2. Observation Alphabet (Σ) or Continuous Space: Either a discrete symbol alphabet Σ = {v₁, ..., v_M} (such as the 4 nucleotides {A, C, G, T} or 20 canonical amino acids) or a continuous vector space ℝ^d.
3. Initial State Distribution (π₀): A probability vector a_0k = P(π₁ = s_k) satisfying ∑_k=1^K a_0k = 1.
4. State Transition Probability Matrix (A): A K × K row-stochastic matrix with entries a_kl = P(π_t = s_l | π_t−1 = s_k), obeying the first-order Markov memoryless property P(π_t | π₁, ..., π_t−1) = P(π_t | π_t−1) and row normalization ∑_l=1^K a_kl = 1.
5. Emission Probability Distributions (E): For each latent state k ∈ S, a probability distribution e_k(x_t) = P(x_t | π_t = s_k) governing the observed symbol (categorical multinomial probabilities ∑_{b ∈ Σ} e_k(b) = 1, or continuous Gaussian mixture densities). Conditional on the current latent state π_t, the observation x_t is statistically independent of all previous states and observations.

Under these two conditional independence assumptions, the joint probability of observing sequence X = (x₁, ..., x_T) together with a specific latent state path π = (π₁, ..., π_T) factorizes into a product of T transitions and T emissions:

P(X, π | θ) = a_0,π1 · e_π1(x₁) · ∏_t=2^T [ a_πt−1,πt · e_πt(x_t) ]

Evaluating or training a Hidden Markov model requires solving three computational problems, each of which would require summing or maximizing over K^T = exp(T ln K) exponential state paths if evaluated by brute force, but reduces to O(T K²) polynomial time via dynamic programming.

1. Likelihood Evaluation: The Forward and Backward Algorithms

To evaluate the total marginal probability P(X | θ) = ∑_π P(X, π | θ) that model θ generated sequence X across all K^T possible hidden paths, the Forward algorithm defines the forward variable α_t(k) = P(x₁, x₂, ..., x_t, π_t = s_k | θ) as the joint probability of observing prefix x_1..t and arriving in state s_k at step t:

• Initialization (t = 1): α₁(k) = a_0k · e_k(x₁) for 1 ≤ k ≤ K.
• Induction (t = 2, ..., T): α_t(l) = e_l(x_t) · ∑_k=1^K α_t−1(k) a_kl.
• Termination: P(X | θ) = ∑_k=1^K α_T(k).

Symmetrically, the Backward algorithm defines the backward variable β_t(k) = P(x_t+1, ..., x_T | π_t = s_k, θ) as the conditional probability of observing the remaining suffix from t + 1 to T given state s_k at time t, initialized at β_T(k) = 1 and propagated backward from t = T − 1 down to 1 via:

β_t(k) = ∑_l=1^K a_kl · e_l(x_t+1) · β_t+1(l)

Because multiplying hundreds or thousands of probabilities smaller than 1 rapidly underflows IEEE 754 floating-point precision (P(X | θ) ∼ exp(−cT)), production implementations either rescale α_t(k) by 1 / ∑_j α_t(j) at every step t or accumulate log-probabilities using the numerically stable identity:

ln(exp(u) + exp(v)) = max(u, v) + ln(1 + exp(−|u − v|))

2. State Sequence Decoding: Viterbi versus Posterior Decoding

When the goal is to uncover the hidden biological structure along a sequence—such as locating the boundaries of transmembrane α-helices or coding exons—two distinct decoding criteria are used:

• Viterbi Path Decoding (Global Path Optimum): Formulated by Andrew J. Viterbi in 1967 for convolutional channel codes, the Viterbi algorithm replaces the sum in the Forward step with a maximum operator to find the single globally most probable state path π^* = arg max_π P(X, π | θ). Operating in log-space to avoid underflow, let V_t(l) denote the log-probability of the best prefix path ending in state l at step t:

V_t(l) = ln e_l(x_t) + max_{1 ≤ k ≤ K} [ V_t−1(k) + ln a_kl ] Saving the maximizing predecessor pointer ptr_t(l) = arg max_k [V_t−1(k) + ln a_kl] allows a linear O(T) traceback from arg max_k V_T(k) to recover π^*.

• Posterior (Forward–Backward) Decoding (Pointwise Marginal Optimum): When many near-optimal paths share similar probabilities, the single Viterbi path may be less informative than the marginal posterior probability γ_t(k) = P(π_t = s_k | X, θ) that position t belongs to state k, summed over all paths passing through (t, k). Combining the Forward and Backward variables yields:

γ_t(k) = [ α_t(k) · β_t(k) ] / P(X | θ) Choosing π̂_t = arg max_k γ_t(k) at each position maximizes the expected number of correctly labeled positions, and the posterior probability γ_t(k) ∈ [0, 1] provides a calibrated pointwise confidence score for each residue alignment column.

3. Parameter Estimation: The Baum–Welch (EM) Algorithm

When the latent state labels are unobserved in training data, the transition matrix A and emission matrix E cannot be estimated by direct frequency counting. Leonard Baum, Ted Petrie, George Soules, and Norman Weiss (1970) developed the Baum–Welch algorithm, a special case of the Expectation–Maximization (EM) framework later generalized by Arthur Dempster, Nan Laird, and Donald Rubin (1977):

• E-Step (Expectation): Using the current parameter estimate θ^(r), run the Forward and Backward passes to compute the posterior state occupancy γ_t(k) = P(π_t = s_k | X, θ^(r)) and the posterior transition responsibility ξ_t(k, l) = P(π_t = s_k, π_t+1 = s_l | X, θ^(r)):

ξ_t(k, l) = [ α_t(k) · a_kl · e_l(x_t+1) · β_t+1(l) ] / P(X | θ^(r)) Summing ξ_t(k, l) over t = 1..T − 1 yields the expected number of transitions Â_kl from state k to state l, and summing γ_t(k) over steps where x_t = b yields the expected emission count Ê_k(b).

• M-Step (Maximization): Normalize the expected counts to obtain updated parameters a_kl^(r+1) = Â_kl / ∑_l′ Â_kl′ and e_k^(r+1)(b) = Ê_k(b) / ∑_b′ Ê_k(b′). By Jensen's inequality and the non-negativity of Kullback–Leibler divergence, each EM iteration monotonically increases (or leaves unchanged) the marginal log-likelihood ln P(X | θ^(r+1)) ≥ ln P(X | θ^(r)) until converging to a local maximum.

Krogh–Haussler–Eddy Profile HMMs (Plan 7 Architecture)

A landmark breakthrough in bioinformatics occurred in 1994 when Anders Krogh, Michael Brown, I. Saira Mian, Kimmen Sjölander, and David Haussler introduced profile Hidden Markov models to convert a multiple sequence alignment of a homologous protein or RNA family into a position-specific probabilistic model. Whereas standard pairwise alignment uses position-invariant substitution matrices (BLOSUM62) and uniform affine gap penalties, a protein domain contains conserved active-site columns (with low Shannon entropy and strong intolerance to insertions) alternating with variable surface loops that readily tolerate insertions and deletions.

A profile HMM of consensus length L models each alignment column j ∈ {1, ..., L} with a repeating trio of states:

1. Match State (M_j): Emits a single residue according to column j's position-specific amino-acid or nucleotide probability vector e_Mj(x) (regularized with Dirichlet mixture priors to handle small sample sizes).
2. Insert State (I_j): Emits one or more extra residues between consensus columns j and j + 1 via a self-loop transition a_Ij,Ij.
3. Silent Delete State (D_j): Emits no symbol, allowing a query sequence to skip consensus column j at a position-specific gap-opening (M_j−1 → D_j) and gap-extension (D_j−1 → D_j) probability.

Implemented in Sean Eddy's HMMER software suite (using the Plan 7 architecture, which forbids direct I_j ↔ D_j transitions and sums over all alignments via SIMD-vectorized Forward scoring) and the Pfam / InterPro databases, profile HMMs detect distant evolutionary homologs across the "twilight zone" (< 20–25% sequence identity) where pairwise BLAST searches fail.

CpG Islands, Generalized Gene-Finding HMMs, and Coalescent HMMs

Beyond homology search, HMMs partition whole genomes into functional annotations:

• CpG Island and Chromatin Segmentation: In vertebrate genomes, 5-methylcytosine spontaneously deaminates to thymine, depleting CG dinucleotides genome-wide except within unmethylated promoter CpG islands. A two-regime HMM (or multi-state ChromHMM over ChIP-seq histone modifications) switches between background and island transition matrices to segment chromosomes via Viterbi or posterior decoding.
• Generalized Hidden Markov Models (GHMMs) for Gene Prediction: Because standard geometric state dwell-time distributions P(length = d) = (1 − p)p^d−1 poorly match empirical eukaryotic exon lengths and 3-base codon periodicity, eukaryotic gene finders such as GENSCAN (Chris Burge and Samuel Karlin, 1997) and Augustus use semi-Markov (generalized) HMMs in which each state emits an entire variable-length feature (5′ UTR, initial exon, donor splice site, intron, acceptor splice site, internal exon, terminal exon) conditioned on explicit length distributions.

Why do sequence analysts use both the Viterbi algorithm and the Forward algorithm on the same profile HMM?

The Viterbi algorithm finds the single best alignment path π^* between a query sequence and a profile HMM, which is required to display residue-to-column correspondences. However, when two distantly related proteins can be aligned in dozens of slightly shifted register variations around gapped loops, summing the probabilities of all valid alignments via the Forward algorithm yields a more sensitive log-odds bit score for homology detection than evaluating any single path alone.

How does the Markov assumption limit standard HMMs on RNA secondary structure?

A first-order Hidden Markov model conditions only on the immediately preceding state π_t−1 (equivalent to a regular grammar in the Chomsky hierarchy). Because RNA base pairing couples distant nucleotides i and j separated by arbitrary hairpin spans (i ≪ j), modeling nested base pairs requires stepping up from regular HMMs to Stochastic Context-Free Grammars (SCFGs) (such as the Infernal covariance models used by the Rfam database), whose Inside and CYK dynamic programming algorithms run in O(L³) rather than O(L) time.

1. Rabiner, L. R. (1989). A tutorial on hidden Markov models and selected applications in speech recognition. Proceedings of the IEEE, 77(2), 257–286. https://doi.org/10.1109/5.18626
2. Baum, L. E., Petrie, T., Soules, G., & Weiss, N. (1970). A maximization technique occurring in the statistical analysis of probabilistic functions of Markov chains. The Annals of Mathematical Statistics, 41(1), 164–171. https://doi.org/10.1214/aoms/1177697196
3. Krogh, A., Brown, M., Mian, I. S., Sjölander, K., & Haussler, D. (1994). Hidden Markov models in computational biology: Applications to protein modeling. Journal of Molecular Biology, 235(5), 1501–1531. https://doi.org/10.1006/jmbi.1994.1104
4. Eddy, S. R. (1998). Profile hidden Markov models. Bioinformatics, 14(9), 755–763. https://doi.org/10.1093/bioinformatics/14.9.755
5. Burge, C., & Karlin, S. (1997). Prediction of complete gene structures in human genomic DNA. Journal of Molecular Biology, 268(1), 78–94. https://doi.org/10.1006/jmbi.1997.0951