The coronavirus disease 2019 (COVID-19) pandemic placed epidemic modeling at the forefront of worldwide public policy making. Nonetheless, modeling and forecasting the spread of COVID-19 remain a challenge. This talk begins with a review of the historical use of epidemic models and addresses the challenges of choosing a model in the early stages of a worldwide pandemic. The spread of COVID-19 has illustrated the heterogeneity of disease spread at different population levels. With finite size populations, chance variations may lead to significant departures from the mean. In real-life applications, finite size effects arise from the variance of individual realizations of an epidemic course about its fluid limit. I will illustrate how to model this variance with a martingale formulation consisting of a deterministic and a stochastic component, along with estimates for the size of the variance compared to real world data and simulations. Another cause of heterogeneities is the differing attitudes at the population level for control measures such as mask-wearing and physical distancing. Often, individuals form opinions about their behaviors from social network opinions. I will show some results from a two-layer multiplex network for the coupled spread of a disease and conflicting opinions. We model each process as a contagion. We develop approximations of mean-field type by generalizing monolayer pair approximations to multilayer networks; these approximations agree well with Monte Carlo simulations for a broad range of parameters and several network structures. We find that lengthening the duration that individuals hold an opinion may help suppress disease transmission, and we demonstrate that increasing the cross-layer correlations or intra-layer correlations of node degrees may lead to fewer individuals becoming infected with the disease.