Multiple zigzag chains Z(m, n) of length n and width m constitute an important class of regular graphene flakes of rectangular shape. The physical and chemical properties of these basic pericondensed benzenoids can be related to their various topological invariants, conveniently encoded as the coefficients of a combinatorial polynomial, usually referred to as the ZZ polynomial of multiple zigzag chains Z(m, n). The current study reports a novel method for determination of these ZZ polynomials based on a hypothesized extension to John–Sachs theorem, used previously to enumerate Kekulé structures of various benzenoid hydrocarbons. We show that the ZZ polynomial of the Z(m, n) multiple zigzag chain can be conveniently expressed as a determinant of a Toeplitz (or almost Toeplitz) matrix of size ⌈ m/ 2 ⌉×⌈m/2 ⌉ consisting of simple hypergeometric polynomials. The presented analysis can be extended to generalized multiple zigzag chains Zk (m, n), i.e., derivatives of Z(m, n) with a single attached polyacene chain of length k. All presented formulas are accompanied by formal proofs. The developed theoretical machinery is applied for predicting aromaticity distribution patterns in large and infinite multiple zigzag chains Z(m, n) and for computing the distribution of spin densities in biradical states of finite multiple zigzag chains Z(m, n).