Circulation problems are not just a generalization of max-flow, there is a reduction backwards as well. Suppose we have some directed graph $G = (V, E)$ with edge costs, capacities, and lower bounds.
Any edge $u \to v$ in $G$ with we can replace with two nodes $s, t$ and two edges $s \to v$ and $u \to t$ where one of the edges has the original cost, capacities, and lower bounds and the other is free and unlimited. Call this graph $G'(e)$, where $e = u\to v$ is the edge that was replaced.
Then if a flow with a certain cost exists in $G'(\cdot)$, it must also exist as a circulation in $G$ with the same cost. Vice versa, if a circulation exists in $G$ and it uses edge $u \to v$, then that flow also exists in $G'(u\to v)$ with the same cost.
Therefore to solve the circulation problem we can pick an arbitrary edge $e$, calculate $G'(e)$ and use a traditional network flow algorithm to find the optimal flow. By the traditional arguments, this optimal flow is integral. We then pick another edge (avoiding edges that were part of a previous optimal flow) and repeat, keeping the best solution, until no more unknown edges are left.
Since in the worst case this adds a factor of $|E|$ to the complexity of the polynomial complexity, this is still polynomial. And of course the optimum from all integral flows found is itself integral.
To handle the lower bounds on the edges of $G'$, one can either note that the linear programming constraint matrix is unimodular (see these MIT lecture notes), from which it follows that there exists an integral solution if there is any solution; or one can use a standard reduction to eliminate the lower bounds (see, e.g., Ahuja et al, Network Flows, page 39) and then solve with a standard algorithm for network flow.