Performance-based earthquake engineering (PBEE) has emerged as a powerful method of analysis and design philosophy in earthquake engineering and is leading the way to a new generation of seismic design guidelines. PBEE requires a comprehensive understanding of the earthquake response of Soil-Foundation-Structure-Interaction (SFSI) systems when damage occurs in the structural system during the earthquake. In the context of PBEE, this research combines finite element (FE) modeling and seismic response analysis of SFSI systems with state-of-the-art methods in response sensitivity and reliability analysis. New analytical and numerical methods are developed and existing algorithms adopted for studying the propagation of uncertainties in nonlinear static and dynamic analyses of SFSI systems and for probabilistic performance assessment of these systems. This research makes several contributions to reliability analysis of structural and SFSI systems. For the purpose of accurately and efficiently computing the response gradients, an 'exact' FE response sensitivity computation algorithm based on the Direct Differentiation Method (DDM) and available in the widely used FE analysis software framework OpenSees is further extended to various types of material models, finite elements and multi-point constraint equations used in modeling large-scale realistic SFSI systems. As a main contribution to this research, this sensitivity algorithm is extended to a multi-yield surface J₂ plasticity model used extensively to model clay soil materials in seismic response analysis. Related to response sensitivity analysis of SFSI systems, several issues are studied, such as discontinuities in response sensitivities, the relative importance of various soil and structural material parameters in regards to a specified aspect of the system response (i.e., response parameters). As contributions to the reliability analysis of structural and SFSI systems, several existing solution tools such as first-order reliability method (FORM), second-order reliability method (SORM), and various sampling techniques, such as importance sampling (IS) and orthogonal plane sampling (OPS), are implemented in OpenSees and/or further improved to solving reliability analysis problems of structural and SFSI systems. A powerful general-purpose optimization toolbox SNOPT, developed by Professor Philip Gill at UCSD, is integrated into the reliability analysis framework in OpenSees and customized for efficiently finding the design point(s) of structural and SFSI systems. For time variant reliability analysis, an existing mean upcrossing rate analysis algorithm is implemented in OpenSees and improved. It is found that the FORM approximation for mean upcrossing rate is significantly inaccurate, especially in cases of highly nonlinear response behavior of the system analyzed. In such case, the OPS method based on the design point(s) of the reliability problem significantly improves the FORM approximation of the mean upcrossing rate and therefore of the upper bound of the failure probability. In order to study the topology of limit-state surfaces (LSS) for reliability problems, a new visualization method called Multi-dimensional Visualization in Principal Plane (MVPP) is developed and implemented in OpenSees. The geometrical insight gained from the MVPP has led to the development of a new hybrid computational reliability method, called the DP-RS-Sim method, which combines the design point (DP) search, the response surface methodology (RS), and simulation techniques (Sim). This method is applied for the time invariant reliability analysis of a realistic nonlinear structural system. Several other closely related topics are studied. A simplified probabilistic response analysis method is developed taking advantage of DDM-based response sensitivity analysis. This method is then applied to a nonlinear structural and SFSI system. It is much more efficient than the crude Monte Carlo Simulation method and provides, at low computational cost, good estimates of the mean and standard deviation of the response for low to moderate level of material nonlinearity in the response. A general-purpose OpenSees- SNOPT based optimization framework was developed and applied to soil model updating problems using numerically simulated data. It is found that the optimization process is significantly more efficient when using the DDM-based over the FDM-based sensitivities. Additionally, nonlinear FE model updating is performed for an actual site, the Lotung downhole array in Taiwan, and based on data recorded during a 1986 earthquake