**Abstract** : The quantitative reliability assessment of a thermal-hydraulic (T-H) passive safety system of a nuclear power plant can be obtained by (i) Monte Carlo (MC) sampling the uncertainties of the system model and parameters, (ii) computing, for each sample, the system response by a mechanistic T-H code and (iii) comparing the system response with pre-established safety thresholds, which define the success or failure of the safety function. The computational effort involved can be prohibitive because of the large number of (typically long) T-H code simulations that must be performed (one for each sample) for the statistical estimation of the probability of success or failure. In this work, Line Sampling (LS) is adopted for efficient MC sampling. In the LS method, an "important direction" pointing towards the failure domain of interest is determined and a number of conditional one-dimensional problems are solved along such direction; this allows for a significant reduction of the variance of the failure probability estimator, with respect, for example, to standard random sampling. Two issues are still open with respect to LS: first, the method relies on the determination of the "important direction", which requires additional runs of the T-H code; second, although the method has been shown to improve the computational efficiency by reducing the variance of the failure probability estimator, no evidence has been given yet that accurate and precise failure probability estimates can be obtained with a number of samples reduced to below a few hundreds, which may be required in case of long-running models. The work presented in this paper addresses the first issue by (i) quantitatively comparing the efficiency of the methods proposed in the literature to determine the LS important direction; (ii) employing artificial neural network (ANN) regression models as fast-running surrogates of the original, long-running T-H code to reduce the computational cost associated to the determination of the LS "important direction" and (iii) proposing a new technique for identifying the LS "important direction", based on the genetic algorithm (GA) minimization of the variance of the LS failure probability estimator. In addition, this work addresses the second issue by assessing the performance of the LS method in estimating small failure probabilities with a reduced (e.g., lower than one hundred) number of samples. The issues are investigated within two case studies: the first one deals with the estimation of the failure probability of a nonlinear structural system subject to creep and fatigue damages [1] and [2]; the second one regards a passive decay heat removal system in a gas-cooled fast reactor (GFR) of literature, [3].