I think the subscript $\text{IN}$ refers to operators in interaction picture rather than "incoming" fields (what are those anyway?).
There exists a way to morph the second expression into some interaction-picture operators sandwiched between the free theory vacuum:
$$ \langle \Omega|T[\phi(x_1)\phi(x_2)...\phi(x_n)]|\Omega\rangle = \langle \Omega|T[\phi_{IN}(x_1)\phi_{IN}(x_2)...\phi_{IN}(x_n) \exp(-i\int H_1^{IN}(z)d^4z)]|\Omega\rangle = $$
$$ = \frac{\langle 0|T[\phi_{IN}(x_1)\phi_{IN}(x_2)...\phi_{IN}(x_n) \exp(-i\int H_1^{IN}(z)d^4z)]|0\rangle}{\langle 0|T[\exp(-i\int H_1^{IN}(z)d^4z)]|0\rangle}. $$
This justifies the use of Wick's theorem. Diagrammatically it is equivalent to the following statement: we only have to consider connected Feynman diagrams (the ones which don't contain vacuum bubbles).
By the way, we could only use the free theory vacuum state in the context of interacting theory in the interaction-picture calculations. This is a mathematical trick. Physically, there is no place for a free theory state in the interacting QFT.
Since OP asked for it in the comments, I've decided to provide a detailed derivation of the relation between the interacting vacuum and the free vacuum.
Lets consider the following exponentiation of the total Hamiltonian:
$$ e^{-i\hat{H}(T+t_{0})}=\sum_{N}e^{-iE_{N}(T+t_{0})}\left|N\right>\left<N\right|. $$
Now comes a dirty trick: for large enough $T$
all the interactions would “die out” and what remains left is the vacuum mode:
$$ \lim_{T\rightarrow\infty}e^{-i\hat{H}(T+t_{0})}\sim e^{-iE_{\Omega}(T+t_{0})}\left|\Omega\right>\left<\Omega\right|. $$
We could try to make this argument a little more mathematically precise at the cost of losing physical intuition by setting $T\rightarrow\infty(1-i\epsilon)$.
In this case only the lowest-energy mode (the vacuum state) of the Hamiltonian operator survives exponentiation since the exponential now contains a small real part.
It could also be understood as follows: in the large $T$
limit all the oscillations in the exponential become much more rapid than the vacuum mode, which gives the leading behavior. The last explanation provides some physical insights, but lacks mathematical precision.
Anyways, we use this trick to act with the mentioned above exponential on the free theory vacuum state $\,\left|0\right>$:
$$ \lim_{T\rightarrow\infty}e^{-i\hat{H}(T+t_{0})}\,\left|0\right>\sim e^{-iE_{\Omega}(T+t_{0})}\left<\Omega|0\right>\,\left|\Omega\right>.$$
On the other hand, we choose the free theory Hamiltonian $\hat{H}_{0}$
to be normal-ordered, which means that the free theory vacuum energy is zero:
$$\hat{H}_{0}\left|0\right>=0;\quad\Longrightarrow\quad e^{i\hat{H}_{0}(T+t_{0})}\,\left|0\right>=\left|0\right>.$$
Therefore,
$$e^{-i\hat{H}(T+t_{0})}\left|0\right>=e^{-i\hat{H}(T+t_{0})}e^{i\hat{H}_{0}(T+t_{0})}\left|0\right>=e^{i\hat{H}((-T)-t_{0})}e^{-i\hat{H}_{0}((-T)-t_{0})}\left|0\right>=\hat{U}^{\dagger}(-T,t_{0})\,\left|0\right>=\hat{U}(t_{0},-T)\,\left|0\right>.$$
Combining these two results, we get:
$$\left|\Omega\right>\sim\frac{\hat{U}(t_{0},-T)\left|0\right>}{e^{-iE_{\Omega}(T+t_{0})}\left<\Omega|0\right>},$$
where we implicitly assume that $T\rightarrow\infty$
. Unless the Hilbert product of the two vacuums is zero (which would mean that $\hat{H}$
can in no way be treated as a perturbation of $\hat{H}_{0}$
), the denominator is just a numerical constant.
Similarly, by acting with the Hermitian conjugate exponential (in order to still take the $T\rightarrow\infty$
limit in the end) on the bra vacuum, we obtain:
$$\left<\Omega\right|\sim\frac{\left<0\right|\,\hat{U}(T,t_{0})}{e^{-iE_{\Omega}(T-t_{0})}\left<0|\Omega\right>}.$$
The normalization condition dictates:
$$\left<\Omega|\Omega\right>=\frac{\left<0\right|\hat{U}(T,t_{0})\,\hat{U}(t_{0},-T)\left|0\right>}{e^{-iE_{\Omega}(T-t_{0})}e^{-iE_{\Omega}(T+t_{0})}\left<\Omega|0\right>\left<0|\Omega\right>}=\frac{\left<0\right|\hat{U}(T,-T)\left|0\right>}{e^{-2iE_{\Omega}T}\cdot\left|\left<0|\Omega\right>\right|^{2}}=1;$$
$$e^{-2iE_{\Omega}T}\cdot\left|\left<0|\Omega\right>\right|^{2}=\left<0\right|\hat{U}(T,-T)\left|0\right>.$$
Now we are ready to derive the expression for the correlations in the interaction picture.
$$\left<\phi_{1}(x_{1})\dots\phi_{n}(x_{n})\right>=\left<\Omega\right|\text{T}\,\hat{\phi}_{1}(x_{1})\dots\hat{\phi}_{n}(x_{n})\left|\Omega\right>=$$
$$=\frac{1}{e^{-2iE_{\Omega}T}\cdot\left|\left<0|\Omega\right>\right|^{2}}\cdot\left<0\right|\hat{U}(T,t_{0})\,\text{T}\left\{ \hat{U}(t_{0},x_{1}^{0})\,\hat{\phi}_{I\,1}(x_{1})\,\hat{U}(x_{1}^{0},t_{0})\dots\right\} \,\hat{U}(t_{0},-T)\left|0\right>.$$
We can glue together the evolution operators between interaction-picture field operators (inside the chronological ordering symbol) by using the composition law:
$$\dots\hat{\phi}_{I\,k}(x_{k})\,\hat{U}(x_{k}^{0},t_{0})\hat{U}(t_{0},x_{k+1}^{0})\,\hat{\phi}_{I\,{k+1}}(x_{k+1})\dots=\dots\hat{\phi}_{I\,k}(x_{k})\,\hat{U}(x_{k}^{0},x_{k+1}^{0})\,\hat{\phi}_{I\,{k+1}}(x_{k+1})\dots.$$
The next step would be to deal carefully with time-ordering. Each of the evolution operators is already time-ordered, because there is a time-ordered exponential in the Dyson formula. There is more: we can shuffle the operators inside the chronological ordering symbol any way we want (since they will get ordered chronologically after all). The last observation would be that since the only two evolution operators outside the chronological ordering symbol are (because of Dyson formula) also time-ordered, we could move them inside without changing anything. After all, all the evolution operators (reshuffled the way we want, inside the chronological ordering brackets) can be nicely glued together to give $\hat{U}(T,-T)$:
$$\left<\phi_{1}(x_{1})\dots\phi_{n}(x_{n})\right>=\frac{1}{N}\cdot\left<0\right|\text{T}\left\{ \hat{U}(T,-T)\,\hat{\phi}_{I1}(x_{1})\dots\hat{\phi}_{In}(x_{n})\right\} \left|0\right>,$$
where $N=e^{-2iE_{\Omega}T}\cdot\left|\left<0|\Omega\right>\right|^{2}$ is the normalization factor which we know (it was derived above) is equal to
$$N=\left<0\right|\hat{U}(T,-T)\left|0\right>=\left<0\right|\text{T}\,\hat{U}(T,-T)\left|0\right>.$$
Therefore, correlations of interacting quantum fields can be expressed formally in the interaction picture:
$$\left<\phi_{1}(x_{1})\dots\phi_{n}(x_{n})\right>=\frac{\left<0\right|\text{T}\left\{ \hat{S}\,\cdot\,\hat{\phi}_{I1}(x_{1})\dots\hat{\phi}_{In}(x_{n})\right\} \left|0\right>}{\left<0\right|\text{T}\,\hat{S}\left|0\right>},$$
where the scattering operator $\hat{S}$ is defined to be
$$\hat{S}=\lim_{T\rightarrow\infty}\,\hat{U}(T,-T)=\text{T}\,\exp\left\{ -i\intop_{-\infty}^{+\infty}dt\,\hat{H}_{I}(t)\right\} .$$
