We calculate four-loop QCD corrections to the electroweak $\rho$ parameter with a non-vanishing $b$ quark mass. At three loops, it was observed that elliptic integrals contribute to this observable. This prompts the question of which classes of functions appear at the next order. We report on the status of our calculation with a focus on the mathematical structures that emerge at four loops.