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 ContrerasPe\~{n}a et al.(2016) has been used to search for highly-variable($\Delta$K$\ge$1\,mag) sources coinciding with dense clumps mapped using the 850\mum continuum emission by the ATLASGAL survey. 18 variable sources are centred on the sub-mm clump peaks, and coincide ($<$1") with a 24$\mu$m point or compact ($<$10") source. 13 of these 18 sources can be fit by YSO models. The 13 variable YSOs(VYSO) have luminosities of $\sim$10$^3$ L$_{\odot}$, an average mass of 8 M$_{\odot}$ and a range of ages up to 10$^6$ yr. 11 of these 13 VYSOs are located in the midst of infrared dark clouds. 9 of the 13 sources have $\Delta$K$>$2 mag, significantly higher compared to the mean variability of the entire VVV sample. The light curves of these objects sampled between 2010-2015 display rising, declining, or quasi-periodic behaviour but no clear periodicity. Light-curve analysis using Plavchan method show that the most prominent phased signals have periods of a few hundred days. The nature and time-scale of variations found in 6.7 Ghz methanol maser emission (MME) 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 spiralling disk feeds dense gas to the young massive star.