Circadian rhythmicity lies at the center of various important physiological and behavioral processes in mammals, such as sleep, metabolism, homeostasis, mood changes and more. It has been shown that this rhythm arises from self-sustained biomolecular oscillations of a neuronal network located in the Suprachiasmatic Nucleus (SCN). Under normal circumstances, this network remains synchronized to the day-night cycle due to signaling from the retina. Misalignment of these neuronal oscillations with the external light signal can disrupt numerous physiological functions and take a long-lasting toll on health and well-being. In this work, we study a modern computational neuroscience model to determine the limits of circadian synchronization to external light signals of different frequency and duty cycle. We employ a matrix-free approach to locate periodic steady states of the high-dimensional model for various driving conditions. Our algorithmic pipeline enables numerical continuation and construction of bifurcation diagrams w.r.t. forcing parameters. We computationally explore the effect of heterogeneity in the circadian neuronal network, as well as the effect of corrective therapeutic interventions, such as that of the drug molecule Longdaysin. Lastly, we employ unsupervised learning to construct a data-driven embedding space for representing neuronal heterogeneity.