Introduction: Determining the status of axillary lymph nodes is extremely important for prognosis and correct therapy selection in breast cancer patients. In patients designated for neoadjuvant systemic therapy (NST), clinical examination and radiological methods are the primary means of determining the status of axillary lymph nodes. Such methods lead to as many as 30% of misdiagnosed cases, which indicates the need for supplementary methods. Due to the increasing availability of data in oncology, machine learning methods are an ideal candidate for improving the diagnosis of axillary lymph nodes in breast cancer. The goal of this dissertation is the training and validation of machine learning models (Random Forest, XGBoost, TabNet Linear Regression and KNN) on clinico-pathological data, with the purpose of accurate axillary lymph node classification. Subjects and Methods: Anonymized data from a total of 13580 patients treated for breast cancer were included in this research. The sample consists of subjects treated for breast cancer in 25 Croatian hospitals in the period from January 2017, to January 2022. Data were collected from hospital information systems during the mentioned period, by searching all patients with ICD code 50 (for breast cancer). All subjects with complete data (PHD, lymph node status, ER, PR, HER2, Ki67, tumor size, age) were included. Subjects with incomplete data (missing some of the input variables) was excluded, after which 9705 subjects remained. Among these, 1324 of them initially received NST, which is why they were also excluded from the analysis. Finally, after applying the criteria for NST, 719 patients were identified for final analysis, optimization and validation of the predictive models. Also, models were trained and validated on the entire population with complete data (n= 8381). Hyperparameter optimization was performed using the 5-fold cross validation technique. Validation of the model was performed using the "bootstrap" technique on the "test" data set. The feature importance analysis and the analysis of the model’s decision-making process was performed using the SHAP technique. Results: In the group meeting the NST criteria, the best performance achieved by the Random forest model with a mean AUC of 0.793 (95% CI: 0.713-0.865), while XGBoost had the best performance on the total population with a mean AUC of 0.762 (95% CI: 0.726-0.794). Shapley value analysis showed that the most important predictor is tumor size, followed by Ki67 index, patient age and tumor grade. The association between axillary lymph node metastases and tumor size grows linearly up to a size of 5 cm, after which a plateau is reached. The association of Ki67 with LN metastasis increases up to Ki67 of 75%, after which a decline is observed, but mostly for the triple-negative type of breast cancer. Random forest and XGBoost achieve significant performance gains over the baseline univariate logistic regression model (trained only on tumor size). Conclusion: Explainable tree-based machine learning methods trained on patient and tumor features obtained during regular pre-operative/pre-NST procedures achieve a good performance in predicting breast cancer axillary lymph node metastasis. Such models can lead to more accurate diagnosis and better treatment selection, especially for NST patients, where radiological and clinical findings are often the only way of axillary lymph node assessment. Potential upstage of diagnosis based on machine learning models for some patients would result in NST and, consecutively, potentially more adjuvant therapy with noncross resistant treatments and better patient outcomes. The addition of genetic and biomarker data and subsequent validation in multinational/multicenter studies is expected from future studies.