The modeling of non-covalent interactions, solvation effects, and chemical reactions in complex molecular environment is a challenging task. Current state-of-the-art approaches often rely on static computations using implicit solvent models and harmonic approximations. Owing to the rich interaction network characteristic of molecular species in the condensed phase, the neglect of the dynamic fluctuations and of the mutual competition with the solvent molecules often leads to an incomplete picture. While ab initio molecular dynamics can, in principle, capture all these effects, the associated computational cost requires compromising the accuracy of the electronic structure method and the statistical convergence of the simulation. The main objective of this thesis is to develop a robust simulation workflow based on neural network potentials and enhanced sampling techniques to describe molecules and reactions in complex environment. The validity and usefulness of the proposed workflow is illustrated on three applications: (a) the solvation environment and protonation state of a mixture of phenol, hydrogen peroxide and methanesulfonic acid relevant to catalyzed industrial phenol hydroxylation, (b) the extension of the above to model the hydroxylation reaction, and (c) the description of a chalcogen-bonded complex in tetrahydrofuran featuring multiple and heavy chemical elements as a potential route for anion transport and catalysis. The developed workflow is illustrated in detail through the first example, for which the solvation environment and protonation states are dominated by strongly competing hydrogen bonds. We demonstrate the relevance of relying upon surrogate neural network potentials, enhanced sampling techniques, and nuclear quantum effects to achieve an accurate description of the complex mixture. The second example extends the simulation pipeline, combining the neural network potentials with metadynamics to enable the modeling of the catalyzed phenol hydroxylation reaction with a special emphasis on its regioselectivity. Finally, the third example broadens the application of neural network potentials to the description of molecular systems containing up to seven chemical elements, which participate in very subtle and specific types of non-covalent interactions. Overall, this work highlights the potential of using neural network-based potentials to model complex reactive systems in condensed phase and to accurately describe competing non-covalent interactions. The developed workflow extends the range of applicability of ab initio molecular dynamics toolbox to real-life chemical systems.