Theoretical analysis and results of calculations are put forward to interpret the energies -epsilon(k) of the occupied Kohn-Sham (KS) orbitals as approximate but rather accurate relaxed vertical ionization potentials (VIPs) I-k. Exact relations between epsilon(k) and I-k are established with a set of linear equations for the epsilon(k), which are expressed through I-k and the matrix elements epsilon(k)(resp) of a component of the KS exchange-correlation (xc) potential v(xc), the response potential v(resp). Although -I-k will be a leading contribution to epsilon(k), other I-jnot equalk do enter through coupling terms which are determined by the overlaps between the densities of the KS orbitals as well as by overlaps between the KS and Dyson orbital densities. The orbital energies obtained with "exact" KS potentials are compared with the experimental VIPs of the molecules N-2, CO, HF, and H2O. Very good agreement between the accurate -epsilon(k) of the outer valence KS orbitals and the corresponding VIPs is established. The average difference, approaching 0.1 eV, is about an order of magnitude smaller than for HF orbital energies. The lower valence KS levels are a few eV higher than the corresponding -I-k, and the core levels some 20 eV, in agreement with the theoretically deduced upshift of the KS levels compared to -I-k by the response potential matrix elements. Calculations of 64 molecules are performed with the approximate v(xc) obtained with the statistical averaging of (model) orbitals potentials (SAOP) and the calculated epsilon(k) are compared with 406 experimental VIPs. Reasonable agreement between the SAOP -epsilon(k) and the outer valence VIPs is found with an average deviation of about 0.4 eV.