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.