We present a quantum-classical methodology for propagating the density matrix of a system coupled to a polyatomic (large molecular or solvent) environment. The system is treated via a full path integral, while the dynamics of the environment is approximated in terms of classical trajectories. We obtain quantum-classical path integral (QCPI) expressions in which the trajectories can undergo transitions to other quantum states at regular time intervals, but the cumulative probability of these transitions is governed by the local strength of the state-to-state coupling as well as the magnitude of the solvent reorganization energy. If quantum effects in the coordinates of the environment are relatively weak, an inexpensive random hop approximation leads to accurate descriptions of the dynamics. We describe a systematic iterative scheme for including quantum mechanical corrections for the solvent by gradually accounting for nonlocal "quantum memory" effects. As the length of the included memory approaches the decoherence time of the environment, the iterative QCPI procedure converges to the full QCPI result. The methodology is illustrated with application to dissipative symmetric and asymmetric two-level systems.