Abstract. This paper presents an approximation method for performing reliability analysis with high fidelity computer codes at a reasonable computational cost. Complex models are common in science and engineering due to their ability to substitute costly and some times infeasible practical experiments. These models, however, suffer from high computational cost. This causes problems when performing sampling -based reliability analysis, since the failure modes of the system typically occupy a small region of the input space and thus relatively large sample sizes are required for the accurate estimation of their characteristics. The sequential sampling method proposed in this article, combines Gaussian process-based optimisation and Subset Simulation. Gaussian process emulators construct a statistical approximation to the output of the original code, which is both affordable to use and has its own measure of predictive uncertainty. Subset Simulation is used to efficiently populate those regions of the initial approximation which are likely to lead to the performance function exceeding a predefined critical threshold. Among all samples, the ones that are likely to contribute most to increasing the quality of the surrogate in the vicinity of the failure regions are selected, using Bayesian optimisation methods. The iterative nature of the method ensures that an arbitrarily accurate approximation of the failure region in performance space is developed at a reasonable computational cost. The presented method is applied to a number of benchmark problems.