I imagine the authors keep filter order as a design parameter, since dynamically changing order is a rather tricky proposition.

For IIRs the approach (imo) would be to use transformed analog biquads (such as an SVF topology, via the TPT/Zavalishin's method) which is designed to handle the time variant issues with state evolution. From there your model wouldn't synthesize filter coefficients directly but either the design parameters of the filter or the gains of the outputs or multi mode filters built from cascaded biquad sections. A good example of this in practice is the Oberheim "pole mixing" filter (analog) or the Shruthi's variation on it, which is a 4th order filter that achieves a wide array of frequency responses by mixing the outputs of each stage. Those gains in the mixer can be varied for very cool results.