Importance sampling (IS) is recognized as a potentially powerful method for reducing simulation run times when estimating the probabilities of rare events in communication systems using Monte Carlo simulation. Of special interest is the probability of buffer overflow in networks of queues. When simulating networks of queues, regenerative techniques make the application of IS feasible and efficient. The application of regenerative techniques is also crucial in obtaining correct confidence intervals for the estimates involved. However, using the most favorable IS settings very often makes the length of regeneration cycles infinite or impractically long. We discuss here a methodology that uses IS dynamically within each regeneration cycle, in order to drive the system back to the regeneration state, after an accurate estimate has been obtained. To obtain large speed-up factors in simulation run time using IS, the modification or bias of the underlying probability measures must be carefully chosen. Analytically or numerically minimizing the variance of the IS estimator with respect to the biasing parameters or finding the optimal exponential change of measure is only possible under certain conditions. In this paper, we formulate a statistically bused technique for optimizing IS parameter values for simulations of queueing systems, including complex systems with bursty arrival processes. We use a deterministic variant of stochastic simulated annealing (SA) called mean field annealing (MFA) to minimize statistical estimates of the IS estimator variance. We demonstrate the technique by evaluating blocking probabilities for the M/M/1/K, M/D/1/K, GVD/1/K, Geo/Geo/1/K, and IBP/Geo/1/K queues (where Geo stands for Geometric and IBP for Interrupted Bernoulli Process), a 16 x 16 synchronous Clos Asynchronous Transfer Mode (ATM) switch, and a 4 x 4 ATM switch with priority and push-out. Run time speed-up factors of two to eleven orders of magnitude over conventional Monte Carlo are obtained for these examples.