The main stream of non-invasive neuroimaging studies has been shifted from identifying functionally specialized regions in the human brain, to elucidating the mechanism of dynamic integration among these regions. To approach the mechanism of integration, discovering brain networks has received growing interest in recent years. Structural and functional brain networks have been studied using diffusion and functional MRI (dMRI and fMRI), respectively, and more recently functional networks are identified by MEG with high temporal resolution. However the brain networks are estimated using either dMRI, fMRI, or MEG in the literature. Here we propose a new method to effectively combine dMRI, fMRI, and MEG for discovering structurally-constrained functional brain networks. The method is based on MEG source reconstruction, jointly estimating source amplitudes and interactions over the whole brain. A multivariate autoregressive (MAR) model are used for formulating directed source interactions in the networks. The MAR model are constrained by structural networks from dMRI data, where the MAR parameters for structurally disconnected source pairs are set to zero and the delay parameters are determined based on structural connectivity lengths. In addition, an fMRI prior on spatial source activities is employed to improve the accuracy and robustness of the method. The dynamic model and prior knowledge are unified within a Bayesian framework using a state-space representation. To jointly estimate the model parameters, the variational Bayesian technique is used for deriving a tractable and computationally efficient algorithm. The proposed method is evaluated on simulated and experimental data. Compared with its non-dynamic counterpart, where the source amplitudes and interactions are separately estimated in a two-stage manner, the proposed method improves quantitative performance measures in simulations, and provides better physiological interpretation and inter-subject consistency in application to face perception data.