and fluctuating interaction network, respectively. In general, the two white noises are considered uncorrelated. However, they may even correspond to identical time-series as soon as they are shifted with a time lag greater than the correlation time. In such case, the endogenous noise is expected to act with a delay due to its very nature of feedback process, whereas, the additive noise is incorporated immediately, free of signal processing.

By including purely exogenous fluctuations, in the process defined by Eq. (\[ILE0\]), we obtain the following Itô-Langevin equation (ILE) $$\label{ILE}
{\rm d} x \;=\;-\gamma[x-\theta]x^{r-1}{\rm d}t + \mu x^s{\rm d}w + \alpha {\rm d}z \;,$$ where $w,z$ are two [*independent*]{} standard Wiener processes, defined as above, and $\mu,\alpha$ their respective amplitudes. The corresponding Fokker-Planck equation reads $$\label{FP}
\partial_t\rho 
=  \gamma \partial_x( [x-\theta]x^{r-1}\rho) + \partial_{xx}( [\mu^2 x^{2s} + \alpha^2]\rho) \; .$$ Its steady state solution with reflecting boundary conditions is $$\label{FPss}
\rho(x) = \frac{\rho_o}{1+\lambda^2x^{2s}}
\exp[-\gamma_\alpha \int^x \frac{ y^{r-1}(y-\theta) }{1+\lambda^2y^{2s}}{\rm d}y]  \; ,$$ with $\rho_o$ a normalization constant, $\gamma_\alpha\equiv \gamma/\alpha^2$, $\lambda^2\equiv (\mu/\alpha)^2\equiv\gamma_\alpha/\gamma_\mu$. In most cases the integral can be written in terms of hypergeometric functions $_2F_1$[@abram], through $$\label{int}
\int^x \frac{y^{\beta-1}}{1+\lambda^2y^{2s}}{\rm d}y \equiv
\frac{\beta}{x^\beta}
\,_2F_1(c,1,c+1,-\lambda^2x^{2s})$$ with $c\equiv\beta/(2s)\neq -1,-2,\ldots$, whereas, in the marginal case $\beta=0$, we 