I'm assuming you're doing your interpolation something like the following way: Given a set of $(x,y,z)$ coordinates $(x_i,y_i,z_i)$ and corresponding value $f_i$ for $i=1,2,...,N$, find the coefficients $c_j$ of the multivariate polynomial
$f(x,y,z) = \sum_{j=1}^M c_j x^{\alpha_j} y^{\beta_j} z^{\gamma_j}$
such that
$f(x_i, y_i, z_i) = f_i$
Here $\alpha_j$, $\beta_j$, and $\gamma_j$ define the orders of each monomial term.
This yields a matrix equation of the following form
$A c = f$
where $A_{ij} = x_i^{\alpha_j} y_i^{\beta_j} z_i^{\gamma_j}$. This is a multivariate version of the Vandermonde matrix. If $A$ is square, then solving for $c$ gives your solution. Obviously, if there are multiple coordinates at the same location, $A$ is singular. I think there are other requirements to guarantee non-singular $A$.
If you have more unknowns ($M$) than constraints ($N$) then your system is underdetermined, which means that there are an infinite number of solutions that satisfy your constraints. It's not that you can't use this polynomial order, it's just that you can't find a unique solution.
If you have more constraints ($N$) than unknowns ($M$) then your system is overdetermined, so there is in general no solution which satisfies all of the constraints. You could solve for $c$ in the least-squares sense, by solving
$c = (A^T A)^{-1} A^T f$
If you need to guarantee that $f(x_i,y_i,z_i)=f_i$ exactly, then this won't work for you. But if you just need it to be close, then this could be the ticket.
You could also start with the underdetermined case and pick some of the highest order terms to remove until $M=N$. Then (assuming $A$ isn't singular) your solution is unique and your constraints are exactly satisfied.
Some nitpicking:
You say this is interpolation in "4D space". I suppose you're considering $f$ to be the fourth dimension, but that's confusing. It's really interpolation of a field $f$ in 3D $(x,y,z)$ space.
You say that a 6th order multivariate polynomial in 3 dimensions has 343 degrees of freedom. Technically, a $p^{th}$ other polynomial contains no terms with order higher than $p$, or in other words $\alpha_j+\beta_j+\gamma_j \le p$ for all $j$. There are ${p+d \choose d}$ degrees-of-freedom, which for $d=3$ is $\frac{(p+3)(p+2)(p+1)}{6}$, and for $p=6,d=3$ is 84.
What you're talking about is something akin to trilinear or tricubic functions, but for higher order. I don't know what you would actually call them, and I also don't really know why you would use them over the normal multivariate polynomial complete functions.