We present a general formalism to treat slowly rotating general relativistic superfluid neutron stars. As a first approximation, their matter content can be described in terms of a two-fluid model, where one fluid is the neutron superfluid, which is believed to exist in the core and inner crust of mature neutron stars, and the other fluid represents a conglomerate of all other constituents (crust nuclei, protons, electrons, etc.). We obtain a system of equations, good to second-order in the rotational velocities, that determines the metric and the matter variables, irrespective of the equation of state for the two fluids. In particular, allowance is made for the so-called entrainment effect, whereby the momentum of one constituent (e.g. the neutrons) carries along part of the mass of the other constituent. As an illustration of the developed framework, we consider a simplified equation of state for which the two fluids are described by different polytropes. We determine numerically the effects of the two fluids on the rotational frame-dragging, the induced changes in the neutron and proton densities and the inertial mass, as well as the change in shape of the star. We further discuss issues regarding conservation of the two baryon numbers, the mass-shedding (Kepler) limit and chemical equilibrium.