The aim of this study was to design physics-preserving and precise surrogate models of the non-linear elastic behaviour of an intervertebral disc (IVD). Based on artificial force-displacement data sets from detailed finite element (FE) disc models, weused greedy kernel and polynomial approximations of second, third and fourth order to train surrogate models for the scalar force-torque-potential. Doing so, the resulting models of the elastic IVD responses ensured the conservation of mechanical energy through their structure. At the sametime, they were capable of predicting disc forces for the full physiological range of motion andfor the coupling of all six degrees of freedom of an intervertebral joint. The performance of allsurrogate models for a subject-specific L4|5 disc geometry was evaluated both on training and test data obtained from uncoupled (one-dimensional), weakly coupled (two-dimensional),and random movement trajectories in the entire six-dimensional (6d) physiological displacement range, as well as on synthetic kinematic data. We observed highest precisions for the kernel surrogate followed by the fourth order polynomial model. Both clearly outperformed the second order polynomial model which is equivalent to the commonly used stiffness matrix in neuro-musculoskeletal simulations.Hence, the proposed model architectures have the potential to improve the accuracy, and, therewith, validity of load predictions in neuro-musculoskeletal spine models.