This paper proposes a gradient domain deformation for wrapping surface models of muscles around bones as they move during a simulation of physiological activities. Each muscle is associated with one or more poly-lines that represent the muscle skeleton to which the surface model of the muscle is bound so that transformation of the skeleton (caused by the movement of bones) produces transformation of the vertices of the mesh subject to Laplacian linear constraints to preserve the local shape of the mesh and non-linear volume constraints to preserve the volume of the mesh. All these constraints form a system of equations that is solved using the iterative Gauss-Newton method with Lagrange multipliers. Our C++ implementation can wrap a muscle of medium size in about a couple of ms up to 400 ms on commodity hardware depending on the type of parallelization, whilst it can keep the change in volume below 0.04%. A preliminary biomechanical assessment of the proposed technique suggests that it can produce realistic results and thanks to its rapid processing speed, it might be an attractive alternative to the methods that are used in clinical practise at present.