Proposition: For any given $T>0$, $y(t,\epsilon)\rightarrow x(t)$ in probability as $\epsilon\rightarrow0$ uniformly for $t\in[0,T]$
Proof: $(\text{Eq}.(3)-\text{Eq}.(1))/\epsilon$ gives $$dy=-(k_0y+k_1(x_0-1)+\epsilon k_1y)\,dt+(\eta_0y+\eta_1x_0+\epsilon \eta_1y)\,dB. \tag5$$ $\text{Eq}.(5)-\text{Eq}.(4)$ gives $$dz = -(k_0z+\epsilon k_1y)\,dt+(\eta_0z+\epsilon \eta_1y)\,dB \tag6$$ where $z=y-x_1$.
Here I can use Duhamel's principle or the linearity of Eq.(4) to obtain $z$ as an integral of $\epsilon y$. And the convergence is pathwise. However, I would like to find a technique that is generalizable to any SDE with factors that are Lipschitz continuous $C^1$ functions. So I will proceed with the following approach.
I am going to prove for any given $T>0$, \begin{align} \mathbf E\big[y(t,\epsilon)^2\big]&<M, \tag{7.1} \\ \mathbf E\big[z(t,\epsilon)^2\big]&<\epsilon^2Be^{AT}, \tag{7.2}\\ \forall\epsilon\in(0,\epsilon_0),&\ t\in[0,T] \end{align} for some positive constants $\epsilon_0,\,M,\,A,\,B$. Therefore by the Chebyshev-Markov inequality, we have $z(t,\omega,\epsilon)\rightarrow0$ in probability as $\epsilon\rightarrow0$ uniformly for $t\in[0,T]$.
We will show the derivation of Eq.(7.2) given Eq.(7.1). A similar technique applies to Eq.(7.1) without premising on Eq.(7.2).
Take the integral form of Eq. (6), square it and apply the Cauchy-Schwartz inequality, $$\frac14z(t)^2\le\Big(\int_0^t k_0z\,ds\Big)^2+\Big(\int_0^t\epsilon k_1y\,ds\Big)^2+\Big(\int_0^t\eta_0z\,dB_s\Big)^2+\Big(\int_0^t\epsilon \eta_1y\,dB_s\Big)^2. \tag8 $$
We make use of the following inequality. For a deterministic function $a(t)$ and a stochastic function $u(t,\omega)$ where $\omega$ is an element of the sample space, \begin{align} &\mathbf E\Big[\Big(\int_0^t a(s)u(s,\omega)dB(s,\omega)\Big)^2\Big] \\ =&\mathbf E\Big[\int_0^t (a(s)u(s,\omega))^2\,ds\Big]\quad \text{(Ito isometry)} \\ \le&\int_0^t a(s)^2ds\int_0^t \mathbf E[u(s,\omega)^2]ds.\quad\text{(Cauchy-Schwartz inequality)}\tag9 \end{align}
Take expectation on Eq.(8) and apply Eq.(9) and again Cauchy-Schwartz, we have $$\mathbf E[z(t)^2]\le \alpha(t)\int_0^t\mathbf E[z(s)^2]\,ds +\epsilon^2\beta(t)\int_0^t\mathbf E[y(s)^2]\,ds$$ for some positive nondecreasing continuous deterministic funtions $\alpha(t)$ and $\beta(t)$. Let $v(t):=\int_0^t\mathbf E[z(s)^2]\,ds$, from Eq.(7.1), we have $$v'(t)\le \alpha(t)v(t)+\epsilon^2\beta(t)< Av(t)+\epsilon^2B,\ \forall t\in[0,T]$$ for some positive $A, B$ as increasing functions of $T$. That implies $$\mathbf E[z(t,\epsilon)^2]=v'(t)< \epsilon^2 Be^{At},\ \forall t\in[0,T].$$
QED