We enhance the Bayesian optimization (BO) approach for simulation-based design of engineering systems consisting of multiple interconnected expensive simulation models. The goal is to find the global optimum design with minimal model evaluation costs. A commonly used approach is to treat the whole system as a single expensive model and apply an existing BO algorithm. This approach is inefficient due to the need to evaluate all the component models in each iteration. We propose a multi-model BO approach that dynamically and selectively evaluates one component model per iteration based on the uncertainty quantification of linked emulators (metamodels) and the knowledge gradient of system response as the acquisition function. Building on our basic formulation, we further solve problems with constraints and feedback couplings that often occur in real complex engineering design by penalizing the objective emulator and reformulating the original problem into a decoupled one. The superior efficiency of our approach is demonstrated through solving two analytical problems and the design optimization of a multidisciplinary electronic packaging system.