We study the generation of EEG rhythms by means of realistically coupled neural mass models. Previous neural mass models were used to model cortical voxels and the thalamus. Interactions between voxels of the same and other cortical areas and with the thalamus were taken into account. Voxels within the same cortical area were coupled (short-range connections) with both excitatory and inhibitory connections, while coupling between areas (long-range connections) was considered to be excitatory only. Short-range connection strengths were modeled by using a connectivity function depending on the distance between voxels. Coupling strength parameters between areas were defined from empirical anatomical data employing the information obtained from probabilistic paths, which were tracked by water diffusion imaging techniques and used to quantify white matter tracts in the brain. Each cortical voxel was then described by a set of 16 random differential equations, while the thalamus was described by a set of 12 random differential equations. Thus, for analyzing the neuronal dynamics emerging from the interaction of several areas, a large system of differential equations needs to be solved. The sparseness of the estimated anatomical connectivity matrix reduces the number of connection parameters substantially, making the solution of this system faster. Simulations of human brain rhythms were carried out in order to test the model. Physiologically plausible results were obtained based on this anatomically constrained neural mass model.