Lagrangian bias factors are calculated as responses, to account for modifications to the halo mass function in modified gravity theories and the environmental dependence of screening effects. We show how this formalism should be extended in the case of each screening mechanism. Finally, we cross-check and validate our results in terms of the 2-point statistics against state-of-the-art cosmological N-body simulations that allow us to assess the LPT predictions in the nonlinear, quasi non-linear and baryon acoustic oscillation peak regimes. Comparison with simulations is required to ensure that the predictions will be sufficiently robust for upcoming surveys such as DESI, Euclid, LSST and WFIRST.

We note that a recent paper, submitted to the arXiv while this work was being finalized, also studied biased tracers in MG using CLPT [@Aviles:2018saf] but used a different bias scheme and results were not compared against full simulations.

Our paper is structured as follows: in Section \[sec:Models\], we present the MG models we studied and the details of the N-body simulations employed to test our LPT implementation. We summarize the formalism for developing biased tracer statistics for GR in Section \[sec:LPTDM\]. In Section \[sec:Results:LPT\], we discuss the modifications required to the perturbative schemes to predict the 