I consider the following equation: $$\frac{dy_t}{dt} = a \frac{dx_t}{dt} + x_t +y_t, \tag{1}$$ where $a=$ constant and where $x_t$ follows an Ornstein Uhlenbeck process (see here under Alternative representation for nonstationary processes): $$x_t=x_0 e^{-\theta t} +\mu (1-e^{-\theta t})+ {\sigma\over\sqrt{2\theta}}e^{-\theta t}W_{e^{2\theta t}-1}. \tag{2} $$
For $\sigma = 0$ there is no problem => By substituting Eq.(2) into Eq.(1), we get an ODE.
For $\sigma >0$ there is a problem. The problem is that when I substitute Eq.(2) into Eq.(1), I will have to differentiate the Brownian motion term. Is it possible to do this? If so, how?
Can someone please help me out? I would like to present the solution eventually in for instance in some plot. If useful, numerical methods to solve are also fine for me.
Your relationship between $y_t$ and $x_t$ satisfies: $$\frac{dy_t}{dt} = a \frac{dx_t}{dt} , \tag{1}$$
If $x_t$ were differentiable (which it isn't) then $$y_t = a x_t + b$$ where $b$ is a constant with respect to $t$ would satisfy that equation.
As I mentioned in the comments, you can make your relationship make sense by instead turning it into an integral equation: $$y_t = y_0 + a\int_0^t dx_s$$