The debate regarding the possible molecular origins of the Mullins effect has been ongoing since its discovery. Molecular dynamics (MD) simulations were carried out to elucidate the underlying mechanism of the Mullins effect. For the first time, the key characteristics associated with the Mullins effect, including (a) the majority of stress softening occurring in the first stretch, (b) continuous softening with stress increase, (c) a permanent set, and (d) recovery with heat treatment, are captured by molecular modeling. It is discovered that the irreversible disentanglement of polymer chains is physically sufficient to interpret these key characteristics, providing molecular evidence for this long-controversial issue. Our results also reveal that filled polymers exhibit three distinct regimes, i.e., the polymer matrix, the interface, and the filler. When subjected to external strain, the polymer matrix suffers from excess deformation, indicating strong heterogeneity within the filled polymer, which offers molecular insight for the formulation of physics-based constitutive relations for filled polymers.