High-amplitude variability in young stellar objects (YSOs) is usually associated with episodic accretion events. It has not been observed so far in massive YSOs. Here, the high-amplitude variable star sample of Contreras Peña et al. has been used to search for highly variable (ΔK 1 mag) sources coinciding with dense clumps mapped using the 850μm continuum emission by the ATLASGAL survey. A total of 18 variable sources are centered on the submillimeter clump peaks and coincide (<1″) with a 24μm point or compact (<10″) source. Of these 18 sources, 13 can be fit by YSO models. The 13 variable YSOs (VYSOs) have luminosities of ∼10 3 L e , an average mass of 8M e , and a range of ages up to 10 6 yr. A total of 11 of these 13 VYSOs are located in the midst of infrared dark clouds. Nine of the 13 sources have ΔK>2 mag, significantly higher compared to the mean variability of the entire VVV sample. The light curves of these objects sampled between 2010 and 2015 display rising, declining, or quasi-periodic behavior but no clear periodicity. Light-curve analysis using the Plavchan method shows that the most prominent phased signals have periods of a few hundred days. The nature and timescale of variations found in 6.7 Ghz methanol maser emission in massive stars are similar to that of the VYSO light curves. We argue that the origin of the observed variability is episodic accretion. We suggest that the timescale of a few hundred days may represent the frequency at which a spiraling disk feeds dense gas to the young massive star.