Constructing accurate global potential energy surfaces (PESs) describing chemically reactive molecule-molecule collisions of alkali metal dimers presents a major challenge. To be suitable for quantum scattering calculations, such PESs must represent accurately three- and four-body interactions, describe conical intersections, and have a proper asymptotic form at the long range. Here, we demonstrate that such global potentials can be obtained by Gaussian Process (GP) regression merged with the analytic asymptotic expansions at the long range. We propose an efficient sampling technique, which allows us to construct an accurate global PES accounting for different chemical arrangements with <2500 ab initio calculations. We apply this method to (NaK)2 and obtain the first global PES for a system of four alkali metal atoms. The resulting surface exhibits a complex landscape including a pair and a quartet of symmetrically equivalent local minima and a seam of conical intersections. The dissociation energy found from our ab initio calculations is 4534 cm−1. This result is reproduced by the GP models with an error of less than 3%. The GP models of the PES allow us to analyze the features of the global PES, representative of general alkali metal four-atom interactions. Understanding these interactions is of key importance in the field of ultracold chemistry.