A model operator approach to calculations of the QED corrections to energy levels in relativistic many-electron atomic systems is developed. The model Lamb shift operator is represented by a sum of local and nonlocal potentials which are defined using the results of ab initio calculations of the diagonal and nondiagonal matrix elements of the one-loop QED operator with H-like wave functions. The model operator can be easily included in any calculations based on the Dirac-Coulomb-Breit Hamiltonian. Efficiency of the method is demonstrated by comparison of the model QED operator results for the Lamb shifts in many-electron atoms and ions with exact QED calculations.