Recently we succeeded, by introducing an interaction between the gluon mean field (presented by the a function of the Polyakov loop) and quarks, to reproduce the lattice equation of state for zero chemical potential with the Polyakov-Nambu-Jona-Lasinio model. Also, entropy density, interaction measure, energy density and the speed of sound are quite nicely reproduced. Even the first coefficient of the Taylor expansion of the lattice data with respect to the chemical potential is in the error bars of the lattice calculations. These findings are of great importance for future studies of heavy ion reactions because the Polyakov-Nambu-Jona-Lasinio model can be extended to finite chemical potentials (where lattice calculations are not possible) without introducing any new parameter. In addition, it shows at large chemical potentials a first order phase transition. It provides therefore a basis for theoretical studies in the energy range of the future FAIR and NICA facilities where one expects that heavy ion collisions are characterized by a large chemical potential. It may also serve as a equation of state for gravitational wave studies.