This study investigates the interpretability of machine learning (ML) models applied to cumulative damage prediction during a sequence of earthquakes, emphasizing the use of techniques such as SHapley Additive exPlanations (SHAP), Partial Dependence Plots (PDPs), Local Interpretable Model-agnostic Explanations (LIME), Accumulated Local Effects (ALE), Permutation and Impurity-based technique. The research explores the cumulative damage during seismic sequences, aiming to identify critical predictors and assess their influence on the cumulative damage. Moreover, the predictors contribution in respect with the range of final damage is evaluated. Nonlinear time history analyses are applied to extract the seismic response of an eight-story Reinforced Concrete (RC) frame. The regression problem’s input variables are divided into two distinct physical classes: pre-existing damage from the initial seismic event and seismic parameters representing the intensity of the subsequent earthquake, expressed by Park and Ang damage index (DIPA) and Intensity Measures (IMs), respectively. The study offers a comprehensive review of cutting-edge ML methods, hyperparameter tuning, and ML method comparisons. A LightGBM model emerges as the most efficient, among 15 different ML methods examined, with critical predictors for final damage being the initial damage caused by the first shock and the IMs of the subsequent shock: IFVF and SIH. The importance of these predictors is supported by feature importance analysis and local/global explanation methods, enhancing the interpretability and practical utility of the developed model.