We derive equations for the time evolution of the reduced density matrix of a collection of heavy quarks and antiquarks immersed in a quark gluon plasma. These equations, in their original form, rely on two approximations: the weak coupling between the heavy quarks and the plasma, the fast response of the plasma to the perturbation caused by the heavy quarks. An additional semi-classical approximation is performed. This allows us to recover results previously obtained for the abelian plasma using the influence functional formalism. In the case of QCD, specific features of the color dynamics make the implementation of the semi-classical approximation more involved. We explore two approximate strategies to solve numerically the resulting equations in the case of a quark-antiquark pair. One involves Langevin equations with additional random color forces, the other treats the transition between the singlet and octet color configurations as collisions in a Boltzmann equation which can be solved with Monte Carlo techniques.