Context. It has been claimed that NGC7538 IRS1 is a high-mass young stellar object (YSO) with 30 M-circle dot, surrounded by a rotating Keplerian disk, probed by a linear distribution of methanol masers. The YSO is also powering a strong compact HII region or ionized wind, and is driving at least one molecular outflow. The axis orientations of the different structures (ionized gas, outflow, and disk) are, however, misaligned, which has led to the different competing models proposed to explain individual structures. Aims. We investigate 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 what is considered the best high-mass accretion disk candidate around an O-type young star in the northern hemisphere. Methods. We used high-angular resolution observations of 6.7 GHz CH3OH masers with the EVN, NH3 inversion lines with the JVLA B-Array, and radio continuum with the VLA A-Array. In particular, we employed four different observing epochs of EVN data at 6.7 GHz, spanning almost eight years, which enabled us to measure line-of-sight (l.o.s.) accelerations and proper motions of CH3OH masers, besides l.o.s. velocities and positions (as done in previous works). In addition, we imaged highly excited NH3 inversion lines, from (6,6) to (13,13), which enabled us to probe the hottest molecular gas very close to the exciting source(s). Results. We confirm previous results that five 6.7 GHz maser clusters (labeled from "A" to "E") are distributed over a region extended N-S across approximate to 1500 AU, and are associated with three components of the radio continuum emission. We propose that these maser clusters identify three individual high-mass YSOs in NGC7538 IRS1, named IRS1a (associated with clusters "B" and "C"), IRS1b (associated with cluster "A"), and IRS1c (associated with cluster "E"). We find that the 6.7 GHz masers distribute along a line, with a regular variation in V-LSR with position, along the major axis of the distribution of maser cluster "A" and the combined clusters "B" + "C". A similar V-LSR gradient (although shallower) is also detected in the NH3 inversion lines. Interestingly, the variation in V-LSR with projected position is not linear but quadratic for both maser clusters. We measure proper motions for 33 maser features, which have an average amplitude (4.8 +/- 0.6 km s(-1)) similar to the variation in V-LSR across the maser cluster, and are approximately parallel to the clusters' elongation axes. By studying the time variation in the maser spectrum, we also derive l.o.s. accelerations for 30 features, with typical amplitude of similar to 10(-3)-10(-2) km s(-1) yr(-1). We modeled the masers in both clusters "A" and "B" + "C" in terms of an edge-on disk in centrifugal equilibrium. Based on our modeling, masers of clusters "B" + "C" may trace a quasi-Keplerian similar to 1 M-circle dot, thin disk, orbiting around a high-mass YSO, IRS1a, of up to approximate to 25 M-circle dot. This YSO dominates the bolometric luminosity of the region. The disk traced by the masers of cluster "A" is both massive (less than or similar to 16 M-circle dot inside a radius of approximate to 500 AU) and thick (opening angle approximate to 45 degrees), and the mass of the central YSO, IRS1b, is constrained to be at most a few M-circle dot. Towards cluster "E", NH3 and 6.7 GHz masers trace more quiescent dynamics than for the other clusters. The presence of a radio continuum peak suggests that the YSO associated with the cluster "E", IRS1c, may be an ionizing, massive YSO as well. Conclusions. We present compelling evidence that NGC 7538 IRS1 is forming not 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.
