Using an approach based upon a set of auxiliary many-electron wavefunctions we present a rigorous derivation of spin-orbit coupling (SOC) within the framework of linear-response time-dependent density functional theory (LR-TDDFT). Our method is based on a perturbative correction of the non-relativistic collinear TDDFT equations using a Breit-Pauli spin-orbit Hamiltonian. The derivation, which is performed within both the Casida and Sternheimer formulations of LR-TDDFT, is valid for any basis set. The requirement of spin noncollinearity for the treatment of spin-flip transitions is also discussed and a possible alternative solution for the description of these transitions in the collinear case is also proposed. Our results are validated by computing the SOC matrix elements between singlet and triplet states of two molecules, formaldehyde and acetone. In both cases, we find excellent agreement with benchmark calculations performed with a high level correlated wavefunction method.