This paper studies computer simulation methods for estimating the sensitivities (gradient, Hessian etc.) of the expected steady-state performance of a queueing model with respect to the vector of parameters of the underlying distribution (an example is the gradient of the expected steady-state waiting time of a customer at a particular node in a queueing network with respect to its service rate). It is shown that such a sensitivity can be represented as the covariance between two processes, the standard output process (say the waiting time process) and what we call the score function process which is based on the score function. Simulation procedures based upon such representations are discussed, and in particular a control variate method is presented. The estimators and the score function process are then studied under heavy traffic conditions. The score function process, when properly normalized, is shown to have a heavy traffic limit involving a certain variant of two-dimensional Brownian motion for which we describe the stationary distribution. From this, heavy traffic (diffusion) approximations for the variance constants in the large sample theory can be computed and are used as a basis for comparing different simulation estimators. Finally, the theory is supported by numerical results.