Mathematical and computational modeling frameworks play the leading role in the analysis and prediction of the dynamics of gene regulatory networks. The literature abounds in various approaches, all of which characterized by strengths and weaknesses. Among the others, Ordinary Differential Equations (ODE) models give a more general and detailed description of the network structure. But, analytical computations might be prohibitive as soon as the network dimension increases, and numerical simulation could be nontrivial, time-consuming and very often impracticable as precise and quantitative information on model parameters are frequently unknown and hard to estimate from experimental data. Last but not least, they do not account for the intrinsic stochasticity of regulation. In the present paper we consider a class of ODE models with stochastic parameters. The proposed method separates the deterministic aspects from the stochastic ones. Under specific assumptions and conditions, all possible trajectories of an ODE model, where incomplete knowledge of parameter values is symbolically and qualitatively expressed by initial inequalities only, are simulated in a single run from an initial state. Then, the occurrence probability of each trajectory, characterized by a sequence of parameter inequalities, is computed by associating probability density functions with network parameters. As demonstrated by its application to the gene repressilator system, the method seems particularly promising in the design of synthetic networks.