To obtain this form of the Hamiltonian you just need to rewrite the sum
$$
\sum_{k,k',q}V(q) a^\dagger_{k+q}a^\dagger_{k'-q}a_{k'}a_k
$$
in the form
$$
V(0)a^\dagger_0a^\dagger_0a_0a_0 + 2\sum_{k\neq 0}(V(0)+V(k))a^\dagger_k a_k a^\dagger_0 a_0 + \sum_{k\neq 0} V(k)(a^\dagger_k a^\dagger_{-k} a_0 a_0 + a^\dagger_0 a^\dagger_0 a_{-k} a_k)+
$$
$$
+ 2\sum_{k,q\neq0}V(q)(a^\dagger_0 a^\dagger_{k+q}a_ka_q + a^\dagger_{k+q}a^\dagger_{-q}a_ka_0) + \sum_{k,k',q\neq0}V(q)a^\dagger_{k+q}a^\dagger_{k'-q}a_{k'}a_k
$$
and then to change operators $a^\dagger_0$ and $a_0$ by c-number $\sqrt{N_0}$.
By the way, I think there is a misprint in the article in a term containing only one $a_0$ or $a^\dagger_0$ operator.
Upd. I just made use of simple formulas like:
$$
\sum_{k} f_k = f_0 + \sum_{k\neq0} f_k, \qquad
\sum_{k,k'} F_{k,k'} = F_{0,0} + \sum_{k\neq0} F_{k,0} + \sum_{k'\neq0}F_{0,k'} +\sum_{k,k'\neq0}F_{k,k'}
$$
and so on. You need to select terms containing different combinations of $a^\dagger_0$ and $a_0$ operators.