The model is y_i = x_i'*beta + n_i'*f+epsilon_i, where beta is the parametric coefficient regression of vaiables X1,..., Xp; f is the functional form of the nonparametric curve (an spline); and epsilon_i is the error term.
If f contains no estimated parameters (as usual for splines), it is a special case of (multiple) linear regression. If your epsilon_i is Gaussian, you can use usual results to build this bands. Like predict(md, int="conf") in R, for instance.