We present a bi-orthogonal approach for modeling the response of localized electromagnetic resonators using quasinormal modes, which represent the natural, dissipative eigenmodes of the system with complex frequencies. For many problems of interest in optics and nanophotonics, the quasinormal modes constitute a powerful modeling tool, and the bi-orthogonal approach provides a coherent, precise, and accessible derivation of the associated theory, enabling an illustrative connection between different modeling approaches that exist in the literature.