A methodology termed the ‘filtered mass density function’ (FMDF) is developed and implemented for large-eddy simulation (LES) of variable-density chemically reacting turbulent flows at low Mach numbers. This methodology is based on the extension of the ‘filtered density function’ (FDF) scheme recently proposed by Colucci et al. (1998) for LES of constant-density reacting flows. The FMDF represents the joint probability density function of the subgrid-scale (SGS) scalar quantities and is obtained by solution of its modelled transport equation. In this equation, the effect of chemical reactions appears in a closed form and the influences of SGS mixing and convection are modelled. The stochastic differential equations (SDEs) which yield statistically equivalent results to those of the FMDF transport equation are derived and are solved via a Lagrangian Monte Carlo scheme. The consistency, convergence, and accuracy of the FMDF and the Monte Carlo solution of its equivalent SDEs are assessed. In non-reacting flows, it is shown that the filtered results via the FMDF agree well with those obtained by the ‘conventional’ LES in which the finite difference solution of the transport equations of these filtered quantities is obtained. The advantage of the FMDF is demonstrated in LES of reacting shear flows with non-premixed reactants. The FMDF results are appraised by comparisons with data generated by direct numerical simulation (DNS) and with experimental measurements. In the absence of a closure for the SGS scalar correlations, the results based on the conventional LES are significantly different from those obtained by DNS. The FMDF results show a closer agreement with DNS. These results also agree favourably with laboratory data of exothermic reacting turbulent shear flows, and portray several of the features observed experimentally.