Binary stars with pulsating components are ubiquitous. Eclipsing binary stars among them are of special interest. Given that eclipsing binaries provide us with accurate fundamental stellar parameters for their components (most notably their masses and radii, along with temperatures and luminosities), finding and studying intrinsically pulsating components allows us to study principal physical processes behind pulsations. Due to mostly the lack of tools to study pulsations and eclipsing binary signatures simultaneously, researchers have resorted to sequential, iterative fitting: they would first fit a binary light curve, then subtract it and fit a pulsation model to the residuals; then correct for correlated residuals and repeat the process until the two signals were disentangled. Of course, when there is commensurability between the two signals, either incidental or a consequence of the physical processes, this procedure fails to produce robust results. Instead, we need to model binarity and pulsations simultaneously. In this contribution we present the latest advancement of the PHOEBE eclipsing binary modeling code that includes support for pulsations. We showcase two synthetic cases: one where the binary signal dominates and the other where the pulsational signal dominates, and then we focus on a real-world example observed by Kepler: KIC 4150611, a quintuple system with the brightest A-type star exhibiting delta-Sct pulsations, and one of the fainter G-type stars featuring g-mode pulsations. This research is sponsored by the NSF grant #1909109, which we gratefully acknowledge.