NGC7538 IRS1 is considered the best high-mass accretion disk candidate around an O-type young star in the northern hemisphere. We investigated the 3D kinematics and dynamics of circumstellar gas with very high linear resolution, from tens to 1500 AU, with the ultimate goal of building a comprehensive dynamical model for this YSO. We employed four different observing epochs of EVN data at 6.7 GHz, spanning almost eight years, which enabled us to measure, besides line-of-sight (l.o.s.) velocities and positions, also l.o.s. accelerations and proper motions of methanol masers. In addition, we imaged with the JVLA-B array highly-excited ammonia inversion lines, from (6,6) to (13,13), which enabled us to probe the hottest molecular gas very close to the exciting source(s). We found five 6.7 GHz maser clusters which are distributed over a region extended N-S across ~1500 AU and are associated with three peaks of the radio continuum. We proposed that these maser clusters identify three individual high-mass YSOs, named IRS1a, IRS1b, and IRS1c. We modeled the maser clusters in IRS1a and IRS1b in terms of edge-on disks in centrifugal equilibrium. In the first case, masers may trace a quasi-Keplerian thin disk, orbiting around a high-mass YSO, IRS1a, of up to 25 solar masses. This YSO dominates the bolometric luminosity of the region. The second disk is both massive (<16 Msun within ~500 AU) and thick, and the mass of the central YSO, IRS1b, is constrained to be at most a few solar masses. In summary, we present compelling evidence that NGC7538 IRS1 is not forming just one single high-mass YSO, but consists of a multiple system of high-mass YSOs, which are surrounded by accretion disks, and are probably driving individual outflows. This new model naturally explains all the different orientations and disk/outflow structures proposed for the region in previous models.