A two-dimensional (2-D) seasonal model has been developed for simulating the transient response of the climate system to the astronomical forcing. The atmosphere is represented by a zonally averaged quasi-geostrophic model which includes accurate treatment of radiative transfer. The atmospheric model interacts with the other components of the climate system (ocean, sea-ice and land surface covered or not by snow and ice) through vertical fluxes of momentum, heat and humidity. The model explicitly incorporates surface energy balances and has snow and sea-ice mass budgets. The vertical profile of the upper-ocean temperature is computed by an interactive mixed-layer model which takes into account the meridional turbulent diffusion of heat. This model is asynchronously coupled to a model which simulates the dynamics of the Greenland, the northern American and the Eurasian ice sheets. Over the last glacial–interglacial cycle, the coupled model simulates climatic changes which are in general agreement with the low frequency part of deep-sea, ice and sea-level records. However, after 6000 yBP, the remaining ice volume of the Greenland and northern American ice sheets is overestimated in the simulation. The simulated climate is sensitive to the initial size of the Greenland ice sheet, to the ice-albedo positive feedback, to the precipitation-altitude negative feedback over the ice sheets, to the albedo of the aging snow and to the insolation increase, particularly at the southern edge of the ice sheets, which is important for their collapse or surge.