Gaussian process modeling (also known as kriging) is an empirical modeling approach that has been widely applied in engineering for the approximation of deterministic functions, due to its flexibility and ability to interpolate observed data. Despite its statistical properties, Gaussian process models (GPM) have not been employed to describe the dynamics of stochastic systems with multiple outputs. Our paper presents a methodology to construct approximate models for multivariate stochastic dynamic simulations using GPM, by combining ideas from design of experiments, spatial statistics and dynamic systems modeling. The methodology is the first application in dynamic systems modeling that combines parameter and state uncertainty propagation in Gaussian process models. We apply the methodology in the prediction of a dynamic size distribution during the synthesis of nanoparticles. The method is robust to the simulation noise, and is able to learn the dynamics using a small number of sequentially designed samples of the nanoparticle simulation. (C) 2010 Elsevier Ltd. All rights reserved.