The brain is a high-dimensional directional network system consisting of many regions as network nodes that exert influences on each other. We aim to reveal whole-brain directional networks based on many subjects' resting-state functional magnetic resonance imaging (fMRI) data. However, it is both statistically and computationally challenging to produce scientifically meaningful estimates of whole-brain directional networks. To address the statistical modeling challenge, we assume modular brain networks, which reflect the brain's functional specialization and integration. We address the computational challenge by developing a variational Bayesian method for the new model. We apply our approach to resting-state fMRI data of many subjects and identify modules and directed connections. The identified modules are accordant with functional brain systems specialized for various functions. In addition, we reveal directed connections between functionally specialized modules. In summary, this paper presents a new computationally efficient and flexible method for directional network studies of the brain and new scientific findings regarding the human brain's functional organization.