28 Hydrophobic Contact Maximization on a Lattice

Source and adaptation

Based on Protein folding, Section 12.28 of Williams (2013). New eight-residue, two-dimensional square-lattice HP toy model. Explicitly replaces the source 50-residue non-lattice example. The results below solve the stated instance and are not presented as the numerical answer to an unmodified textbook problem.

28.1 Problem Statement

Fold the sequence H–P–H–H–P–P–H–H on a two-dimensional square lattice. Consecutive residues occupy adjacent lattice sites, no site may contain two residues, and a contact is a pair of H residues at Manhattan distance one that are not consecutive along the chain. Maximize the number of contacts, equivalent to minimizing energy of −1 per contact. Fix residue 1 at (0,0) and residue 2 at (1,0) to remove translation and rotation. Reflections remain allowed. This model illustrates a discrete optimization idea; it is not a realistic protein-structure predictor.

28.2 Model Formulation

Enumerate all self-avoiding lattice conformations of this eight-residue chain with the first bond fixed. Let \(K\) be this finite set and \(c_k\) the number of nonbonded H–H contacts in conformation \(k\).

\[\sum_{k\in K}y_k=1,\qquad y_k\in\{0,1\},\qquad\max\sum_{k\in K}c_ky_k.\]

The enumeration enforces chain connectivity and excluded volume before the optimization is built. The one-hot model selects a configuration; it is an exact but deliberately elementary formulation for this small lattice problem. It is not a scalable residue-placement MILP and makes no claim about the source non-lattice model.

All decision variables and units refer to the problem statement above. Continuous variables are nonnegative unless explicitly stated otherwise; binary and integer domains are specified in the equations and code.

28.3 Pyomo Implementation

The following Python implementation uses Pyomo and HiGHS. Run from the workbook root so that the models package and shared helpers are importable. The shared solver and audit functions check optimal termination before loading values, then verify every active constraint, variable bound, and integer domain. Any imported earlier-chapter model supplies the data and balances already explained there.

The full source is problem28.py. To solve and print this problem independently:

~/.venvs/optim/bin/python -m models.common 28

import pyomo.environ as pyo
from models.common import frame

SEQUENCE = "HPHHPPHH"
CONFORMATIONS = []


def enumerate_paths(path):
    if len(path) == len(SEQUENCE):
        CONFORMATIONS.append(tuple(path))
        return
    x, y = path[-1]
    for dx, dy in [(1, 0), (-1, 0), (0, 1), (0, -1)]:
        point = (x + dx, y + dy)
        if point not in path:
            enumerate_paths(path + [point])


enumerate_paths([(0, 0), (1, 0)])


def contacts(path):
    return [
        (i, j)
        for i in range(len(path))
        for j in range(i + 2, len(path))
        if SEQUENCE[i] == SEQUENCE[j] == "H"
        and sum(abs(path[i][k] - path[j][k]) for k in [0, 1]) == 1
    ]


def build():
    m = pyo.ConcreteModel()
    m.K = pyo.RangeSet(0, len(CONFORMATIONS) - 1)
    m.y = pyo.Var(m.K, domain=pyo.Binary)
    m.one = pyo.Constraint(expr=sum(m.y[k] for k in m.K) == 1)
    m.obj = pyo.Objective(
        expr=sum(len(contacts(path)) * m.y[k] for k, path in enumerate(CONFORMATIONS)),
        sense=pyo.maximize,
    )
    return m


def chosen(m):
    return CONFORMATIONS[next(k for k in m.K if pyo.value(m.y[k]) > 0.5)]


def check(m):
    path = chosen(m)
    assert len(set(path)) == 8
    assert all(
        sum(abs(a[k] - b[k]) for k in [0, 1]) == 1 for a, b in zip(path, path[1:])
    )
    assert len(contacts(path)) == max(map(lambda p: len(contacts(p)), CONFORMATIONS))


def tables(m):
    path = chosen(m)
    return {
        "conformation": frame(
            [[i + 1, SEQUENCE[i], *point] for i, point in enumerate(path)],
            ["Residue", "Type", "x", "y"],
        ),
        "contacts": frame(
            [[i + 1, j + 1] for i, j in contacts(path)], ["Residue i", "Residue j"]
        ),
        "enumeration": frame(
            [[len(CONFORMATIONS), pyo.value(m.obj)]],
            ["Conformations tested", "Maximum contacts"],
        ),
    }


def plot(m):
    path = chosen(m)
    return (
        [str(i + 1) for i in range(8)],
        [sum(i in pair for pair in contacts(path)) for i in range(8)],
        "Nonbonded H contacts per residue",
    )

# Solve, audit constraints, and run domain-specific checks.
from models.common import solve, audit
model = solve(build())
audit(model)
check(model)
for name, result_table in tables(model).items():
    print(name)
    print(result_table.to_string(index=False))

28.4 Optimal Solution

The solver reports optimal termination. The objective is 3 nonbonded hydrophobic contacts (maximize). The largest violation across active constraints, variable bounds, and integer domains is 0.00e+00 in the corresponding model units. The problem-specific checks also pass. These checks establish numerical consistency with the stated model, not the validity of its assumptions for a real facility.

28.4.1 Conformation

Residue	Type	x	y
1	H	0	0
2	P	1	0
3	H	1	-1
4	H	1	-2
5	P	1	-3
6	P	0	-3
7	H	0	-2
8	H	0	-1

28.4.2 Contacts

Residue i	Residue j
1	8
3	8
4	7

28.4.3 Enumeration

Conformations tested	Maximum contacts
543	3

Figure 28.1: One optimal lattice fold. Amber residues are hydrophobic, teal residues are polar, and dashed lines identify nonbonded H contacts.

Tables round numerical values for reading; feasibility checks use the original solver values. Multiple optimal decisions may exist. Machine-readable result records the solver status and package versions. Figure-generation code is in figures/workbook.py.

28.5 Brief Discussion

Enumeration tests 543 self-avoiding conformations. The best has 3 nonbonded H–H contacts, giving energy -3 in the stated contact-energy units.

Excluding covalently adjacent residues prevents every H–H bond from being counted as a folding benefit. Several conformations can share the best contact score. Enumeration grows rapidly with chain length; the small example makes the feasible set and contact definition auditable. Experiment: change one polar residue to H and compare both the optimal score and geometry.