Persistent bacteremia caused by Staphylococcus aureus (SA), especially methicillin-resistant SA (MRSA), is a significant cause of morbidity and mortality. Despite susceptibility phenotypes in vitro, persistent MRSA strains fail to clear with appropriate anti-MRSA therapy during bacteremia in vivo. Thus, identifying the factors that cause such MRSA persistence is of direct and urgent clinical relevance. To address the dynamics of MRSA persistence in the face of host immunity and typical antibiotic regimens, we developed a mathematical model based on the overarching assumption that phenotypic heterogeneity is a hallmark of MRSA persistence. First, we applied an ensemble modeling approach and obtained parameter sets that satisfied the condition of a minimum inoculum dose to establish infection. Second, by simulating with the selected parameter sets under vancomycin therapy which follows clinical practices, we distinguished the models resulting in resolving or persistent bacteremia, based on the total SA exceeding a detection limit after five days of treatment. Third, to find key determinants that discriminate resolving and persistent bacteremia, we applied a machine learning approach and found that the immune clearance rate of persister cells is a key feature. But, fourth, when relapsing bacteremia was considered, the growth rate of persister cells was also found to be a key feature. Finally, we explored pharmacological strategies for persistent and relapsing bacteremia and found that a persister killer, but not a persister formation inhibitor, could provide for an effective cure the persistent bacteremia. Thus, to develop better clinical solutions for MRSA persistence and relapse, our modeling results indicate that we need to better understand the pathogen-host interactions of persister MRSAs in vivo.