  Halo occupation distribution (HOD) models describe the number of galaxies
that reside in different haloes, and are widely used in galaxy-halo connection
studies using the halo model (HM). Here, we introduce and study HOD response
functions $R_\mathcal{O}^g$ that describe the response of the HODs to
long-wavelength perturbations $\mathcal{O}$. The linear galaxy bias parameters
$b_\mathcal{O}^g$ are a weighted version of $b_\mathcal{O}^h +
R_\mathcal{O}^g$, where $b_\mathcal{O}^h$ is the halo bias, but the
contribution from $R_\mathcal{O}^g$ is routinely ignored in the literature. We
investigate the impact of this by measuring the $R_\mathcal{O}^g$ in separate
universe simulations of the IllustrisTNG model for three types of
perturbations: total matter perturbations, $\mathcal{O}=\delta_m$; baryon-CDM
compensated isocurvature perturbations, $\mathcal{O}=\sigma$; and potential
perturbations with local primordial non-Gaussianity, $\mathcal{O}\propto f_{\rm
NL}\phi$. Our main takeaway message is that the $R_\mathcal{O}^g$ are not
negligible in general and their size should be estimated on a case-by-case
basis. For stellar-mass selected galaxies, the responses $R_\phi^g$ and
$R_\sigma^g$ are sizeable and cannot be neglected in HM calculations of the
bias parameters $b_\phi^g$ and $b_\sigma^g$; this is relevant to constrain
inflation using galaxies. On the other hand, we do not detect a strong impact
of the HOD response $R_1^g$ on the linear galaxy bias $b_1^g$. These results
can be explained by the impact that the perturbations have on
stellar-to-total-mass relations. We also look into the impact on the bias of
the gas distribution and find similar conclusions. We show that a single extra
parameter describing the overall amplitude of $R_\mathcal{O}^g$ recovers the
measured $b_\mathcal{O}^g$ well, which indicates that $R_\mathcal{O}^g$ can be
easily added to HM/HOD studies as a new ingredient.
