We present a theoretical framework based on a higher order density correlation function, analogous to that used to investigate spin glasses, to describe dynamical heterogeneities in simulated glass-forming liquids. These higher order correlation functions are a four-point, time-dependent density correlation function g₄(r,t) and a corresponding 'structure factor' S₄(q,t) which measure the spatial correlations between the local liquid density at two points in space, each at two different times. g₄(r,t) and S₄(q,t) were extensively studied via molecular dynamics simulations of a binary Lennard-Jones mixture approaching the mode coupling temperature from above in Franz et al (1999 Phil. Mag.  B 79  1827), Donati et al (2002 J. Non-Cryst. Solids  307  215), Glotzer et al (2000 J. Chem. Phys.  112  509), Lacević et al (2002 Phys. Rev.  E 66  030101), Lacević et al (2003 J. Chem. Phys.  submitted) and Lacević (2003 Dissertation  The Johns Hopkins University). Here, we examine the contribution to g₄(r,t), S₄(q,t) and the corresponding dynamical correlation length, as well as the corresponding order parameter Q(t) and generalized susceptibility χ₄(t), from localized particles. We show that the dynamical correlation length ξ₄^SS(t) of localized particles has a maximum as a function of time t, and the value of the maximum of ξ₄^SS(t) increases steadily in the temperature range approaching the mode coupling temperature from above.