In a complex system it is desirable to reduce the number of expensive function evaluations required for an accurate estimation of the probability of failure. An efficient reliability estimation method is presented for engineering systems with multiple failure regions and potentially multiple most probable points. The method can handle implicit nonlinear limit state functions with correlated or noncorrelated random variables, which can be described by any probabilistic distribution. It uses a combination of approximate or “accurate-on-demand,” global and local metamodels, which serve as indicators to determine the failure and safe regions. Sample points close to limit states define transition regions between safe and failure domains. A clustering technique identifies all transition regions, which can be, in general, disjoint, and local metamodels of the actual limit states are generated for each transition region. Importance sampling generates sample points only in the identified transition and failure regions, thus, allowing the method to focus on the areas near the failure region and not expend computational effort on the sample points in the safe domain. A robust maximin “space-filling” sampling technique is used to construct the metamodels. Two numerical examples highlight the accuracy and efficiency of the method.