We consider interaction splines which model a multivariate regression function $f$ as a constant plus the sum of functions of one variable (main effects), plus the sum of functions of two variables (two-factor interactions), and so on. The estimation of $f$ by the penalized least squares method and the asymptotic properties of the models are studied in this article. It is shown that, under some regularity conditions on the data points, the expected squared error averaged over the data points converges to zero at a rate of $O(N^{-2m/(2m + 1)})$ as the sample size $N \rightarrow \infty$ if the smoothing parameters are appropriately chosen, where $m$ is a measure of the assumed smoothness of $f.$