Apical constriction during neural tube closure is driven by cell contractions which are preceded by asynchronous and cell-autonomous Ca2+ flashes, as demonstrated in recent experiments. Disruption of these Ca2+ signals and contractions leads to neural tube defects, such as anencephaly and spina bifida. A good understanding of the two-way mechanochemical coupling of Ca2+ signalling and mechanics remains elusive, while live-cell imaging is difficult. Thus, mathematical modelling is essential but existing models do not exhibit good agreement with experiments. We present two new mechanochemical vertex models of apical constriction during neural tube closure and simulate them using CelluLink, a new user-friendly open-source Python package for vertex modelling. The first one-way mechanochemical model only studies the effect of Ca2+ signalling on cell mechanics. It improves previous models, reproducing some key experimental observations, such as the reduction of the neural plate size to 2%-8% of its initial area. Other novel features of the one-way model is the incorporation of the surface ectoderm and of the experimental amplitude and frequency profiles of the Ca2+ flashes. Furthermore, guided by experiments, the damping coefficient of the vertices and cell-cell adhesion are modelled as functions of the actomyosin concentration and cell size. Furthermore, we develop a two-way model which improves the one-way model by capturing the two-way coupling between Ca2+ signalling and cell mechanics, through the incorporation of stretch-sensitive Ca2+ channels. These channels enable cells to sense mechanical stimuli and encode them into Ca2+ signals. In the two-way model, the Ca2+ flash frequency and amplitude profiles are model outputs and are not inputs as in the one-way model. Finally, we use both models to propose a series of hypotheses for future experiments.