We design physics-informed loss functions for training Artificial Neural Network (ANN) models to forecast the ionospheric vertical Total Electron Content (vTEC) from 1 to 24 hours in advance. The ANN models exploit our physics-informed loss functions, data provided by the Global Navigation Satellite Systems (GNSS) receiver installed at Tsukuba (36.06o\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$$^{o}$$\\end{document} N, 140.05o\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$$^{o}$$\\end{document} E), Japan, and external drivers (solar and geomagnetic indices). The time series used span from January 1, 2006, to December 31, 2018, i.e., a full solar cycle. A proper set of external drivers for the ANN models training are selected by ranking their importance in relation to the vTEC dynamics at different forecasting horizons. They result to be the 10.7 cm Solar Flux (F10.7), the magnitude of the Interplanetary Magnetic Field (IMF) BT\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$$\\boldsymbol{B}_{\\boldsymbol{T}}$$\\end{document}, and the Auroral Electrojet (AE) index. Moreover, a second set of indices among those available has been considered as constraints in the design of the physics-informed loss functions. They are the Disturbance Storm Time (Dst) index, the solar wind speed v,BT\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$$\\boldsymbol{v}, \\boldsymbol{B}_{\\boldsymbol{T}}$$\\end{document}, and the By and Bz components of the interplanetary magnetic field. To assess the performance of the resulting ANN models, we use the statistical parameter coefficient of determination (R2\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$$\\boldsymbol{R}^{{\ extbf {2}}}$$\\end{document}), the standard deviation (SD), and the Wilcoxon non-parametric signed ranked test. We show that, in the testing period analyzed (from 2017-09-13, at 04:40:00, to 2018-12-31, at 23:55:00), one of our physics-informed loss functions provides a better performance of the ANN with regard to the standard loss function commonly adopted. In particular, when the new loss function is used in the ANN model, the average SD is minimized across all forecasting horizons in the training, validation and test datasets. SD is 0.2560 TECU, 0.3183 TECU and 0.4240 TECU for the training, validation and test dataset respectively, where 1 TECU = 1016\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$${\ extbf {10}}^{{\ extbf {16}}}$$\\end{document} electrons/m2\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$$\\boldsymbol{m}^{{\ extbf {2}}}$$\\end{document}. The ANN model, incorporating the new loss function and applied to the test dataset, shows a significant improvement according to the Wilcoxon signed ranked test. In fact by selecting a significance level α=0.05\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$$\\boldsymbol{\\alpha }~ \\mathbf {= 0.05}$$\\end{document}, the probability to obtain results by chance with the new loss function as compared to the standard loss function is 0.01504 (i.e., <α\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$$<\\alpha $$\\end{document}), which implies that the new loss function gives a statistical improvement to the forecasting capability of the ANN model. To the best of our knowledge, this is the first time a physics-informed loss function has been designed for the task of forecasting the ionospheric vTEC.