  An orthogonal product basis (OPB) of a finite-dimensional Hilbert space
$H=H_1\otimes H_2\otimes\cdots\otimes H_n$ is an orthonormal basis of $H$
consisting of product vectors $x_1\otimes x_2\otimes\cdots\otimes x_n$. We show
that the problem of classifying the OPBs of an $n$-qubit system can be reduced
to a purely combinatorial problem. We solve this combinatorial problem in the
case of four qubits and obtain 33 multiparameter families of OPBs. Each OPB of
four qubits is equivalent, under local unitary operations and qubit
permutations, to an OPB belonging to at least one of these families.
