Black holes are the most massive, energetic, and brightest sources in the Universe. Massive black holes (MBHs), which are up to many billions of times more massive than our sun, are also believed to shape the fate of entire galaxies. How these behemoths form and grow in the early Universe, and how they couple to their host galaxies, remains a mystery. Following galaxy mergers, binary MBHs can form which are expected to produce low-frequency gravitational waves (GW). The simulation framework that I have developed represents the most advanced and comprehensive suite of models available for MBH binaries and their host galaxies. In this talk, I will describe my predictions for both GW and electromagnetic signatures, and highlight the unique discovery space opened up by their combination as multimessenger signals. In the next decade, my research program is poised to revolutionize our understanding of MBH formation in the early Universe, and their co-evolution with galaxies over cosmic time. In the NANOGrav collaboration, we have recently identified a nanohertz signal that is consistent with the first hints of a low-frequency GW signal. I will discuss work that I am leading as the chair of NANOGrav's astrophysics working group in developing a comprehensive framework for predicting and interpreting the upcoming GW and electromagnetic signatures. Finally, I will describe the new generation of cosmological simulations that we are developing to model high-redshift MBHs that will soon be observable by JWST and eventually with the LISA mission.