The yield stress and shear thinning properties of mucus are identified as critical for ciliary coordination and mucus transport in human airways. We use here numerical simulations to explore the hydrodynamic coupling of cilia and mucus with these two properties using the Herschel–Bulkley model, in a lattice Boltzmann solver for the fluid flow. Three mucus flow regimes, i.e. a poorly organized regime, a swirly regime, and a fully unidirectional regime, are observed and analysed by parametric studies. We systematically investigate the effects of ciliary density, interaction length, Bingham number and flow index on the mucus flow regime formation. The underlying mechanism of the regime formation is analysed in detail by examining the variation of two physical quantities (polarization and integral length) and the evolution of the flow velocity, viscosity and shear-rate fields. Mucus viscosity is found to be the dominant parameter influencing the regime formation when enhancing the yield stress and shear thinning properties. The present model is able to reproduce the solid body rotation observed in experiments (Loiseau et al. , Nat. Phys. , vol. 16, 2020, pp. 1158–1164). A more precise prediction can be achieved by incorporating non-Newtonian properties into the modelling of mucus as proposed by Gsell et al. ( Sci. Rep. , vol. 10, 2020, 8405).