We present a novel hybrid finite element (FE) - spectral boundary integral (SBI) scheme that enables efficient simulation of earthquake cycles. This combined FE-SBI approach captures the benefits of finite elements in modelling problems with nonlinearities, as well as the computational superiority of SBI. The domain truncation enabled by this scheme allows us to utilize high-resolution finite elements discretization to capture inhomogeneities or complexities that may exist in a narrow region surrounding the fault. Combined with an adaptive time stepping algorithm, this framework opens new opportunities for modeling earthquake cycles with high-resolution fault zone physics. In this initial study, we consider a two dimensional (2-D) anti-plane model with a vertical strike-slip fault governed by rate and state friction in the quasi-dynamic limit under the radiation damping approximation. The proposed approach is first verified using the benchmark problem BP-1 from the Southern California Earthquake Center (SCEC) sequence of earthquake and aseismic slip (SEAS) community verification effort. The computational framework is then utilized to model the earthquake sequence and aseismic slip of a fault embedded within a low-velocity fault zone (LVFZ) with different widths and compliance levels. Our results indicate that sufficiently compliant LVFZs contribute to the emergence of sub-surface events that fail to penetrate to the free surface and may experience earthquake clusters with nonuniform inter-seismic time. Furthermore, the LVFZ leads to slip rate amplification relative to the homogeneous elastic case. We discuss the implications of our results for understanding earthquake complexity as an interplay of fault friction and bulk heterogeneities.