Abstract The aim of the present work is to use one of the machine learning techniques named the genetic programming (GP) to model the p-p interactions through discovering functions. In our study, GP is used to simulate and predict the multiplicity distribution of charged pions (P(n ch)), the average multiplicity (〈n ch〉) and the total cross section (σ tot) at different values of high energies. We have obtained the multiplicity distribution as a function of the center of mass energy ($$ \sqrt s $$) and charged particles (n ch). Also, both the average multiplicity and the total cross section are obtained as a function of $$ \sqrt s $$. Our discovered functions produced by GP technique show a good match to the experimental data. The performance of the GP models was also tested at non-trained data and was found to be in good agreement with the experimental data.