Context. It has been claimed that NGC 7538 IRS1 is a high-mass young stellar object (YSO) with 30 M , 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 CH 3 OH masers with the EVN, NH 3 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 CH 3 OH masers, besides l.o.s. velocities and positions (as done in previous works). In addition, we imaged highly excited NH 3 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 ≈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 NGC 7538 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 NH 3 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 ∼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 ∼1 M , thin disk, orbiting around a high-mass YSO, IRS1a, of up to ≈25 M . This YSO dominates the bolometric luminosity of the region. The d...