We consider nonlinear mixed effects models including high-dimensional covariates to model individual parameters variability. The objective is to identify relevant covariates among a large set under sparsity assumption and to estimate model parameters. To face the high dimensional setting we consider a regularized estimator namely the maximum likelihood estimator penalized with the l1-penalty. We rely on the use of the eBIC model choice criteria to select an optimal reduced model. Then we estimate the parameters by maximizing the likelihood of the reduced model. We calculate in practice the maximum likelihood estimator penalized with the l1-penalty though a weighted proximal stochastic gradient descent algorithm with an adaptive learning rate. This choice allows us to consider very general models, in particular models that do not belong to the curved exponential family. We demonstrate first in a simple linear toy model through a simulation study the good convergence properties of this optimization algorithm. We compare then the performance of the proposed methodology with those of the \glmmLasso procedure in a linear mixed effects model in a simulation study. We illustrate also its performance in a nonlinear mixed-effects logistic growth model through simulation. We finally highlight the beneficit of the proposed procedure relying on an integrated single step approach regarding two others two steps approaches for variable selection objective.