How can I design Nyquist interpolation filters with the Parks-McClellan algorithm?

We can easily design interpolation filters that obey certain frequency-domain constraints using the Parks-McClellan algorithm. However, it's not immediately clear how to enforce time-domain constraints; in particular, I'm interested in generating Nyquist filters. So if I'm oversampling by a factor of N, I want the filter to have zero-crossings at kN, for non-zero integer k (this ensures that the input samples to my interpolator will appear in the output sequence).

I've seen Harris¹ talk about a technique for designing half-band filters, i.e. the special case where N=2. Is there a general solution for this? (I know that we can easily design filters with the window method, but that doesn't give us the same control.)

_{[1] Multirate Signal Processing for Communication Systems, pp. 208-209}