In the recent COVID-19 pandemic, mathematical modeling constitutes an important tool to evaluate the prospective effectiveness of non-pharmaceutical interventions (NPIs) and to guide policy-making. Most research is, however, centered around characterizing the epidemic based on point estimates like the average infectiousness or the average number of contacts. In this work, we use stochastic simulations to investigate the consequences of a population's heterogeneity regarding connectivity and individual viral load levels. Therefore, we translate a COVID-19 ODE model to a stochastic multi-agent system. We use contact networks to model complex interaction structures and a probabilistic infection rate to model individual viral load variation. We observe a large dependency of the dispersion and dynamical evolution on the population's heterogeneity that is not adequately captured by point estimates, for instance, used in ODE models. In particular, models that assume the same clinical and transmission parameters may lead to different conclusions, depending on different types of heterogeneity in the population. For instance, the existence of hubs in the contact network leads to an initial increase of dispersion and the effective reproduction number, but to a lower herd immunity threshold (HIT) compared to homogeneous populations or a population where the heterogeneity stems solely from individual infectivity variations.