The Lagrangian spectral relaxation ~LSR! model is extended to treat turbulent mixing of two passive scalars (fa and fb) with different molecular diffusivity coefficients ~i.e., differential-diffusion effects!. Because of the multiscale description employed in the LSR model, the scale dependence of differential-diffusion effects is described explicitly, including the generation of scalar decorrelation at small scales and its backscatter to large scales. The model is validated against DNS data for differential diffusion of Gaussian scalars in forced, isotropic turbulence at four values of the turbulence Reynolds number (Rl538, 90, 160, and 230! with and without uniform mean scalar gradients. The explicit Reynolds and Schmidt number dependencies of the model parameters allows for the determination of the Re ~integral-scale Reynolds number! and Sc ~Schmidt number! scaling of the scalar difference z5fa2fb . For example, its variance is shown to scale like ^z2& ;Re20.3. The rate of backscatter (bD) from the diffusive scales towards the large scales is found to be the key parameter in the model. In particular, it is shown that bD must be an increasing function of the Schmidt number for Sc<1 in order to predict the correct scalar-to-mechanical time-scale ratios, and the correct long-time scalar decorrelation rate in the absence of uniform mean scalar gradients.