Solution (source code)

= Solution

Write the source's noise terms consistently as $\eta_R,\eta_P$. Subtract the steady <mean> to get $\dot{\delta R}=-\gamma_R\delta R+\eta_R$. Use the Fourier convention $f(t)=\int d\omega\,e^{-i\omega t}\widetilde f(\omega)/(2\pi)$. After initial transients decay,
$$
\widetilde{\delta R}(\omega)=\frac{\widetilde\eta_R(\omega)}{\gamma_R-i\omega}.
$$
Fourier transforming the stationary white-noise <covariance> gives
$$
\langle\widetilde\eta_R(\omega)\widetilde\eta_R(\omega')\rangle
=2\pi q_R\delta(\omega+\omega').
$$
Therefore the <RNA> fluctuation <power spectral density> is $S_R(\omega)=q_R/(\gamma_R^2+\omega^2)$. Its integral is the equal-time <variance>:
$$
\boxed{\langle\delta R^2\rangle=\int_{-\infty}^\infty\frac{d\omega}{2\pi}
\frac{q_R}{\gamma_R^2+\omega^2}=\frac{q_R}{2\gamma_R}.}
$$
The positivity of $\gamma_R$ ensures a stationary limit. The same calculation gives the <autocovariance> $\langle\delta R(t)\delta R(0)\rangle=[q_R/(2\gamma_R)]e^{-\gamma_R|t|}$: the decay rate filters the white forcing into a finite correlation time. Here the <covariance> convention is $q_R\delta(t-t')$; replacing it by $2D_R\delta(t-t')$ would identify $q_R=2D_R$.