HP Lattice Folding: Simplified Protein Structure Prediction

Proteins fold into specific 3D structures, and that structure determines their function. Predicting structure from sequence is one of the grand challenges of computational biology. While AlphaFold has revolutionized the field, the simplified HP lattice model remains a useful testbed for optimization algorithms.

The HP model strips protein folding down to its combinatorial essence: a sequence of hydrophobic (H) and polar (P) amino acids, a discrete lattice, and an energy function that counts contacts. Simple rules, NP-hard optimization.

The HP Model

Ken Dill proposed the HP model in the 1980s as a minimal model of protein folding. The idea: proteins fold to bury hydrophobic residues in a core, away from water. Polar residues end up on the surface.

In the 2D lattice version:

• The protein sequence is a string like “HPHPPHHPHP”
• Each amino acid occupies one cell on a square lattice
• Adjacent residues in the sequence must be adjacent on the lattice
• Self-intersection is forbidden
• Energy = -(number of non-sequential H-H contacts)

A “contact” means two H residues are adjacent on the lattice but not adjacent in the sequence. More contacts means lower energy, which is better.

This is a drastic simplification. Real proteins have 20 amino acid types, fold in 3D continuous space, and have complex energy functions involving electrostatics, hydrogen bonding, and solvent effects. But the HP model captures something essential: the hydrophobic effect drives folding, and finding the optimal fold is hard.

Why It’s Hard

For a sequence of length N, the genome has N-1 moves. Each move is one of three choices: forward, left, or right (relative to the current direction). That’s 3^(N-1) possible folds.

But most of those are invalid because they self-intersect. The valid folds form a complex, irregular subset of the search space. You can’t just enumerate them efficiently.

Even among valid folds, the fitness landscape is rugged. Similar move sequences can produce very different structures. A single move change might either slightly adjust the shape or cause a collision that invalidates the whole fold.

For longer sequences, the number of local optima grows rapidly. Gradient-based methods don’t apply because the space is discrete. Exhaustive search is impossible. This is exactly where genetic algorithms are useful.

Encoding the Problem

ParGA uses real-valued genomes, but HP folding needs discrete moves. The solution: encode moves as continuous values and discretize during fitness evaluation.

def decode_moves(genes: np.ndarray) -> list[int]:
    """Convert continuous genes to discrete moves (0, 1, 2)."""
    moves = np.clip(np.floor(genes * 3).astype(int), 0, 2)
    return moves.tolist()

Genes in [0, 1) map to moves 0, 1, or 2. Crossover and mutation work on the continuous values, which then get discretized. This lets us use standard real-valued GA operators on what is really a discrete problem.

Folding Algorithm

Given a move sequence, the folding is deterministic:

def fold_protein(sequence: str, moves: list[int]) -> list[tuple[int, int]] | None:
    positions = [(0, 0)]  # First amino acid at origin
    occupied = {(0, 0)}
    direction = 0  # Start facing "up"

    for move in moves:
        # Apply turn
        if move == Move.LEFT:
            direction = (direction - 1) % 4
        elif move == Move.RIGHT:
            direction = (direction + 1) % 4
        # FORWARD keeps same direction

        # Move in current direction
        dx, dy = DIRECTIONS[direction]
        new_pos = (positions[-1][0] + dx, positions[-1][1] + dy)

        # Check for self-intersection
        if new_pos in occupied:
            return None

        positions.append(new_pos)
        occupied.add(new_pos)

    return positions

The function returns None if self-intersection occurs. Invalid folds get heavily penalized in fitness.

Counting Contacts

The energy calculation counts non-sequential H-H contacts:

def count_hh_contacts(sequence: str, positions: list[tuple[int, int]]) -> int:
    contacts = 0
    n = len(sequence)
    pos_to_idx = {pos: i for i, pos in enumerate(positions)}

    for i in range(n):
        if sequence[i] != "H":
            continue

        x, y = positions[i]
        for dx, dy in DIRECTIONS:
            neighbor = (x + dx, y + dy)
            if neighbor in pos_to_idx:
                j = pos_to_idx[neighbor]
                if abs(i - j) > 1 and sequence[j] == "H":
                    contacts += 1

    return contacts // 2  # Each contact counted twice

The position lookup dictionary makes neighbor checks O(1). For longer sequences, this matters.

Fitness Function

Putting it together:

def hp_fitness(genes: np.ndarray, sequence: str) -> float:
    moves = decode_moves(genes)
    positions = fold_protein(sequence, moves)

    if positions is None:
        return -1000.0  # Invalid fold penalty

    contacts = count_hh_contacts(sequence, positions)
    return float(contacts)

Higher fitness is better, so we return contact count directly. Invalid folds get -1000, which is bad enough to be eliminated but not so extreme that it dominates the fitness landscape.

Running the Optimization

from parga import minimize

sequence = "HPHPPHHPHPPHPHHPPHPH"  # 20 residues

def fitness(genes):
    return -hp_fitness(genes, sequence)

result = minimize(
    fitness,
    genome_length=len(sequence) - 1,  # N-1 moves
    bounds=(0.0, 1.0),
    population_size=200,
    generations=300,
    mutation_rate=0.1,  # Higher for discrete problems
    seed=42,
)

best_moves = decode_moves(result.best_genes())
best_positions = fold_protein(sequence, best_moves)
best_contacts = count_hh_contacts(sequence, best_positions)

The mutation rate is higher than for continuous problems. With discretization, small continuous changes often have no effect (they map to the same move), so we need more aggressive mutation to explore effectively.

Benchmark Sequences

The literature has standard test sequences:

Longer sequences are harder because the search space grows exponentially while valid folds become rarer.

Results

Five standard sequences, island model, three independent restarts each:

The obvious expectation is that this degrades with length, since the number of self-avoiding walks grows exponentially. It mostly does, and then seq24 breaks the pattern by being solved exactly while the shorter seq20 is not.

That is worth sitting with, because it is the more useful lesson. seq20 and seq24 differ by four residues and the longer one is easier, so the difficulty is a property of the sequence rather than of its length. What makes an HP sequence hard is how its H residues are distributed: a sequence whose hydrophobic residues can be brought together by several different folds gives the search a gradient to follow, and one whose optimum requires a single specific compact arrangement does not.

The search space size tells you how long the enumeration would take. It does not tell you whether the fitness function will lead you anywhere.

The benchmark checks itself

Small instances have a property worth exploiting: they are exhaustively enumerable, so you can find out what the answer actually is rather than trusting a reference value.

seq10 is ten residues, which is 3^9 = 19,683 move encodings, of which 12,201 fold without self-intersecting. That is nothing to enumerate:

best = -1
for moves in itertools.product([0, 1, 2], repeat=len(seq) - 1):
    coords = fold_protein(seq, list(moves))
    if coords is not None:
        best = max(best, count_hh_contacts(seq, coords))

Running that over the whole space gives the true optimum for the sequence, which is a fact about the problem rather than an opinion about the optimizer. If your GA returns that number, it is done, and any percentage reported against a different reference is telling you about the reference.

I would recommend that check on any small instance in any benchmark suite. It costs a few seconds, it is unambiguous, and it converts “the optimizer got 50%” into either “the optimizer is leaving something on the table” or “the optimizer is finished” without any interpretation in between.

Where the method stops

The two sequences it does not solve, seq20 at 78% and seq48 at 78%, fail in different ways. seq48 is a budget problem: 133 seconds, a space of roughly 3^47 encodings, and a population of a few hundred is sampling an unimaginably small fraction of it. More compute would help.

seq20 is not a budget problem. It is twenty residues, the runs converge quickly and repeatedly to 7 or 8 contacts, and more generations do not move it. That is the signature of a landscape where the optimum sits in a basin the search does not reach from where crossover and mutation can take it, and the fix is a different operator or a different encoding rather than a bigger budget.

Telling those two failures apart is most of the practical skill in applying a GA. Both report the same thing: a converged run below the known optimum.

What HP Folding Teaches

Several lessons from this problem:

Discrete spaces via continuous encoding: The gene-to-move mapping lets us use real-valued GAs on discrete problems. It’s not always the best approach, but it’s widely applicable.

Constraint handling via penalty: Invalid folds get penalized rather than repaired. This is simpler than trying to constrain the search space directly.

Problem-specific mutation rates: The 10% mutation rate would be too high for most continuous problems but works well here because discretization dampens the effect of small changes.

Search space structure matters: The geometry of valid folds is complex. Some move sequences are “close” in Hamming distance but produce completely different structures. The GA doesn’t know this, but it works anyway because the population maintains diversity.

Relation to Real Protein Folding

The HP lattice model is a toy. Real protein folding involves:

• 20 amino acid types with different properties
• 3D continuous space, not 2D lattice
• Backbone angles (phi, psi) as degrees of freedom
• Complex energy functions with many terms
• Solvent effects and temperature dependence

AlphaFold and related methods use deep learning on massive datasets of known structures. They don’t “fold” proteins so much as predict what fold a sequence will adopt based on patterns in similar sequences.

But the HP model captures something real: hydrophobic collapse is a key driving force in folding, and the combinatorial complexity of finding optimal configurations persists even in simplified models. If your algorithm can solve HP folding, it can probably handle other combinatorial optimization problems too.

Next: optimizing rocket motor grain geometry, where the physics gets more complex but the optimization strategy remains similar.