We have developed a necklace model of hydrophobic polyelectrolytes in which the necklace structure consisting of polymeric globules (beads) connected by extended sections of the chain (strings of monomers) appears as a result of the counterion condensation and is caused by the balance of the correlation-induced attraction of condensed counterions to charged monomers and electrostatic repulsion between uncompensated charges. The size of the beads increases with polymer concentration while their number per chain decreases. We predict coexistence of necklaces with different number of beads on a polymer backbone at any polymer concentration. To test this necklace model, we performed molecular dynamics simulations of polyelectrolyte chains with degree of polymerization N varying from 25 to 373 and with fraction of charged monomers f = (3)\(1), (2)\(1), and 1 in poor solvent conditions for polymer backbone. The observed concentration dependence of the bead size supports the assumption of the counterion condensation origin of the necklace structure. The overlap concentration is almost independent of the degree of polymerization for weakly charged chains (f = (3)\(1)). For strongly charged chains with f = 1 the overlap concentration follows the normal N dependence observed for polyelectrolyte solutions in circle minus and good solvent regimes for polymer backbone. In semidilute solutions the correlation length of fully charged chains is inversely proportional to the square root of polymer concentration. The osmotic coefficient of the solutions of polyelectrolytes in poor solvent conditions for polymer backbone exhibits nonmonotonic concentration dependence in agreement with two-zone model predictions. It decreases with increasing polymer concentration in dilute solutions, while it is an increasing function of polymer concentration in the semidilute regime. The Kratky plot of the chain form factor is in excellent agreement with the neutron scattering experiments.