Different equations were employed to directly model individual lactation curves of Tunisian Holstein-Friesian cows using 260 241 test day records for milk yield, fat, and protein contents. Eleven mathematical models were compared. Parametric curves (Legendre polynomials, Ali and Schaeffer (AS), Wood, and Wilmink models) and regression splines were tested. Goodness of fit was assessed by considering the adjusted R-square ranked according to five classes, statistical criteria, and residuals analysis. Regression splines (quadratic and cubic spline models with three knots) showed better fitting performances and greater flexibility for all milk traits. The sixth-order Legendre orthogonal polynomial model and the Ali and Schaeffer model also gave the best fit for milk traits compared with the three-parameter models (Wood and Wilmink) and with the lower-order polynomial models, but the AS model was less correlated for fat content (R = 0.87), gave a moderate correlation (R = 0.90) for protein content, and had less similarity between the observed and estimated lactation curves. The performance of Legendre orthogonal polynomials and quadratic splines was strongly affected by the models’ order and the number of knots.