The problem of the evolution of a large number of particles due to gravity is crucial to many astrophysical phenomena. An important problem is the dynamical evolution of a dense stellar system, such as a globular cluster (GC), a galactic nucleus (GN) or nuclear star cluster (NSC). Such loci are the breeding grounds of sources of tidal disruptions and gravitational waves. Right in the middle of these regions a massive black hole (MBH) might be lurking, which makes the problem even more interesting, because such massive objects can form a pair and later a binary, which could be powerful source of gravitational radiation for space-borne observatories. The detailed tracking of the dynamical evolution of a set of $N$ stars is a complex problem. Since we lack an analytical solution, it needs to be studied by approximations and numerical methods close to what we might expect from Nature. The close interactions between stars define the core mechanism that determines the global evolution of dense stellar systems. These interactions are responsible for defining the timescale in which catastrophic phenomena happen, such as the core collapse of the system; particularly relevant for the formation of a gravitational capture, that eventually will evolve mostly due to the emission of gravitational radiation. Moreover, depending on the problem we are addressing we might need to add further layers of complexity. For instance, in the case of a GN the presence of gas can play an crucial role, so it needs to be considered, particularly in the massive black hole binary (MBHB) formation process. It has been put forward in the literature that this gas will distribute itself around each MBH in the shape of a disc. The formation of the disc structure around the MBHB is in particular a very important problem which has received very little numerical investigation until the presentation of this work. It is usually assumed that the gas is supplied via the accumulated infall of gaseous clouds on to the binary, and hence this gas is distributed in a disc-like structure around it. Hence, it is relevant to address the formation of binaries taking into account such a gaseous disc around the system in different orbits, and the interaction of the gas with the black holes, not just dynamically, but also via the accretion on to them. Motivated by the complexity and many open question of these fundamental problems, this thesis is (i) a detailed study of the non-linear dynamics that occur in dense stellar systems with state-of-the-art numerical techniques, (ii) a detailed study of the impact of gas on to the binary, in particular to address the role of circumbinary discs on the evolution of a MBHB, and (iii) how repeated infall events of gaseous clouds distribute and shape around such massive binaries, as well as the impact on the dynamical evolution of the binary itself. All of these topics are intertwined and I have worked in them in a parallel way during my PhD. The most remarkable findings of my work are that (i) the use of a softening parameter is critical to analyse the long-term evolution of a dense stellar system, with an important impact on the timescale in which crucial events happen, including the formation of binaries, (ii) the way binaries of MBHs accrete gas in counter-rotating circumbinary discs, will determine the evolution of the massive binary, and (iii) the formation of disc-like structures around these binaries in a GN is, to say the least, challenging. Also, episodic circumbinary structures will modify the orbital evolution of MBHBs, altering their associated gravitational merger timescale.