Seismic isolation has become increasingly popular due to the excellent performance observed in buildings during past earthquakes. The seismic response of isolated structures to earthquake input is strongly controlled by the force-deformation constitutive behavior of the isolators. High Damping Rubber Bearings (HDRBs) are one of the most widely manufactured and used isolation systems in practice. Because of the large shear flexibility, the stress-strain constitutive behavior of the elastomeric material controls the overall behavior of the device. Thus, the behavior of HDRBs is highly nonlinear and characterized by the same phenomena as the elastomeric material, which is challenging to model analytically. Consequently, this article describes a simple but sufficiently accurate numerical model for simulating the bi-directional shear behavior of HDRBs under large shear deformations. A brief description of the experimental test data of HDRBs is presented, as well as the importance of including the observed phenomena in the numerical model proposed. The mathematical formulation of the proposed model is summarized, based on a hyperelastic spring and dissipative element connected in parallel. The former considers anisotropic degradation (scragging and Mullins effect), while the latter includes the isotropic hardening phenomenon. The proposed model was validated using experimental results of bi-directional shear tests of cylindrical disks of high damping natural rubber, unidirectional shear cyclic tests, and a bi-directional shear loading history applied to an HDRBs. Since the proposed model is capable of simulating the bi-directional cyclic behavior of HDRBs by considering anisotropic degradation instead of the classical isotropic behavior, which improves the numerical predictions, it can be used in dynamic analyses of buildings isolated with HDRBs.