The friction variable is treated as another dynamical variable and integrated using symplectic integrators like the Verlet methods. There are many discussions of this, one I found via Google is here: www2.ph.ed.ac.uk/~dmarendu/MVP/MVP03.pdf

The Verlet methods provide a second-order symplectic method. Higher orders are more accurate but come at the expense of greater computational complexity, i.e., they cost more per time step but are more stable for somewhat longer time steps. Most codes live with simpler integrators that can be readily parallelized and live with smaller time steps for stability.