Functional connectivity (FC) has been widely used to study brain network interactions underlying the emerging cognition and behavior of an individual. FC is usually defined as the correlation or partial correlation between brain regions. Although FC is proved to be a good starting point to understand the brain organization, it fails to tell the causal relationship or the direction of interactions. Many directed acyclic graph (DAG) based methods were applied to study the directed interactions using functional magnetic resonance imaging (fMRI) data but the performance was severely limited by the small sample size and high dimensionality, hindering its applications. To overcome the obstacles, we propose a score based joint directed acyclic graph model to estimate the directed FC in fMRI data. Instead of using a combinatorial optimization framework, the structure of DAG is characterized with an algebra equation and further regularized with sparsity and group similarity terms. The simulation results have demonstrated the improved accuracy of the proposed model in detecting causality as compared to other existing methods. In our case-control study of the MIND Clinical Imaging Consortium (MCIC) data, we have successfully identified decreased functional integration, disrupted hub structures and characteristic edges (CtEs) in schizophrenia (SZ) patients. Further comparison between the results from directed FC and undirected FC illustrated the their different emphasis on selected features. We speculate that combining the features from undirected graphical model and directed graphical model might be a promising way to do FC analysis.