Ultra-large scales close to the cosmological horizon will be probed by the upcoming observational campaigns. They hold the promise to constrain single-field inflation as well as general relativity, but in order to include them in the forthcoming analyses, their modelling has to be robust. In particular, general relativistic effects may be mistaken for primordial signals, and no consensus has emerged either from analytical modelling nor from the numerical route, obstructed by the large dynamical range to be simulated. In this work, we present a numerical technique to overcome the latter limitation: we compute the general relativistic displacement field with the N-body relativistic code gevolution and combine it with the accurate Newtonian simulation Gadget-4. This combination leads to an effective simulation reproducing the desired behaviour at the level of the matter power spectrum and bispectrum. We then measure, for the first time in a simulation, the relativistic scale-dependent bias in Poisson gauge; at redshift z = 0, we find b 1 GR = -8.1 ± 2.8. Our results at the field level are only valid in the Poisson gauge and need to be complemented with a relativistic ray tracing algorithm to compute the number count observable.