Your equation (1) describes approximately the centre-of-mass (COM) coordinates of every atom = (some nucleons + some electrons) system. Of course there are many other degrees of freedom that are not taken into account in this description. But those degrees of freedom can always be ignored unless they become correlated with the centre-of-mass coordinates. This will happen if, for example, molecules form, or if inelastic collisions result in transitions between internal electronic states. At low densities and ultracold temperatures (the normal regime for alkali atom BEC experiments), such processes are strongly suppressed. This is why the internal degrees of freedom are normally not written in the wavefunction explicitly. If you want to write them in, you would have something like
$$ \Psi^{\prime}(\mathbf{r},\mathbf{s}) = \Psi(r_1,r_2,\ldots) \times \Phi(s_1,s_2,\ldots), $$
where $\Psi(\mathbf{r})$ is given by your equation (1) and $\Phi(\mathbf{s})$ describes all the (relative) coordinates $s_1,s_2,\ldots$ of the electrons and nucleons. The lack of correlations between COM and internal states means that $\Phi(\mathbf{s})$ appears as a boring common factor multiplying all expressions, so it can just be thrown away*.
The internal degrees of freedom are nevertheless extremely important since they result in interactions between the atoms. This interaction can be modelled fairly well by a Lennard-Jones type potential. This can be understood intuitively as a sum of two parts: 1) a long-range attraction between the nucleons and electrons, 2) a short-range Coulomb/Pauli repulsion once the electron orbitals start to overlap. However, this intuition only makes sense within the Born-Oppenheimer approximation, in which the quantum correlations between nuclear and electronic degrees of freedom are neglected. In principle, this interaction is a horrendously complicated many-fermion problem. Even for just a single pair of $^{87}$Rb atoms, the full scattering problem taking all the nuclear and electronic degrees of freedom couldn't be solved on even the biggest supercomputer in the world. Fortunately, there is never enough motional energy (at low temperatures) in the centre-of-mass degrees of freedom to cause transitions between the internal states, meaning that there is no chance for the internal structure to become correlated with the external motion. It is therefore perfectly fine in practice to treat the atoms as structureless blobs of matter, with some effective classical potential replacing the complicated quantum motions and interactions of their internal parts.
The important role of the interatomic scattering interaction in cold atom experiments is that it correlates the centre-of-mass coordinates $r_1,r_2,\ldots$ of the atoms. Since your equation (1) neglects these correlations, it can only be an approximation. What you have written is called the Gross-Pitaevskii (GP) approximation for the ground state. This is the zeroth-order term in a systematic expansion of the ground state in terms of the parameter $(na_s^3)^{1/2}$, where $n$ is the density and $a_s$ is the $s$-wave scattering length. The next order of approximation is called the "Bogoliubov vacuum", and takes the form
$$\Psi(r_1,r_2,\ldots) = \prod_{i<j} \psi(r_i-r_j).$$
The reason that the GP approximation normally gives good results is that the pair wave function $\psi(r)$ only has a non-trivial behaviour for very small $r$, i.e. less than 100 nm. If the spatial resolution of your measuring device is larger than this, which is normally the case for diffraction-limited light scattering, the GP and Bogoliubov ground states give equivalent predictions. A good list of references on these issues can be found in the Leggett RMP is a good source. Leggett also wrote an interesting open-access paper on similar topics.
*Of course, if you assume that the positions of the electrons or nucleons could be measured, then $\Phi(\mathbf{s})$ would need to be taken into account.
