The paper demonstrates the use of the response surface method (RSM) to carry out probabilistic assessment of the seismic bearing capacity of shallow footings seated near naturally occurring heterogeneous slopes. To this end, a pseudo-static loading is applied to a randomly uniform slope, which is homogeneous in each case but random between realizations. The method substantially reduces the number of Monte Carlo simulations required to carry out cumbersome probabilistic slope stability analyses. A finite element limit analysis model based on the lower bound theorem is developed, which is then used to generate a large synthetic database of numerical results for the seismic bearing capacity of shallow foundations resting on inherently variable natural slopes. To this end, a permutation of the key parameters is formed and lower bound FELA-based limit loads are sought through optimization in MATLAB. A closed-form solution is formulated using RSM-based polynomials. The RSM equations, which are acquired from least squares regression analyses, are used to carry out probabilistic Monte Carlo simulations and the results are presented in forms of cumulative distribution functions. Results from the probabilistic analyses are introduced into some reliability-based design approach to render design loads for different reliability levels.