Protease inhibitors (PIs) target the protease (PR) enzyme to suppress viral replication. Their efficacy in human immunodeficiency virus treatment is compromised by the emergence of drug-resistant strains. Therefore, forecasting drug-resistance during viral evolution would help in the design of effective treatment strategies. We develop a probabilistic large-deviation model to infer epistatic interactions in genotypes observed in different treatment regimens and compute transition probabilities of point-mutations conditioned on the genotype and the treatment regimen. We simulate stochastic evolutionary paths weighted by such transition probabilities and show that low probable mutations are required for the viral population to evolve to diverse fit genotypes. We train classification models using a clinical data set of in vitro susceptibility tests to learn to infer the drug resistance of a genotype. We infer drug resistance along simulated evolutionary paths and predict that the combination PI-therapy of Atazanavir (ATV) and Ritonavir (RTV) is the least drug resistant. Without prior knowledge of PI-associated mutations, our model predicts known primary and secondary PI-resistant mutations as critical to drug resistance. This validates that our model learned mechanistic relations in the small data sets, tackling the challenge of sparse sequence data compared to the large combinatorial complexity of protein evolution and changing functionality in dynamic environments.