The matter power spectrum P(k) is one of the main quantities connecting observational and theoretical cosmology. Although for a fixed redshift this can be numerically computed very efficiently by Boltzmann solvers, an analytical description is always desirable. However, accurate fitting functions for P(k) are only available for the concordance model. Taking into account that forthcoming surveys will further constrain the parameter space of cosmological models, it is also of interest to have analytical formulations for P(k) when alternative models are considered. Here, we use the genetic algorithms, a machine learning technique, to find a parametric function for P(k) considering several possible effects imprinted by modifications of gravity. Our expression for the P(k) of modified gravity shows a mean accuracy of around 1-2% when compared with numerical data obtained via modified versions of the Boltzmann solver CLASS, and thus it represents a competitive formulation given the target accuracy of forthcoming surveys.