First consider this system of ODEs. Say I have two variables $u$ and $a$, following

$$ \dot u = -u + f(a) $$

$$ \tau \dot a = -a + gu $$

Where $f$ is nonlinear. If $\tau\ll1$, the timescale of $a$ is fast, and $a$ can be approximated by its steady-state value $gu$, giving a simpler reduced system.

$$ \dot u = -u + f(gu) $$

Now, what happens if the equation for $a$ has additive Gaussian noise? Here $\xi$ denotes the time derivative of a unit-variance Wiener process:

$$ \dot a = [gu-a]/\tau + \sigma \xi $$

If $u$ can be assumed constant relative to the timescales of $a$, this can be approximated as an Ornstein–Uhlenbeck process, and has the steady-state distribution

$$ a \sim \mathcal{N}(gu,\tau\sigma^2/2) $$

I don't have training in stochastic calculus, so I'm wondering how (if it is possible), one might "plug in" this expression for the steady-state distribution of $a$ into the equation for $u$? Since $f$ is nonlinear, we can't simply turn this in to additive or multiplicative noise, and it doesn't seem like Itô's lemma is applicable either.

I think this question is equivalent to asking if there is a good way to treat and reason about the continuous time analogue of this stochastic difference equation:

$$ \frac{\Delta u}{\Delta t} = -u + f(gu + \mathcal{N}(0,\tau\sigma^2/2)) $$

If it matters, I'm particularly interested in the case where the nonlinearity $f$ is the logistic function $f(x) = 1/(1+\exp(-x))$. Approximations are OK, but I would prefer that $u$ be confined to the range of the logistic function $u\in [0,1]$.