In this paper, we present constrained spline dynamics (CSD) as a unified framework for the elastodynamic simulation of elastic rods subjected to constraints at interactive rates. The geometry of the rod and its kinematics are discretized using smooth spline functions and the rod's centerline co-ordinates as degrees of freedom (DOF). Interpolating B-spline shape functions are used to take advantage of the smooth basis and the Kronecker delta property. The formulation is developed from Hamilton's principle with bending and twisting energies represented as compliant constraints. The bend-twist coupled behavior is modeled using the concept of holonomy of curves utilizing the smooth and accurate curvature and bi-normal vector fields, eliminating rotational director frames as degrees of freedom. By enforcing uniform arc-length parametrization, high accuracy is achieved in modeling bend, twist, and bend-twist coupling. Several numerical examples are presented that demonstrate the convergence behavior, computational accuracy and efficiency of the formulation.