Cervical spine behavior for generalized loading is often characterized using a full three-dimensional flexibility matrix. While experimental studies have been aimed at determining cervical motion segment behavior, their accuracy and utility have been limited both experimentally and analytically. For example, the nondiagonal terms, describing coupled motions, of the matrices have often been omitted. Flexibility terms have been primarily represented as constants despite the known nonlinear stiffening response of the spine. Moreover, there is presently no study validating the flexibility approach for predicting vertebral motions;​ nor have the effects of approximations and simplifications to the matrix representations been quantified. Yet, the flexibility matrix currently forms the basis for all multibody dynamics models of cervical spine motion. Therefore, the purpose of this study is to fully quantify the flexibility relationships for cervical motion segments, examine the diagonal and nondiagonal components of the flexibility matrix, and determine the extent to which multivariable relationships improve cervical spine motion prediction. To that end, using unembalmed human cervical spine motion segments, a full battery of flexibility tests were performed for a neutral orientation and also following an axial pretorque. Primary and coupled matrix components were described using linear and piecewise nonlinear incremental constants. An additional approach utilized multivariable incremental relationships to describe matrix terms. Measured motions were predicted using structural flexibility methods and were evaluated using RMS error of the difference between the predicted and measured responses. Results of this study provide a full set of flexibility relationships describing primary and coupled motions for C3-C4 and C5-C6 motion segment levels. Analysis of these data indicates that a flexibility matrix using incremental responses describing primary and coupled motions offers improved predictions over using linear methods (p<0.01). However, there is no significant improvement using more generalized nonlinear terms represented by the multivariable functional approach (p<0.2). Based on these findings, it is suggested that a multivariable approach for flexibility is more demanding experimentally and analytically while not offering improved motion prediction.