Sorry if that phrasing was too mean; I've described many models that I judged to be uselessly complex as "numerical sludge" in the past, including my own less successful efforts. Do you at least understand my specific objection above? You're almost certainly a better biologist than me, but differential equation models get a lot closer to my real expertise.
The output of their model is the structure of the resulting phylogenetic tree. I believe the primary determinant of that will be the extent of overdispersion of transmission, i.e. the shape of the pdf of how many people you expect each patient to infect. For example, let's say R0 = 2 and each patient always infects exactly two others. Then all lineages will survive (at least initially, until a significant fraction of the population is recovered and immune). The phylogenetic tree will be as perfect a binary tree as the mutation rate and our sampling permit. If we cut off the beginning, then we'll see many descendant lineages.
Instead, let's say the spread is absurdly overdetermined, so that each patient infects either zero others (with p = 0.99) or 200 (with p = 0.01). Then (after trying the simulation enough times for the pandemic to finally get started) our tree will just be a succession of super-spreading events, with far fewer lineages.
I believe that by manipulating that pdf, we can reshape that phylogenetic tree to almost anything we want. I see no description from the authors of why they believe the pdf resulting from their choice of a contact network is the right one, nor any study of robustness to changes in the shape of that pdf. Am I missing something?