The modelling of proton-proton collisions at high energies using both group method data handling (GMDH) and gene expression programming (GEP) approaches is investigated over a wide range of center-of-mass energy (from \( \sqrt{s} = 23.6\) GeV to 7 TeV). We have used GMDH and GEP models to obtain two different mathematical formulae expressing the complex relation of the charged-particle pseudorapidity distribution \((\frac{d N_{ch}}{d \eta})\) of proton-proton interactions as a function of the center-of-mass energy \((\sqrt{s})\) in a simple explicit form. We have used the obtained formulae to simulate and predict \( \frac{d N_{ch}}{d \eta}\). Predictions are made at \( \sqrt{s}= 10\) TeV and 14 TeV. Compared to the available experimental data and to some widely used Monte Carlo generators (PYTHIA, PHOJET, and QGSM models), it is found that our models can accurately simulate and predict the charged-particle pseudorapidity distribution of proton-proton collisions.