In this chapter, we overview recent developments of a simulation framework capable of capturing the highly nonequilibrium physics of the strongly coupled electron and phonon systems in quantum cascade lasers (QCLs). In midinfrared (mid-IR) devices, both electronic and optical phonon systems are largely semiclassical and described by coupled Boltzmann transport equations, which we solve using an efficient stochastic technique known as ensemble Monte Carlo. The optical phonon system is strongly coupled to acoustic phonons, the dominant carriers of heat, whose dynamics and thermal transport throughout the whole device are described via a global heat-diffusion solver. We discuss the roles of nonequilibrium optical phonons in QCLs at the level of a single stage , anisotropic thermal transport of acoustic phonons in QCLs, outline the algorithm for multiscale electrothermal simulation, and present data for a mid-IR QCL based on this framework.