The preterm phenotype results from the interplay of multiple disorders affecting the brain and cognitive outcomes. Accurately characterising these interactions can reveal prematurity markers. Bayesian Networks (BNs) are powerful tools to disentangle these relationships, as they inherently measure associations between variables while mitigating confounding factors. We present Modified PC-HC (MPC-HC), a Bayesian Network (BN) structural learning algorithm. MPC-HC employs statistical testing and search-and-score techniques to explore equivalent classes. We employ MPC-HC to estimate BNs for extremely preterm (EP) young adults and full-term controls. Using MRI measurements and cognitive performance markers, we investigate predictive relationships and mutual influences through predictions and sensitivity analysis. We assess the confidence in the estimated BN structures using bootstrapping. Furthermore, MPC-HC’s validation involves assessing its ability to recover benchmark BN structures. MPC-HC achieves an average prediction accuracy of 72.5% compared to 62.5% of PC, 64.5% of MMHC, and 71.5% of HC, while it outperforms PC, MMHC, and HC algorithms in reconstructing the true structure of benchmark BNs. The sensitivity analysis shows that MRI measurements mainly affect EP cognitive scores. Our work has two key contributions: first, the introduction and validation of a new BN structure learning method. Second, demonstrating the potential of BNs in modelling variable relationships, predicting variables of interest, modelling uncertainty, and evaluating how variables impact each other. Finally, we demonstrate this by characterising complex phenotypes, such as preterm birth, and discovering results consistent with literature findings.