Background Neuroblastoma is a malignant childhood cancer with significant inter- and intrapatient heterogeneity arising from the abnormal differentiation of neural crest cells into sympathetic neurons. The lack of actionable mutations limits therapeutic options, highlighting the need to better understand the molecular mechanisms that drive this differentiation. Although RNA velocity has provided some insights, modeling regulatory relationships is limited. Methods To address this, we applied our integrated gene regulatory network (GRNs) inference (CARDAMOM) and simulation (HARISSA) tools using a published single-cell RNAseq dataset from human sympathoadrenal differentiation. Results Our analysis identified a 97-gene GRN that drives the transition from Schwann cell precursors to chromaffin cells and sympathoblasts, highlighting dynamic interactions such as self-reinforcing loops and toggle switches. The simulation of that GRN was able to reproduce very satisfactorily the experimentally observed gene expression distributions. Conclusions Altogether, these findings demonstrate the utility of our GRN model framework for inferring GRN structure, even in the absence of a time-resolved dataset.