Phase-amplitude coupling (PAC) between neural oscillations in different frequency bands has been observed in population recordings of the hippocampus and neocortex [1, 2]. Specifically, bursts of gamma (>30 Hz) oscillations are often found to be phase-locked to a slower oscillation in the theta range (4-8 Hz), raising questions about the mechanisms and computational consequences of this pattern. Recently, inter-areal gamma synchronization was also found to be systematically modulated by theta rhythms . We study mechanisms for these experimental findings and their potential role in long-range inter-areal coupling using spiking neuron network models. Using dynamical systems theory, bursting has previously been analyzed in detail for single neuron dynamics , requiring interaction of fast and slow adapting variables, driving the neuron through alternations of active (limit cycle) and resting (stable equilibrium) periods. On the population level, bursting can be studied within the same framework, treating average synaptic activity and time courses as slow/fast variables. Our model (Panel A) is based on observations of dynamically distinct interneuron populations in the cortex, namely fast-spiking, non-adapting (FS) and adapting interneurons (low-threshold spiking, LTS) exhibiting fast and slow synaptic time courses . Coupled with a population of principal cells (regular spiking, RS), this network displays spontaneous phase-amplitude coupling as sparsely synchronized population bursts. Feed-forward connections between populations of this type lead to theta modulated gamma synchronization (Panel B and C). We study in detail the parameter range/bifurcations of this model. Considering that in-vivo neural circuits are generally only sparsely synchronized, we contrast oscillatory patterns on the population level with stochastic firing of single neurons.