We propose a method to incorporate the coupling between shape and pairing collective degrees of freedom in the framework of the interacting boson model (IBM), based on the nuclear density functional theory. To account for pairing vibrations, a boson-number non-conserving IBM Hamiltonian is introduced. The Hamiltonian is constructed by using solutions of self-consistent mean-field calculations based on a universal energy density functional and pairing force, with constraints on the axially-symmetric quadrupole and pairing intrinsic deformations. By mapping the resulting quadrupole-pairing potential energy surface onto the expectation value of the bosonic Hamiltonian in the boson condensate state, the strength parameters of the boson Hamiltonian are determined. An illustrative calculation is performed for $^{122}$Xe, and the method is further explored in a more systematic study of rare-earth $N=92$ isotones. The inclusion of the dynamical pairing degree of freedom significantly lowers the energies of bands based on excited $0^+$ states. The results are in quantitative agreement with spectroscopic data, and are consistent with those obtained using the collective Hamiltonian approach.