Long term observations of black-hole X-ray binaries show that these systems exhibit extreme, aperiodic variability on time scales of few milliseconds to seconds. The observed light-curves display various characteristic features like log-normal distribution and linear rms-flux relation, which indicates that the underlying variability process is stochastic in nature and is thought to be intrinsic to the accretion process. Theoretical models explain this variability as the inward propagating fluctuations of mass accretion rate on viscous timescales, although confirmations from the numerical simulations of magnetized accretion flows are required for a better understanding of the underlying variability process. Using a set of five exceptionally long general relativistic magnetohydrodynamic (GRMHD) simulations of geometrically thick, optically thin, black hole accretion flows as test-beds, we look for hints of propagating fluctuations of accretion rate in the simulation data. Indeed, our results from these simulations show evidence for the inward propagating fluctuations. Our further findings on how these results compare with the propagating fluctuations model and the observations will also be discussed.