Analyses of recent helioseismic data have produced ample evidence for substantial dynamical variation of the differential rotation within the solar convection zone. Given the inevitable difficulties in resolving the precise nature of variations in deeper layers, much effort has recently gone into determining theoretically the expected modes of behaviour, using nonlinear dynamo models. Two important limitations of these models are that they have so far included only one form of nonlinearity, and as yet they have not taken into account the density stratification in the solar convection zone. Here we address both of these issues by studying the effects of including density stratification, as well as including an α-quenching nonlinearity in addition to the previously studied effects of the Lorentz force on the differential rotation. We find that observationally important features found in the earlier uniform density models remain qualitatively unchanged, although there are quantitative differences. This is important as it provides more realistic theoretical predictions to be compared with and guide observations, especially in the deeper regions where the uncertainties in the inversions are larger. However the presence of an effective alpha-quenching nonlinearity significantly reduces the amplitudes of the oscillations.