Gas microbubbles stabilised by a surfactant or polymer coating are widely used in echocardiography and, increasingly, for quantitative studies of tissue perfusion. It is their highly nonlinear response to ultrasound excitation that makes them such effective contrast agents. It is this property also, however, that leads to image artefacts and difficulties in obtaining accurate quantitative information. There are few theoretical models describing the response of a contrast agent population that take into account both nonlinear behaviour and multiple bubble interactions. Computational complexity arises when the medium contains a polydisperse bubble population because a nonlinear ordinary differential equation (ODE) governing the bubble response must be computed for the current radius of each bubble size R0 present at each spatial location at every time step, which unfortunately makes the numerical model impractical for real-time clinical use. Further complexity occurs when near-resonant microbubbles shed their lipid coating drastically transforming the suspension’s acoustic response to subsequent pulses. In this work, we investigate approximations that can significantly reduce the computational complexity and demonstrate that, under certain parameter regimes, these approximations can simulate the nonlinear propagation of a repeated ultrasound pulses through a polydisperse contrast agent suspension to reasonably high accuracy.