Models of galaxy and halo clustering commonly assume that the tracers can be treated as a continuous field locally biased with respect to the underlying mass distribution. In the peak model pioneered by Bardeen et al. [Astrophys. J. 304, 15 (1986)], one considers instead density maxima of the initial, Gaussian mass density field as an approximation to the formation site of virialized objects. In this paper, the peak model is extended in two ways to improve its predictive accuracy. First, we derive the two-point correlation function of initial density peaks up to second order and demonstrate that a peak-background split approach can be applied to obtain the k-independent and k-dependent peak bias factors at all orders. Second, we explore the gravitational evolution of the peak correlation function within the Zel’dovich approximation. We show that the local (Lagrangian) bias approach emerges as a special case of the peak model, in which all bias parameters are scale independent and there is no statistical velocity bias. We apply our formulas to study how the Lagrangian peak biasing, the diffusion due to large scale flows, and the mode coupling due to nonlocal interactions affect the scale dependence of bias from small separations up to the baryon acoustic oscillation (BAO) scale. For 2σ density peaks collapsing at z = 0.3, our model predicts a ~5% residual scale-dependent bias around the acoustic scale that arises mostly from first order Lagrangian peak biasing (as opposed to second order gravity mode coupling). We also search for a scale dependence of bias in the large scale autocorrelation of massive halos extracted from a very large N-body simulation provided by the MICE Collaboration. For halos with mass M ≳ 1014M⊙/h, our measurements demonstrate a scale-dependent bias across the BAO feature which is very well reproduced by a prediction based on the peak model.