Multidisciplinary analysis (MDA) is nowadays a powerful tool for analysis and optimization of complex systems. The present study is interested in the case where MDA involves feedback loops between disciplines (i.e., the output of a discipline is the input of another and vice versa). When the models for each discipline involve non-negligible modeling uncertainties, it is important to be able to efficiently propagate these uncertainties to the outputs of the MDA. The present study introduces a polynomial chaos expansion (PCE)-based approach to propagate modeling uncertainties in MDA. It is assumed that the response of each disciplinary solver is affected by an uncertainty modeled by a random field over the design and coupling variables space. A semi-intrusive PCE formulation of the problem is proposed to solve the corresponding nonlinear stochastic system. Application of the proposed method emphasizes an important particular case in which each disciplinary solver is replaced by a surrogate model (e.g., kriging). Three application problems are treated, which show that the proposed approach can approximate arbitrary (non-Gaussian) distributions very well at significantly reduced computational cost.