Realistic simulations of local control of cardiac EC coupling require Ca release unit (CaRU) models generated using Markov chain models of L-type Ca channels and ryanodine receptors (RyRs) as “building blocks.” Because compositionally defined CaRU models result in a combinatorial explosion of release site states, most whole cell simulations to date have utilized ad hoc CaRU models in an effort to maintain computational efficiency (e.g., modeling RyR clusters as a “megachannel”). To overcome this state-space explosion, we have implemented, validated, and benchmarked several methods for automated reduction of mechanistic CaRU models that feature an automated process of state aggregation and evaluation of reduction error through comparison of the jump probability matrices of full and reduced models. When there is separation of time scales in the single channel model (e.g., fast activation and slow inactivation), we perform numerical fast/slow reduction by categorizing rate constants in the single channel model as either fast or slow, aggregating states in the expanded CaRU model that are connected by fast transitions, and calculating transition rates between lumped states using the conditional probability distribution of states within each group. For large problems where the conditional distributions can not be directly calculated from the full model, we employ iterative aggregation/disaggregation to calculate conditional distributions in a memory-efficient fashion. For problems without time scale separation, how states should be aggregated to yield good reductions can not be determined a priori. Consequently, we implemented a genetic algorithm that evolves potential schemes for state aggregation, ultimately yielding simplified CaRU models with low reduction error. We demonstrate that such automated CaRU reduction procedures can be used to accelerate multiscale models of local control of CICR in cardiac myocytes.