Phytoremediation, is a promising biochemical process which has gained wide acceptance in remediating the contaminants from the soil. Phytoremediation process comprises of biochemical mechanisms such as adsorption, transport, accumulation and translocation. State-of-the-art modelling methods used for studying this process in soil are limited to the traditional ones. These methods rely on the assumptions of the model structure and induce ambiguity in its predictive ability. In this context, the Artificial Intelligence approach of Genetic programming (GP) can be applied. However, its performance depends heavily on the architect (objective functions, parameter settings and complexity measures) chosen. Therefore, this present work proposes a comprehensive study comprising of the experimental and numerical one. Firstly, the lead removal efficiency (%) from the phytoremediation process based on the number of planted spinach, sampling time, root and shoot accumulation of the soil is measured. The numerical modelling procedure comprising of the two architects of GP investigates the role of the two objective functions (SRM and AIC) having two complexity measures: number of nodes and order of polynomial in modelling this process. The performance comparison analysis of the proposed models is conducted based on the three error metrics (RMSE, MAPE and R) and cross-validation. The findings reported that the models formed from GP architect using SRM objective function and order of polynomial as complexity measure performs better with lower size and higher generalization ability than those of AIC based GP models. 2-D and 3-D surface analysis on the selected GP architect suggests that the shoot accumulation influences (non-linearly) the lead removal efficiency the most followed by the number of planted spinach, the root accumulation and the sampling time. The present work will be useful for the experts to accurately determine lead removal efficiency based on the explicit GP model, thus saving the waste of input resources.