A coarse-grained model has been proposed recently in order to describe plastic properties of glassy polymers with space resolution the scale xi of dynamical heterogeneities. This model allows for describing plastic flow by assuming that the elastic energy stored at the scale of dynamical heterogeneities reduces by a similar amount the free energy barriers for alpha-relaxation. The aim of this article is to consider in more details the evolution of the distribution of relaxation times under plastic deformation. This is achieved by taking explicitly into account the so-called facilitation mechanism during plastic flow introduced by Merabia and Long [Eur. Phys. J. E 2002, 9, 195; J. Chem. Phys. 2006, 125, 234901]. This mechanism, which allows for calculating the scale of dynamical heterogeneities, is key for explaining temporal asymmetry regarding the dynamical evolution of a glassy polymer upon cooling or upon heating, as observed experimentally by Kovacs. Upon heating, fast regions appear first and melt the slowest ones. We propose that this same process is at work upon stretching a polymer material. Some subunits get accelerated due to deformation, and their growing number allows for melting the slowest subunits. We show that this facilitation mechanism is responsible for the narrowing of the relaxation time spectrum observed e.g. during plastic strain by using optical probe dynamics measurements by Ediger and co-workers. The physical model is solved numerically in 3D by a method similar to dissipative particle dynamics, in which the strain and stress fields, as well as the local relaxation times, are calculated with a resolution the size of dynamical heterogeneities.