The use of Multivariate Granger Causality (MVGC) in estimating directed Blood-Oxygen-Level- Dependant (BOLD) connectivity is still controversial. This is mostly due to the short data lenghts typically available in functional MRI (fMRI) acquisitions, to the very nature of the BOLD acquisition strategy (which yields extremely low signal- to-noise-ratio) and importantly to the fact that neuronal activity is convolved with a slow-varying haemodynamic response function (HRF) which therefore generates a temporal confound which is arduous to account for when basing MVGC estimates on vector autoregressive models (VAR). In this paper, we employ realistic complex network models based on Izhikevich neuronal populations, interlinked by realistic neuronal fiber bundles which exert compounded directed influences and cascade into Baloon-model-like neurovascular coupling, to explore and validate the MVGC approach to directed connectivity estimation in realistic fMRI conditions and in a complex directed network setting. In particular, we show in silico that the top 1 percentile of a BOLD connectivity matrix estimated with MVGC from BOLD data similar to the one provided by the Human Connectome Project (HCP) has a Positive Predictive Value very close to 1, hence corroborating the evidence that the “strongest” connections can be safely studied with this method in fMRI.