Coherent structures in turbulent flows arise from instabilities during transition or from the dynamics of the turbulent mean flow. At realistic flow regimes, direct numerical simulation of the Navier--Stokes equations remains computationally prohibitive, even on modern high-performance computing platforms. As a result, turbulence modeling relies on reduced-order or lower-fidelity approaches that parameterize unresolved small-scale effects through closure models. In turbulent regimes, frequency-domain formulations of the Navier--Stokes equations provide a computationally efficient alternative to time-domain simulations. Using a Fourier--Galerkin projection, periodic solutions can be obtained in the form of the Harmonic-Balanced Navier--Stokes (HBNS) equations. In this work, we develop physics-informed neural networks (PINNs) to model the nonlocal closure terms arising from frequency-truncated HBNS systems. The HBNS equations are embedded as soft constraints within the loss function, enabling the inference of high-frequency corrections induced by turbulent interactions while preserving the underlying physical structure. The proposed framework is applied to turbulent flow past a circular cylinder at Reynolds numbers $\mathrm{Re}=3900$ using direct numerical simulation data and $\mathrm{Re}=5000$ using experimental measurements. The learned closure models are validated at previously unseen Reynolds numbers, demonstrating strong generalization capability across flow regimes. We further assess model accuracy by systematically varying the number of retained temporal harmonics in the turbulent regime. Finally, we investigate the interpretability of the inferred closures using symbolic regression, enabling gray-box modeling and the discovery of governing expressions that can be directly compared with existing turbulence modeling approaches.