Introduction: Agent-based modeling frequently used modeling method for multi-scale mechanistic modeling. However, the same properties that make agent-based models (ABMs) well suited to representing biological systems also present significant challenges with respect to their construction and calibration, particularly with respect to the large number of free parameters often present in these models. The challenge of dealing with parameters is further exacerbated due to the fact that a great deal of phenotypic and clinical heterogeneity can be attributed to intrinsic genetic/epigenetic variation manifesting as functional parameter variation. As a result, effectively representing biological populations with ABMs requires dealing with very large multi-dimensional parameter spaces that exponentially increase the computational demands for their use. We have proposed that various machine learning (ML) and evolutionary computing approaches (such as genetic algorithms (GAs)) can be used to more effectively and efficiently deal with parameter space characterization; the current work applies GAs to the challenge of calibrating a complex ABM to a specific data set in a fashion that preserves the parameter spaces required to deal with biological heterogeneity. Methods: This project uses a GA to fit a previously validated ABM of acute systemic inflammation, the Innate Immune Response ABM (IIRABM) to clinical time series data of systemic cytokine levels. The genome for the GA is a vector generated from the IIRABM’s Model Rule Matrix (MRM), which is a matrix representation of not only the constants/parameters associated with the IIRABM’s cytokine interaction rules, but also the existence of rules themselves. Capturing heterogeneity is accomplished by a fitness function that incorporates the sample value range (“error bars”) of the clinical data. Results: The GA-enabled parameter space exploration resulted in a set of putative MRM parameterizations which closely (though not perfectly) match the cytokine time course data used to design the fitness function. The number of non-zero elements in the MRM increases significantly as the model parameterizations evolve towards a fitness function minimum, transitioning from a sparse to a dense matrix. This results in a model structure that more closely resembles (at a superficial level) the structure of data generated by a standard differential gene expression experimental study, in that there are a small number of powerful causative correlations and a much larger number of weaker/less significant (individually) connections. Conclusion: We present an HPC-enabled evolutionary computing approach that utilizes a GA to calibrate a complex ABM to clinical data while preserving biological heterogeneity. The integration of machine learning/evolutionary computing, HPC and multi-scale mechanistic modeling provides a pathway forward to more effectively represent the heterogeneity of clinical populations and their data.