Inferring graphical models from high dimensional observations has become an important problem in machine learning and statistics because of its importance in a variety of application domains. One such application is inferring functional connectivity between brain regions from neuroimaging data such as magnetoencephalograpy (MEG) recordings that produce signals with good temporal and spatial resolution. Unfortunately, existing techniques to learn graphical models that have been applied to neuroimaging data have assumed the data to be i.i.d. over time, ignoring key temporal dynamics. Additionally, the signals that arise from neuroimaging data do not exist in isolation as the brain is performing many tasks simultaneously so that most existing methods can introduce spurious connections. We address these issues by introducing a method to learn Gaussian graphical models between multiple time series with latent processes. In addition, we allow for heterogeneity between different groups of MEG recordings by using a hierarchical penalty. The proposed methods are formulated as convex optimization problems that we efficiently solve by developing an alternating directions method of multipliers algorithm. We evaluate the proposed model on synthetic data as well as on global stock index returns and a real MEG data set.