$\displaystyle\tilde{\mathbf{D}}_{\mathbf{k},\nu,\nu^{\prime}}=\frac{1}{M}\sum_{\mathbf{i},\mathbf{j}}\mathrm{e}^{-\mathrm{i}\mathbf{k}\cdot\mathbf{R}_{\mathbf{ij}}}\tilde{\mathbf{D}}_{\mathbf{R}_{\mathbf{ij}},\nu,\nu^{\prime}}.$
(104)

Again, it is worth noting that, in principle, for any _real_ material,
$\tilde{T}_{\mathbf{R}_{\mathbf{ij}},\nu,\nu^{\prime}}$ and
$\tilde{\mathbf{D}}_{\mathbf{R}_{\mathbf{ij}},\nu,\nu^{\prime}}$ can be
obtained as standard outputs of _Wannier90_ code [77], given its interfaces to
a certain number of available DFT codes.

##  Appendix B The Houston approach

One of the methods used to simulate the behavior of pumped semiconductors is
the Houston approach [27, 70]. Such an approach is usually formulated in the
velocity gauge and first quantization, just for the reason that will become
clear in the following. Let us start from the time-independent single-particle
Hamiltonian of the Bloch system under analysis,

$\hat{H}_{0}=\frac{\hat{p}^{2}}{2m}+V\left(\hat{\mathbf{r}}\right),$ (105)

where $\hat{\mathbf{p}}$ and $\hat{\mathbf{r}}$ are the momentum and position
operators, respectively, $m$ is the electron mass, and
$V\left(\hat{\mathbf{r}}\right)$ is the periodic potential of the system under
analysis. The related time-independent Schrödinger equation
$\hat{H}_{0}\left|\psi_{\emph{$\boldsymbol{k}$,n}}\right\rangle=\varepsilon_{\emph{$\boldsymbol{k}$,n}}\left|\psi_{\emph{$\boldsymbol{k}$,n}}\right\rangle$
is solved in terms of the Bloch bands
$\varepsilon_{\emph{$\boldsymbol{k}$,n}}$ and of the Bloch functions
$\left|\psi_{\emph{$\boldsymbol{k}$,n}}\right\rangle=\mathrm{e}^{-\mathrm{i}\boldsymbol{k}\cdot\hat{\mathbf{r}}}\left|u_{\emph{$\boldsymbol{k}$,n}}\right\rangle$
where $\left|u_{\emph{$\boldsymbol{k}$,n}}\right\rangle$ displays the same
periodicity of the potential. Accordingly, we have the following reduced
equation
$\hat{H}_{0,\boldsymbol{k}}\left|u_{\emph{$\boldsymbol{k}$,n}}\right\rangle=\varepsilon_{\emph{$\boldsymbol{k}$,n}}\left|u_{\emph{$\boldsymbol{k}$,n}}\right\rangle$
where
$\hat{H}_{0,\boldsymbol{k}}=\mathrm{e}^{\mathrm{i}\boldsymbol{k}\cdot\hat{\mathbf{r}}}\hat{H}_{0}\mathrm{e}^{-\mathrm{i}\boldsymbol{k}\cdot\hat{\mathbf{r}}}=\frac{\left(\hat{\mathbf{p}}-\hbar\boldsymbol{k}\right)^{2}}{2m}+V\left(\hat{\mathbf{r}}\right)$.
Now, if we have a pump pulse described by the vector potential
$\mathbf{A}\left(t\right)$ impinging on the system, the related time-dependent
minimal-coupling Hamiltonian in the velocity gauge reads as

$\hat{H}\left(t\right)=\frac{\left(\hat{\mathbf{p}}+e\mathbf{A}\left(t\right)\right)^{2}}{2m}+V\left(\hat{\mathbf{r}}\right),$
(106)

where $e>0$ is the electronic charge. It is straightforward to demonstrate
that the eigenfunctions and the eigenvalues of this Hamiltonian are simply
$\left|\varphi{}_{\emph{$\boldsymbol{k}$,n}}\left(t\right)\right\rangle=\mathrm{e}^{-\mathrm{i}\boldsymbol{k}\cdot\hat{\mathbf{r}}}\left|u_{\emph{$\boldsymbol{k}$+$\frac{e}{\hbar}$$\mathbf{A}\left(t\right)$,n}}\right\rangle$
and
$\varepsilon_{\emph{$\boldsymbol{k}$+$\frac{e}{\hbar}$$\mathbf{A}\left(t\right)$,n}}$,
respectively. The set of such eigenfunctions is usually named the
instantaneous or the adiabatic basis because these states would exactly
describe the behavior of the system only if the pump pulse would be so slowly
varying on the characteristic timescales (energies) of the system to allow it
to adjust to the pump pulse at each instant of time (i.e., adiabatically).
Accordingly, they do not solve the general time-dependent Schrödinger equation
$\hat{H}\left(t\right)\left|\phi{}_{\boldsymbol{k}}\left(t\right)\right\rangle=\mathrm{i}\hbar\frac{\partial}{\partial
t}\left|\phi{}_{\boldsymbol{k}}\left(t\right)\right\rangle$, but they can be
used as a basis for expanding
$\left|\phi{}_{\boldsymbol{k}}\left(t\right)\right\rangle=\sum_{n}\lambda{}_{\emph{$\boldsymbol{k}$,n}}\left(t\right)\left|\varphi{}_{\emph{$\boldsymbol{k}$,n}}\left(t\right)\right\rangle$.
The projection coefficients
$\lambda{}_{\emph{$\boldsymbol{k}$,n}}\left(t\right)$ are determined via the
following equation of motion,

$\mathrm{i}\hbar\frac{\partial}{\partial
t}\lambda{}_{\emph{$\boldsymbol{k}$,n}}\left(t\right)=\left(\varepsilon_{\emph{$\boldsymbol{k}$+$\frac{e}{\hbar}$$\mathbf{A}\left(t\right)$,n}}-\theta{}_{\emph{$\boldsymbol{k}$,n}}\left(t\right)\right)\lambda{}_{\emph{$\boldsymbol{k}$,n}}\left(t\right)\\\
+\mathrm{i}\hbar\frac{e}{m}\sum_{n^{\prime}\left(\neq
n\right)}\frac{\mathbf{E}\left(t\right)\cdot\mathbf{p}{}_{\boldsymbol{k},n,n^{\prime}}\left(t\right)}{\Delta\varepsilon_{\boldsymbol{k},n,n^{\prime}}\left(t\right)}\lambda{}_{\boldsymbol{k},n^{\prime}}\left(t\right),$
(107)

where
$\theta{}_{\emph{$\boldsymbol{k}$,n}}\left(t\right)=\left\langle\varphi{}_{\emph{$\boldsymbol{k}$,n}}\left(t\right)\right|\mathrm{i}\hbar\frac{\partial}{\partial
t}\left|\varphi{}_{\emph{$\boldsymbol{k}$,n}}\left(t\right)\right\rangle$ is
connected to the Berry phase of the system and can be neglected if there is no
degeneracy,
$\mathbf{p}{}_{\boldsymbol{k},n,n^{\prime}}\left(t\right)=\left\langle\varphi{}_{\emph{$\boldsymbol{k}$,n}}\left(t\right)\right|\hat{\mathbf{p}}\left|\varphi{}_{\boldsymbol{k},n^{\prime}}\left(t\right)\right\rangle$
is the matrix element of the momentum in the instantaneous basis,
$\Delta\varepsilon_{\boldsymbol{k},n,n^{\prime}}\left(t\right)=\varepsilon_{\emph{$\boldsymbol{k}$+$\frac{e}{\hbar}$$\mathbf{A}\left(t\right)$,n}}-\varepsilon_{\boldsymbol{k}+\frac{e}{\hbar}\mathbf{A}\left(t\right),n^{\prime}}$
and $\mathbf{E}\left(t\right)=-\frac{\partial}{\partial
t}\mathbf{A}\left(t\right)$ is the applied electric field. Defining new
coefficients $\beta{}_{\boldsymbol{k},n}\left(t\right)$ such that
$\lambda{}_{\emph{$\boldsymbol{k}$,n}}\left(t\right)=\beta{}_{\emph{$\boldsymbol{k}$,n}}\left(t\right)\mathrm{e}^{-\frac{\mathrm{i}}{\hbar}\int_{-\infty}^{t}\varepsilon_{\emph{$\boldsymbol{k}$+$\frac{e}{\hbar}$$\mathbf{A}\left(t^{\prime}\right)$,n}}dt^{\prime}}$,
the projection-coefficient equation further simplifies to,

$\mathrm{i}\hbar\frac{\partial}{\partial
t}\beta{}_{\emph{$\boldsymbol{k}$,n}}\left(t\right)=-\theta{}_{\emph{$\boldsymbol{k}$,n}}\left(t\right)\beta{}_{\emph{$\boldsymbol{k}$,n}}\left(t\right)\\\
+\mathrm{i}\hbar\frac{e}{m}\sum_{n^{\prime}\left(\neq
n\right)}\frac{\mathbf{E}\left(t\right)\cdot\mathbf{p}{}_{\boldsymbol{k},n,n^{\prime}}\left(t\right)}{\Delta\varepsilon_{\boldsymbol{k},n,n^{\prime}}\left(t\right)}\mathrm{e}^{\frac{\mathrm{i}}{\hbar}\int_{-\infty}^{t}\Delta\varepsilon_{\boldsymbol{k},n,n^{\prime}}\left(t^{\prime}\right)dt^{\prime}}\beta{}_{\boldsymbol{k},n^{\prime}}\left(t\right)$
(108)

and the corresponding basis
$\left\\{\mathrm{e}^{-\frac{\mathrm{i}}{\hbar}\int_{-\infty}^{t}\varepsilon_{\emph{$\boldsymbol{k}$+$\frac{e}{\hbar}$$\mathbf{A}\left(t^{\prime}\right)$,n}}dt^{\prime}}\left|\varphi{}_{\emph{$\boldsymbol{k}$,n}}\left(t\right)\right\rangle\right\\}$
is the Houston basis, which differs just by a time-dependent phase factor from
the instantaneous basis (actually, both basis are often dubbed in the
literature as Houston basis).

The main appeal of such a procedure resides in the possibility of obtaining
sensible results even if one: (i) focuses only on few bands (e.g., one
valence, one conduction and, if needed, one core band), (ii) supposes that
$\mathbf{p}{}_{\boldsymbol{k},n,n^{\prime}}\left(t\right)$ is approximately
$\boldsymbol{k}$ independent so that
$\mathbf{p}{}_{\boldsymbol{k},n,n^{\prime}}\left(t\right)\approx\mathbf{p}{}_{n,n^{\prime}}$
(the time dependence gets cancelled as well), and (iii) uses the parabolic
approximation
$\varepsilon_{\emph{$\boldsymbol{k}$,n}}\approx\varepsilon_{n}+\frac{\hbar^{2}k^{2}}{2m_{n}}$,
that is, retains only relevant gaps
$\left|\varepsilon_{n}-\varepsilon_{n^{\prime}}\right|$ and effective masses
$m_{n}$ in the proximity of few selected $\boldsymbol{k}$ points.

|  |
---|---|---
|  |
Figure 15: Top and bottom panels: the dimensionless spectral
function,$\bar{A}_{\mathbf{k}}^{R}\left(\bar{\omega},t\right)$, and occupied
spectral function, $\bar{A}_{\mathbf{k}}^{<}\left(\bar{\omega},t\right)$,
respectively, as function of $\omega$ along the _main_ path, for the same
pump-pulse and system parameters as the one of Fig. 6. In the left, middle and
right panels, the time $t$ is chosen to be well before the application of pump
pulse, in the center of the pump pulse ($t=0$), and well after the application
of pump pulse, respectively. The black solid curves show the equilibrium
bands. The numerical value of $0^{+}$ has been chosen to be $0.2\Delta/\hbar$.

##  Appendix C Out-of-Equilibrium Spectral Functions

To obtain the spectral functions, we need the Fourier transformation of the
GFs with respect to time, which we perform as follows,

$\displaystyle G_{\mathbf{k}}^{R,<}\left(\omega,t\right)$
$\displaystyle=\int_{-\infty}^{+\infty}d\tau
e^{i\omega\tau-0^{+}\left|\tau\right|}G_{\mathbf{k}}^{R,<}\left(t+\frac{\tau}{2},t-\frac{\tau}{2}\right),$
(109)

where $0^{+}$ is an infinitesimal convergence factor. Then, the retarded
spectral function is given by

$A_{\mathbf{k}}^{R}\left(\omega,t\right)=-\frac{1}{\pi}\imaginary\left[\Tr
G_{\mathbf{k}}^{R}\left(\omega,t\right)\right],$ (110)

while the lesser spectral function is defined as

$A_{\mathbf{k}}^{<}\left(\omega,t\right)=\frac{1}{2\pi}\imaginary\left[\Tr
G_{\mathbf{k}}^{<}\left(\omega,t\right)\right].$ (111)

In Fig. 15 top and bottom panels, we report the dimensionless retarded and
lesser functions,
$\bar{A}_{\mathbf{k}}^{R,<}\left(\bar{\omega},t\right)=A_{\mathbf{k}}^{R,<}\left(\bar{\omega},t\right)\Delta/\hbar$
, respectively, as function of $\bar{\omega}$ along the _main_ path, for the
same pump-pulse and system parameters as the ones of Fig. 6. In the left,
middle and right panels, the time $t$ is chosen to be well before the
application of the pump pulse, in the center of the pump pulse ($t=0$), and
well after the application of pump pulse, respectively. Clearly, during the
application of the pump pulse, the spectral functions,
$\bar{A}_{\mathbf{k}}^{R,<}\left(\bar{\omega},t\right)$, become negative and
lose their original physical interpretations.

Obviously, in the absence of the pump pulse (both before and after its
application), the spectral function gives correct information about the energy
bands of the system. The left panel of Fig. 15 can be directly compared to
Fig. 3 and the only difference to be acknowledged is that, in the former, the
band broadening originates from the finite numerical value of $0^{+}$, while,
in the latter, it originates from the finite FWHM of the probe pulse. After
the application of the pump pulse, the bands recover their equilibrium shape
as it can be seen by comparing the top-right panel of Fig. 15 to its top-left
panel. On the other hand, the occupied spectral function (see bottom-right
panel of Fig. 15) shows that at some k points we have residual excitations
similarly to what is reported in the top panel of Fig. 7.

## References

  * Brabec and Krausz [2000] T. Brabec and F. Krausz, Reviews of Modern Physics 72, 545 (2000).
  * Krausz and Ivanov [2009] F. Krausz and M. Ivanov, Reviews of modern physics 81, 163 (2009).
  * Krausz and Stockman [2014] F. Krausz and M. I. Stockman, Nature Photonics 8, 205 (2014).
  * Calegari  _et al._ [2016] F. Calegari, G. Sansone, S. Stagira, C. Vozzi, and M. Nisoli, Journal of Physics B: Atomic, Molecular and Optical Physics 49, 062001 (2016).
  * Borrego-Varillas  _et al._ [2022] R. Borrego-Varillas, M. Lucchini, and M. Nisoli, Reports on Progress in Physics 85, 066401 (2022).
  * Zürch  _et al._ [2017a] M. Zürch, H.-T. Chang, L. J. Borja, P. M. Kraus, S. K. Cushing, A. Gandman, C. J. Kaplan, M. H. Oh, J. S. Prell, D. Prendergast, C. D. Pemmaraju, D. M. Neumark, and S. R. Leone, Nature Communications 8, 15734 (2017a).
  * Zürch  _et al._ [2017b] M. Zürch, H.-T. Chang, P. M. Kraus, S. K. Cushing, L. J. Borja, A. Gandman, C. J. Kaplan, M. H. Oh, J. S. Prell, D. Prendergast, C. D. Pemmaraju, D. M. Neumark, and S. R. Leone, Structural Dynamics 4, 044029 (2017b).
  * Kaplan  _et al._ [2018] C. J. Kaplan, P. M. Kraus, A. D. Ross, M. Zürch, S. K. Cushing, M. F. Jager, H.-T. Chang, E. M. Gullikson, D. M. Neumark, and S. R. Leone, Phys. Rev. B 97, 205202 (2018).
  * Perfetti  _et al._ [2008] L. Perfetti, P. A. Loukakos, M. Lisowski, U. Bovensiepen, M. Wolf, H. Berger, S. Biermann, and A. Georges, New Journal of Physics 10, 053019 (2008).
  * Schmitt  _et al._ [2008] F. Schmitt, P. S. Kirchmann, U. Bovensiepen, R. G. Moore, L. Rettig, M. Krenz, J.-H. Chu, N. Ru, L. Perfetti, D. Lu, _et al._ , Science 321, 1649 (2008).
  * von der Linde  _et al._ [1995] D. von der Linde, T. Engers, G. Jenke, P. Agostini, G. Grillon, E. Nibbering, A. Mysyrowicz, and A. Antonetti, Physical Review A 52, R25 (1995).
  * Norreys  _et al._ [1996] P. A. Norreys, M. Zepf, S. Moustaizis, A. P. Fews, J. Zhang, P. Lee, M. Bakarezos, C. N. Danson, A. Dyson, P. Gibbon, P. Loukakos, D. Neely, F. N. Walsh, J. S. Wark, and A. E. Dangor, Phys. Rev. Lett. 76, 1832 (1996).
  * Chin  _et al._ [2001] A. H. Chin, O. G. Calderón, and J. Kono, Physical Review Letters 86, 3292 (2001).
  * Ghimire  _et al._ [2011] S. Ghimire, A. D. DiChiara, E. Sistrunk, P. Agostini, L. F. DiMauro, and D. A. Reis, Nature physics 7, 138 (2011).
  * Luu  _et al._ [2015] T. T. Luu, M. Garg, S. Y. Kruchinin, A. Moulet, M. T. Hassan, and E. Goulielmakis, Nature 521, 498 (2015).
  * Han  _et al._ [2016] S. Han, H. Kim, Y. W. Kim, Y.-J. Kim, S. Kim, I.-Y. Park, and S.-W. Kim, Nature communications 7, 13105 (2016).
  * Borja  _et al._ [2016] L. J. Borja, M. Zürch, C. D. Pemmaraju, M. Schultze, K. Ramasesha, A. Gandman, J. S. Prell, D. Prendergast, D. M. Neumark, and S. R. Leone, J. Opt. Soc. Am. B 33, C57 (2016).
  * You  _et al._ [2017] Y. S. You, D. A. Reis, and S. Ghimire, Nature physics 13, 345 (2017).
  * Liu  _et al._ [2017] H. Liu, Y. Li, Y. S. You, S. Ghimire, T. F. Heinz, and D. A. Reis, Nature Physics 13, 262 (2017).
  * Jiang  _et al._ [2019] S. Jiang, S. Gholam-Mirzaei, E. Crites, J. E. Beetar, M. Singh, R. Lu, M. Chini, and C. Lin, Journal of Physics B: Atomic, Molecular and Optical Physics 52, 225601 (2019).
  * Schultze  _et al._ [2013] M. Schultze, E. M. Bothschafter, A. Sommer, S. Holzner, W. Schweinberger, M. Fiess, M. Hofstetter, R. Kienberger, V. Apalkov, V. S. Yakovlev, _et al._ , Nature 493, 75 (2013).
  * Stojchevska  _et al._ [2014] L. Stojchevska, I. Vaskivskyi, T. Mertelj, P. Kusar, D. Svetin, S. Brazovskii, and D. Mihailovic, Science 344, 177 (2014).
  * Schultze  _et al._ [2014] M. Schultze, K. Ramasesha, C. Pemmaraju, S. Sato, D. Whitmore, A. Gandman, J. S. Prell, L. Borja, D. Prendergast, K. Yabana, _et al._ , Science 346, 1348 (2014).
  * Lucchini  _et al._ [2016] M. Lucchini, S. A. Sato, A. Ludwig, J. Herrmann, M. Volkov, L. Kasmi, Y. Shinohara, K. Yabana, L. Gallmann, and U. Keller, Science 353, 916 (2016).
  * Mashiko  _et al._ [2016] H. Mashiko, K. Oguri, T. Yamaguchi, A. Suda, and H. Gotoh, Nature physics 12, 741 (2016).
  * Zürch  _et al._ [2017] M. Zürch, H.-T. Chang, P. M. Kraus, S. K. Cushing, L. J. Borja, A. Gandman, C. J. Kaplan, M. H. Oh, J. S. Prell, D. Prendergast, _et al._ , Structural Dynamics 4, 044029 (2017).
  * Schlaepfer  _et al._ [2018] F. Schlaepfer, M. Lucchini, S. A. Sato, M. Volkov, L. Kasmi, N. Hartmann, A. Rubio, L. Gallmann, and U. Keller, Nature Physics 14, 560 (2018).
  * Kaplan  _et al._ [2019] C. J. Kaplan, P. M. Kraus, E. M. Gullikson, L. J. Borja, S. K. Cushing, M. Zürch, H.-T. Chang, D. M. Neumark, and S. R. Leone, J. Opt. Soc. Am. B 36, 1716 (2019).
  * Geneaux  _et al._ [2019] R. Geneaux, H. J. Marroux, A. Guggenmos, D. M. Neumark, and S. R. Leone, Philosophical Transactions of the Royal Society A 377, 20170463 (2019).
  * Inzani  _et al._ [2022] G. Inzani, L. Adamska, A. Eskandari-asl, N. D. Palo, G. L. Dolso, B. Moio, L. J. D’Onofrio, A. Lamperti, A. Molle, R. Borrego-Varillas, M. Nisoli, S. Pittalis, C. A. Rozzi, A. Avella, and M. Lucchini, arXiv preprint arXiv:2212.02157 (2022).
  * Neufeld  _et al._ [2023] O. Neufeld, N. Tancogne-Dejean, U. De Giovannini, H. Hübener, and A. Rubio, npj Computational Materials 9, 39 (2023).
  * Rohwer  _et al._ [2011] T. Rohwer, S. Hellmann, M. Wiesenmayer, C. Sohrt, A. Stange, B. Slomski, A. Carr, Y. Liu, L. M. Avila, M. Kalläne, _et al._ , Nature 471, 490 (2011).
  * Smallwood  _et al._ [2012] C. L. Smallwood, J. P. Hinton, C. Jozwiak, W. Zhang, J. D. Koralek, H. Eisaki, D.-H. Lee, J. Orenstein, and A. Lanzara, Science 336, 1137 (2012).
  * Hellmann  _et al._ [2012] S. Hellmann, T. Rohwer, M. Kalläne, K. Hanff, C. Sohrt, A. Stange, A. Carr, M. Murnane, H. Kapteyn, L. Kipp, _et al._ , Nature communications 3, 1069 (2012).
  * Papalazarou  _et al._ [2012] E. Papalazarou, J. Faure, J. Mauchain, M. Marsi, A. Taleb-Ibrahimi, I. Reshetnyak, A. van Roekeghem, I. Timrov, N. Vast, B. Arnaud, and L. Perfetti, Phys. Rev. Lett. 108, 256808 (2012).
  * Wang  _et al._ [2013] Y. Wang, H. Steinberg, P. Jarillo-Herrero, and N. Gedik, Science 342, 453 (2013).
  * Johannsen  _et al._ [2013] J. C. Johannsen, S. Ulstrup, F. Cilento, A. Crepaldi, M. Zacchigna, C. Cacho, I. C. E. Turcu, E. Springate, F. Fromm, C. Raidel, T. Seyller, F. Parmigiani, M. Grioni, and P. Hofmann, Phys. Rev. Lett. 111, 027403 (2013).
  * Rameau  _et al._ [2016] J. Rameau, S. Freutel, A. Kemper, M. A. Sentef, J. Freericks, I. Avigo, M. Ligges, L. Rettig, Y. Yoshida, H. Eisaki, _et al._ , Nature communications 7, 13761 (2016).
  * Reimann  _et al._ [2018] J. Reimann, S. Schlauderer, C. Schmid, F. Langer, S. Baierl, K. Kokh, O. Tereshchenko, A. Kimura, C. Lange, J. Güdde, _et al._ , Nature 562, 396 (2018).
  * Smith and Traum [1975] N. Smith and M. Traum, Physical Review B 11, 2087 (1975).
  * Himpsel and Eastman [1978] F. J. Himpsel and D. E. Eastman, Phys. Rev. B 18, 5236 (1978).
  * Kampf and Schrieffer [1990] A. P. Kampf and J. R. Schrieffer, Phys. Rev. B 42, 7967 (1990).
  * Smith and Kevan [1991] K. E. Smith and S. D. Kevan, Progress in solid state chemistry 21, 49 (1991).
  * Damascelli  _et al._ [2003] A. Damascelli, Z. Hussain, and Z.-X. Shen, Reviews of modern physics 75, 473 (2003).
  * Hüfner [2013] S. Hüfner,  _Photoelectron spectroscopy: principles and applications_ (Springer Science & Business Media, 2013).
  * Sobota  _et al._ [2021] J. A. Sobota, Y. He, and Z.-X. Shen, Reviews of Modern Physics 93, 025006 (2021).
  * Bovensiepen and Kirchmann [2012] U. Bovensiepen and P. S. Kirchmann, Laser & Photonics Reviews 6, 589 (2012).
  * Smallwood  _et al._ [2016] C. L. Smallwood, R. A. Kaindl, and A. Lanzara, EPL (Europhysics Letters) 115, 27001 (2016).
  * Zhou  _et al._ [2018] X. Zhou, S. He, G. Liu, L. Zhao, L. Yu, and W. Zhang, Reports on Progress in Physics 81, 062101 (2018).
  * Weinelt [2002] M. Weinelt, Journal of Physics: Condensed Matter 14, R1099 (2002).
  * Echenique  _et al._ [2004] P. M. Echenique, R. Berndt, E. V. Chulkov, T. Fauster, A. Goldmann, and U. Höfer, Surface Science Reports 52, 219 (2004).
  * Sobota  _et al._ [2013] J. A. Sobota, S.-L. Yang, A. F. Kemper, J. J. Lee, F. T. Schmitt, W. Li, R. G. Moore, J. G. Analytis, I. R. Fisher, P. S. Kirchmann, T. P. Devereaux, and Z.-X. Shen, Phys. Rev. Lett. 111, 136802 (2013).
  * Broers and Mathey [2022] L. Broers and L. Mathey, Physical Review Research 4, 013057 (2022).
  * Avigo  _et al._ [2013] I. Avigo, R. Cortés, L. Rettig, S. Thirupathaiah, H. S. Jeevan, P. Gegenwart, T. Wolf, M. Ligges, M. Wolf, J. Fink, _et al._ , Journal of Physics: Condensed Matter 25, 094003 (2013).
  * Yang  _et al._ [2015] S.-L. Yang, J. A. Sobota, D. Leuenberger, Y. He, M. Hashimoto, D. H. Lu, H. Eisaki, P. S. Kirchmann, and Z.-X. Shen, Phys. Rev. Lett. 114, 247001 (2015).
  * De Giovannini  _et al._ [2013] U. De Giovannini, G. Brunetto, A. Castro, J. Walkenhorst, and A. Rubio, ChemPhysChem 14, 1298 (2013).
  * De Giovannini  _et al._ [2016] U. De Giovannini, H. Hubener, and A. Rubio, Nano letters 16, 7993 (2016).
  * Wopperer  _et al._ [2017] P. Wopperer, U. De Giovannini, and A. Rubio, The European Physical Journal B 90, 1 (2017).
  * De Giovannini  _et al._ [2017] U. De Giovannini, H. Hübener, and A. Rubio, Journal of chemical theory and computation 13, 265 (2017).
  * Pemmaraju  _et al._ [2018] C. D. Pemmaraju, F. D. Vila, J. J. Kas, S. A. Sato, J. J. Rehr, K. Yabana, and D. Prendergast, Computer Physics Communications 226, 30 (2018).
  * Tancogne-Dejean and Rubio [2018] N. Tancogne-Dejean and A. Rubio, Science advances 4, eaao5207 (2018).
  * De Giovannini  _et al._ [2018] U. De Giovannini, A. Castro, _et al._ , Attosecond Molecular Dynamics 13, 424 (2018).
  * Armstrong  _et al._ [2021] G. S. Armstrong, M. A. Khokhlova, M. Labeye, A. S. Maxwell, E. Pisanty, and M. Ruberti, The European Physical Journal D 75, 209 (2021).
  * Schüler  _et al._ [2021] M. Schüler, J. A. Marks, Y. Murakami, C. Jia, and T. P. Devereaux, Physical Review B 103, 155409 (2021).
  * Freericks  _et al._ [2006] J. K. Freericks, V. M. Turkowski, and V. Zlatić, Phys. Rev. Lett. 97, 266408 (2006).
  * Aoki  _et al._ [2014] H. Aoki, N. Tsuji, M. Eckstein, M. Kollar, T. Oka, and P. Werner, Reviews of Modern Physics 86, 779 (2014).
  * Golež  _et al._ [2019] D. Golež, L. Boehnke, M. Eckstein, and P. Werner, Physical Review B 100, 041111(R) (2019).
  * Eckstein [2021] M. Eckstein, Journal of Electron Spectroscopy and Related Phenomena 253, 147108 (2021).
  * DÓnofrio  _et al._ [2022] L. J. DÓnofrio, A. Eskandari-asl, and A. Avella, Journal of Magnetism and Magnetic Materials 546, 168785 (2022).
  * Sato  _et al._ [2018] S. A. Sato, M. Lucchini, M. Volkov, F. Schlaepfer, L. Gallmann, U. Keller, and A. Rubio, Physical Review B 98, 035202 (2018).
  * Mancini and Avella [2004] F. Mancini and A. Avella, Advances in Physics 53, 537 (2004).
  * Avella and Mancini [2007] A. Avella and F. Mancini, Phys. Rev. B 75, 134518 (2007).
  * Avella and Mancini [2012] A. Avella and F. Mancini, The composite operator method (com), in  _Strongly Correlated Systems: Theoretical Methods_, edited by A. Avella and F. Mancini (Springer Berlin Heidelberg, Berlin, Heidelberg, 2012) pp. 103–141.
  * Avella [2014a] A. Avella, The European Physical Journal B 87, 1 (2014a).
  * Avella [2014b] A. Avella, Advances in Condensed Matter Physics 2014, 515698 (2014b).
  * Di Ciolo and Avella [2018] A. Di Ciolo and A. Avella, Condensed Matter Physics 21 (2018).
  * Mostofi  _et al._ [2014] A. A. Mostofi, J. R. Yates, G. Pizzi, Y.-S. Lee, I. Souza, D. Vanderbilt, and N. Marzari, Computer Physics Communications 185, 2309 (2014).
  * Peierls [1933] R. Peierls, Zeitschrift für Physik 80, 763 (1933).
  * Ismail-Beigi  _et al._ [2001] S. Ismail-Beigi, E. K. Chang, and S. G. Louie, Physical Review Letters 87, 087402 (2001).
  * Yue and Gaarde [2022] L. Yue and M. B. Gaarde, JOSA B 39, 535 (2022).
  * Freericks  _et al._ [2009] J. K. Freericks, H. R. Krishnamurthy, and T. Pruschke, Phys. Rev. Lett. 102, 136401 (2009).
  * Sentef  _et al._ [2013] M. Sentef, A. F. Kemper, B. Moritz, J. K. Freericks, Z.-X. Shen, and T. P. Devereaux, Physical Review X 3, 041033 (2013).
  * Freericks  _et al._ [2015] J. Freericks, H. Krishnamurthy, M. Sentef, and T. Devereaux, Physica Scripta 2015, 014012 (2015).
  * Schüler and Sentef [2021] M. Schüler and M. A. Sentef, Journal of Electron Spectroscopy and Related Phenomena 253, 147121 (2021).


#  Non-adiabatic Level Crossing in Resonant Neutrino Oscillations

Stephen J. Parke  Theoretical Physics Department,
Fermi National Accelerator Laboratory,
Batavia, IL 60510, USA

(May 27, 1986 — FERMILAB-Pub-86/67-T)

###### Abstract

Analytic results are presented for the probability of detecting an electron
neutrino after passage through a resonant oscillation region. If the electron
neutrino is produced far above the resonance density, this probability is
simply given by

$\langle\,P_{\nu_{e}}\,\rangle\approx\sin^{2}\theta_{0}+P_{\text{\bf x}}\cos
2\theta_{0}\,,$

where $\theta_{0}$ is the vacuum mixing angle. The probability is averaged
over the production as well as the detection positions of the neutrino and
$P_{\text{\bf x}}$ is the Landau-Zener transition probability between
adiabatic states. Finally, this result is applied to resonance oscillations
within the solar interior.

††preprint: https://inspirehep.net/literature/229055
††preprint:  DOI: 10.1103/PhysRevLett.57.1275

Recently Mikheyev and SmirnovMikheev and Smirnov (1986) and BetheBethe (1986)
have revived interest in the solar-neutrino deficit by demonstrating that
electron neutrinos produced in the sun can be efficiently rotated into muon
neutrinos by passage through a resonant oscillation region. This mechanism may
solve the solar-neutrino puzzle. In this paper, I present an analytic result
for the probability of detecting an electron neutrino after passage through
one or more resonant oscillation regions. This result is then used to show the
regions of parameter space, the difference of the squared masses versus the
vacuum mixing angle, for which the solar-neutrino puzzle is solved.

A neutrino state is assumed to be a linear combination of the two flavor
states $|\nu_{e}\rangle$ and $\ket{\nu_{\mu}}$:

$\ket{\nu(t)}=C_{e}(t)\ket{\nu_{e}}+C_{\mu}(t)\ket{\nu_{\mu}}.$ (1)

If the neutrinos are massive, then the mass eigenstates need not be identical
to the flavor eigenstates, so that the Dirac equation which governs the
evolution of the neutrino state is not necessarily diagonal in the flavor
basis. This leads to the well known phenomena of vacuum neutrino oscillations.
In the presence of matter, the non-diagonal nature of this evolution is
further enhanced by coherent forward scattering which can lead to resonant
neutrino oscillations. WolfensteinWolfenstein (1978, 1979) has derived the
Dirac equation for this process, in the ultra-relativistic limit, in terms of
the vacuum mass eigenstates. Here, I use his result, in the flavor basis,
after discarding a term proportional to the identity matrix, as this term only
contributes an overall phase factor to the state $\ket{\nu(t)}$. The resulting
Schrodinger-type wave equation is

$i\frac{d}{dt}\left(\begin{array}[]{c}C_{e}\\\
C_{\mu}\end{array}\right)=\frac{1}{2}\left(\begin{array}[]{cc}-{\Delta_{0}}\cos
2\theta_{0}+\sqrt{2}G_{F}N&{\Delta_{0}}\sin 2\theta_{0}\\\ {\Delta_{0}}\sin
2\theta_{0}&{\Delta_{0}}\cos
2\theta_{0}-\sqrt{2}G_{F}N\end{array}\right)\left(\begin{array}[]{c}C_{e}\\\
C_{\mu}\end{array}\right),$ (2)

where $\Delta_{0}=(m^{2}_{2}-m^{2}_{1})/2k$, $m_{1},\,m_{2}$ are the neutrino
masses, k is the neutrino momentum, $\theta_{0}$ is the vacuum mixing angle,
$N$ is number density of electrons, and $G_{F}$ is the Fermi constant. The
constraints $\Delta_{0}>0$ and $\theta_{0}<\pi/4$ are assumed. At an electron
density, $N$, the matter mass eigenstates are

$\displaystyle\ket{\nu_{1},N}$ $\displaystyle~{}=~{}$
$\displaystyle\cos\theta_{N}~{}\ket{\nu_{e}}~{}-~{}\sin\theta_{N}~{}\ket{\nu_{\mu}}\,,$
$\displaystyle\ket{\nu_{2},N}$ $\displaystyle~{}=~{}$
$\displaystyle\sin\theta_{N}~{}\ket{\nu_{e}}~{}+~{}\cos\theta_{N}~{}\ket{\nu_{\mu}}\,,$
(3)

which have eigenvalues $\pm\Delta_{N}/2$, where

$\Delta_{N}=[(\Delta_{0}\cos
2\theta_{0}-\sqrt{2}G_{F}N)^{2}+\Delta^{2}_{0}\sin^{2}2\theta_{0}]^{1/2}\,,$
(4)

and $\theta_{N}$ satisfies

$\Delta_{N}\sin 2\theta_{N}=\Delta_{0}\sin 2\theta_{0}\,.$ (5)

These states evolve in time by the multiplication of a phase factor, if the
electron density is a constant. For such a constant density there are three
regions of interest: (i) Well below resonance,
$\sqrt{2}G_{F}N\ll\Delta_{0}\cos 2\theta_{0}$, where the matter mixing angle
is $\theta_{N}\sim\theta_{0}$ and the oscillation length is
$L_{0}=2\pi/\Delta_{0}$. Typically, this is the region that the electron
neutrinos are detected in. (ii) At resonance, $\sqrt{2}G_{F}N=\Delta_{0}\cos
2\theta_{0}$, where the matter mixing angle is $\theta_{N}=\pi/4$ and the
resonant oscillation length is $L_{R}=L_{0}/\sin 2\theta_{0}$, which for small
vacuum mixing angle can be many times the vacuum oscillation length. (iii) Far
above resonance, $\sqrt{2}G_{F}N\gg\Delta_{0}\cos 2\theta_{0}$, where the
matter mixing angle $\theta_{N}\sim\pi/2$, and the oscillation length
$L_{N}=2\pi/\Delta_{N}$ is much smaller than the vacuum oscillation length
$L_{0}$. For the situation of current interest the electron neutrinos are
produced above resonance, pass through resonance, and are detected in the
vacuum.

If the electron density varies slowly, the states which evolve independently
in time (the adiabatic states) are

$\exp(-i\,\frac{1}{2}\int^{t}\Delta_{N}dt)~{}\ket{\nu_{1},N(t)}$

and

$\exp(+i\,\frac{1}{2}\int^{t}\Delta_{N}dt)~{}\ket{\nu_{2},N(t)}\,.$

Therefore, it is convenient to use these states, as the basis states, in the
region for which there are no transitions (away from the resonance region). As
a neutrino goes through resonance these adiabatic states may be mixed, but on
the other side of resonance, the neutrino state can still be written as a
linear combination of these states. That is, a basis state produced at time
$t$, going through resonance at time $t_{r}$, and detected at time
$t^{\prime}$ is described by

$\displaystyle\exp(-i\,\frac{1}{2}\int^{t}_{t_{r}}\Delta_{N}dt)~{}\ket{\nu_{1},N(t)}\rightarrow$
$\displaystyle~{}a_{1}~{}\exp(-i\,\frac{1}{2}\int^{t^{\prime}}_{t_{r}}\Delta_{N}dt)~{}\ket{\nu_{1},N(t^{\prime})}~{}+~{}a_{2}~{}\exp(+i\,\frac{1}{2}\int^{t^{\prime}}_{t_{r}}\Delta_{N}dt)~{}\ket{\nu_{2},N(t^{\prime})}\,,$
$\displaystyle\exp(+i\,\frac{1}{2}\int^{t}_{t_{r}}\Delta_{N}dt)~{}\ket{\nu_{2},N(t)}\rightarrow$
$\displaystyle-a^{\ast}_{2}~{}\exp(-i\,\frac{1}{2}\int^{t^{\prime}}_{t_{r}}\Delta_{N}dt)~{}\ket{\nu_{1},N(t^{\prime})}~{}+~{}a^{\ast}_{1}~{}\exp(+i\,\frac{1}{2}\int^{t^{\prime}}_{t_{r}}\Delta_{N}dt)~{}\ket{\nu_{2},N(t^{\prime})}\,,$

where $a_{1}$ and $a_{2}$ are complex numbers such that
$|a_{1}|^{2}+|a_{2}|^{2}=1$. The relationship between the coefficients, for
these two basis states, is due to the special nature of the wave equation,
Eq.(2). The phase factors have been chosen so that coefficients $a_{1}$ and
$a_{2}$ are characteristics of the transitions at resonance and are not
related to the production and detection of the neutrino state.

Hence, the amplitude for producing, at time $t$, and detecting, at time
$t^{\prime}$, an electron neutrino after passage through resonance is

$\displaystyle
A_{1}(t)~{}\exp{-i\,\frac{1}{2}\int^{t^{\prime}}_{t_{r}}\Delta_{N}dt}~{}+~{}A_{2}(t)~{}\exp{+i\,\frac{1}{2}\int^{t^{\prime}}_{t_{r}}\Delta_{N}dt}\,,$

where

$\displaystyle A_{1}(t)$ $\displaystyle=$
$\displaystyle~{}\cos\theta_{0}~{}\biggr{[}a_{1}\cos\theta_{N}~{}\exp{+i\,\frac{1}{2}\int^{t}_{t_{r}}\Delta_{N}dt}~{}-~{}a^{\ast}_{2}\sin\theta_{N}~{}\exp{-i\,\frac{1}{2}\int^{t}_{t_{r}}\Delta_{N}dt}\biggr{]}\,,$
$\displaystyle A_{2}(t)$ $\displaystyle=$
$\displaystyle~{}\sin\theta_{0}~{}\biggr{[}a_{2}\cos\theta_{N}~{}\exp{+i\,\frac{1}{2}\int^{t}_{t_{r}}\Delta_{N}dt}~{}+~{}a^{\ast}_{1}\sin\theta_{N}~{}\exp{-i\,\frac{1}{2}\int^{t}_{t_{r}}\Delta_{N}dt}\biggr{]}.$

Thus the probability of detecting this neutrino as an electron neutrino is
given by

$\displaystyle P_{\nu_{e}}(t,t^{\prime})$ $\displaystyle~{}=~{}$
$\displaystyle|A_{1}(t)|^{2}~{}+~{}|A_{2}(t)|^{2}~{}+~{}2|A_{1}(t)A_{2}(t)|\cos(\int^{t^{\prime}}_{t_{r}}\Delta_{N}dt~{}+~{}\Omega)$

with $\Omega=\arg(A_{1}^{\ast}A_{2})$. After averaging over the detection
position, the detection averaged probability is given by

$\displaystyle P_{\nu_{e}}(t)$ $\displaystyle=$ $\displaystyle{1\over
2}+{1\over 2}(|a_{1}|^{2}~{}-~{}|a_{2}|^{2})\cos 2\theta_{N}\cos
2\theta_{0}-|a_{1}a_{2}|\sin 2\theta_{N}\cos
2\theta_{0}\cos(\int^{t}_{t_{r}}\Delta_{N}dt~{}+~{}\omega)$

with $\omega=\arg(a_{1}a_{2})$. The last term shows that the phase of the
neutrino oscillation at the point the neutrino enters resonance can
substantially effect this probability. Therefore, we must also average over
the production position, to obtain the fully averaged probability of detecting
an electron neutrino as

$\langle P_{\nu_{e}}\,\rangle=~{}{1\over 2}~{}+~{}\biggr{(}{1\over
2}~{}-~{}P_{x}\biggr{)}\cos 2\theta_{N}\cos 2\theta_{0}\,.$ (6)

where $P_{x}~{}=~{}|a_{2}|^{2}$, the probability of transition from
$|\nu_{2},N\rangle$ to $|\nu_{1},N\rangle$ (or vice versa) during resonance
crossing. The adiabatic case Barger et al. (1986) is trivially obtained by
setting $P_{x}=0$. Also, if the electron neutrinos are produced at a density
much greater than the resonance density, so that $\cos 2\theta_{N}\approx-1$,
then

$\langle P_{\nu_{e}}\,\rangle\approx~{}\sin^{2}\theta_{0}~{}+~{}P_{x}\,\cos
2\theta_{0}\,.$ (7)

Thus, in the very small $\theta_{0}$ limit, the survival probability is just
equal to the probability of level crossing during resonance passage.

Similar calculations can also be performed for the case of double resonance
crossing (neutrinos from the far side of the Sun). Here we must average not
only over the production and detection positions of the neutrino but also over
the separation between resonances. This sensitivity to the separation of the
resonances can be understood as the effect of the phase of the oscillation as
the neutrino enters the second resonance region. The fully averaged
probability of detecting an electron neutrino is the same as Eq.(6) with
$P_{x}$ replaced by $P_{1x}(1-P_{2x})~{}+~{}(1-P_{1x})P_{2x}$ (the classical
probability result). Therefore, the generalization to any number of resonance
regions, suitable averaged, is obvious.

To calculate the probability, $P_{x}$, I make the approximation that the
density of electrons varies linearly in the transition region. That is, a
Taylor series expansion is made about the resonance position and the second
and higher derivative terms are discarded;

$N(t)~{}\approx~{}N(t_{r})~{}+~{}(t-t_{r}){dN\over dt}|_{t_{r}}.$ (8)

In this approximation the probability of transition between adiabatic states
was calculated by LandauLandau (1932) and ZenerZener (1932). This is achieved
by solving the Schrodinger equation, Eq.(2), exactly in this limit. The
solution is in terms of Weber (parabolic cylinder) functions. Application of
the Landau-Zenner result to the current situation gives

$P_{x}=\exp\left[-\frac{\pi}{2}\frac{\sin^{2}2\theta_{0}}{\cos
2\theta_{0}}\,\frac{\Delta_{0}}{|(1/N)dN/dt|_{t_{r}}}\right]\,.$ (9)

This expression, together with Eq.(6) are the main analytical results of this
paper and demonstrate that only the electron number density, at production,
and the logarithmic derivative of this density, at resonance, determine the
probability of detecting an electron neutrino in the vacuum. It should be
emphasized here, that this result assumes that the neutrino state is produced
before significant transitions take place and thus Eq.(9) is not valid for
neutrinos produced in the transition region.

From Eq.(9) the size of the transition region can be determined. There are
significant transitions ($P_{x}~{}>~{}0.01$) if
$\theta_{0}~{}<~{}\theta_{crit}$ where $\theta_{crit}$ satisfies

${\sin^{2}2\theta_{crit}\over\cos
2\theta_{crit}}~{}=~{}3~{}{1\over\Delta_{0}}~{}\biggr{|}\frac{1}{N}\frac{dN}{dt}\biggr{|}_{t_{r}}.$
(10)

Hence, the maximum separation between the eigenstates for which transitions
take place is $\Delta_{0}\sin 2\theta_{crit}$. Therefore, the transition
region is defined by

$\Delta_{N}~{}<~{}\Delta_{0}\sin 2\theta_{crit}.$ (11)

This can only happen if $\theta_{0}<\theta_{crit}$. In this transition region,
the maximum variation of the electron number density from the resonant value
is $\pm\,\delta N$, where

$\displaystyle\delta N/N(t_{r})=\sin 2\theta_{crit}.$

Thus, the size of the transition region is

$\displaystyle|t-t_{r}|$ $\displaystyle=$ $\displaystyle\sin
2\theta_{crit}\biggr{/}\biggr{|}\frac{1}{N}\frac{dN}{dt}\biggr{|}_{t_{r}}.$

This is the maximum $|t-t_{r}|$ for which the linear approximation must be
good, so that Eq.(9) gives a reasonable estimate of the probability of
crossing. For an exponential density profile, the Taylor series expansion is
an expansion in $\sin 2\theta_{crit}$, so that small $\theta_{crit}$ is an
excellent approximation.

For the sun, the density profile is exponential except for the region near the
center. In Fig. 1, I have plotted the probability contours for detection of an
electron neutrino at the Earth in the $\Delta_{0}/\sqrt{2}G_{F}N_{c}$ vs $\sin
2\theta_{0}$ plane for an exponential density profile. $N_{c}$ is the electron
number density at the point at which the electron neutrinos are produced. This
plot depends only on the properties of the sun and this dependency is only
through the combination $R_{s}N_{c}$ where $R_{s}$ is the scale height. For
Fig. 1, I have used an $N_{c}$, corresponding to a density of 140 g/cm3 and
$Y_{e}$ =0.7. The scale height $R_{s}$ is 0.092 times the radius of the sun.

Figure 1: Probability-contour plot for detecting an electron neutrino at the
Earth which was produced in the solar interior.

Above the line $\Delta_{0}/\sqrt{2}G_{F}N_{C}=1/\cos 2\theta_{0}$, the
neutrinos never cross the resonance density on their way out of the sun. Here,
the probability of detecting an electron neutrino is close to the standard
neutrino oscillation result. Below this line, the effects of passing through
resonance come into play. Inside the 0.1 contour, there is only a small
probability of transitions between the adiabatic states as the neutrino passes
through resonance. To the right-hand side of this contour, the probability of
detecting a neutrino grows, not because of transitions, but because both
adiabatic states have a substantial mixture of electron neutrino at zero
density. To the left-hand side and below this contour, the probability grows
because here there are significant transitions between the adiabatic states as
the neutrino crosses resonance. The diagonal lines of these contours have
slope of -2 because of the form of $P_{x}$. It is only the intercept of these
lines which depends on the product $R_{S}N_{c}$. Therefore, if one wishes to
change the production density, which is held fixed in this plot, only these
lines need to be shifted. In fact, a line labeled with $P_{e}$ “crosses”
$\Delta_{0}/\sqrt{2}G_{F}N_{C}=1$, when a small $\theta_{0}$ satisfies

$\frac{\sin^{2}2\theta_{0}}{\cos
2\theta_{0}}=\frac{-\sqrt{2}\,\text{ln}(P_{e})}{\pi G_{F}R_{s}N_{c}}\,.$ (12)

Note that I find the probability of detecting an electron neutrino, which
crosses resonance, to be greater than 0.25 when $\theta_{0}<0.01$.

This iso-probability plot can easily be converted into an approximate iso-SNU
(solar neutrino units) plot for the Davis et al experimentBahcall et al.
(1985). The predicted result for this experimentBahcall et al. (1982) is 6
SNU, with 4.3 SNU coming from the 8B neutrinos and 1.6 SNU from the lower-
energy neutrinos (pep, 7Be, 13N, and 15O), whereas Davis et al observe
2.1$\pm$0.3 SNU. Roughly speaking, the 2 SNU contour, in the
$m^{2}_{2}-m^{2}_{1}$ vs $\sin 2\theta_{0}$ log-log plot, will be a triangle,
similar to the 0.3 contour of Fig. 1, with rounded corners. The three straight
sections of this triangle are approximately given below. The horizontal line
is given by choosing the parameters so that all the low-energy neutrinos and
only 12% of the 8B neutrinos are observed. This gives the constraints obtained
by BetheBethe (1986),

$\displaystyle(m^{2}_{2}-m^{2}_{1})\approx 8\times 10^{-5}eV^{2}\,,$
$\displaystyle 0.03<\sin 2\theta_{0}<0.6\,.$ (13)

For the vertical line, the probability of detecting an electron neutrino)is
nearly independent of energy, if $1>\Delta_{0}/\sqrt{2}G_{F}N>10^{-3}$.
Therefore, we need to reduce all neutrinos by 30% Barger et al. (1986). This
is achieved when

$\displaystyle 8\times 10^{-8}eV^{2}<(m^{2}_{2}-m^{2}_{1})<1\times
10^{-5}eV^{2}\,,$ $\displaystyle\sin 2\theta_{0}\approx 0.9.$ (14)

For the diagonal line, we need to arrange that the Davis experiment only
observed 50% of the 8B neutrinos and none of the lower-energy neutrinosRosen
and Gelb (1986); Kolb et al. (1986). This is achieved when the probability for
the mean 8B neutrino, weighted by the detector cross section (energy $\sim$10
MeV), is 0.5. This gives the following constraint:

$\displaystyle(m^{2}_{2}-m^{2}_{1})\sin^{2}2\theta_{0}=3\times
10^{-8}eV^{2}\,,$ $\displaystyle 0.03<\sin 2\theta_{0}<0.6.$ (15)

To summarize, Eqs. (13)—(15) give regions of parameter space for which the
expected result from the Davis experiment is $\sim$2 SNU.

Since the proposed gallium experiment observes lower-energy neutrinos, from
the pp process, these three regions will be distinguishable by use of the
results of this experiment. More precise iso-SNU plots, for both experiments,
are being generated taking into account the production energy and production
position distributions of the neutrinos from the various processes within the
solar interior.

I would like to acknowledge discussions with T. Walker and R. Kolb. Fermilab
is operated by the Universities Research Association Inc., under contract with
the United States Department of Energy.

## References

  * Mikheev and Smirnov (1986) S. P. Mikheev and A. Y. Smirnov, Nuovo Cim. C 9, 17 (1986).
  * Bethe (1986) H. A. Bethe, Phys. Rev. Lett. 56, 1305 (1986).
  * Wolfenstein (1978) L. Wolfenstein, Phys. Rev. D 17, 2369 (1978).
  * Wolfenstein (1979) L. Wolfenstein, Phys. Rev. D 20, 2634 (1979).
  * Barger et al. (1986) V. D. Barger, R. J. N. Phillips, and K. Whisnant, Phys. Rev. D 34, 980 (1986).
  * Landau (1932) L. D. Landau, Phys. Z. Sowjetunion 2, 46 (1932).
  * Zener (1932) C. Zener, Proc. Roy. Soc. Lond. A 137, 696 (1932).
  * Bahcall et al. (1985) J. N. Bahcall, B. T. Cleveland, R. Davis, Jr., and J. K. Rowley, Astrophys. J. Lett. 292, L79 (1985).
  * Bahcall et al. (1982) J. N. Bahcall, W. F. Huebner, S. H. Lubow, P. D. Parker, and R. K. Ulrich, Rev. Mod. Phys. 54, 767 (1982).
  * Rosen and Gelb (1986) S. P. Rosen and J. M. Gelb, Phys. Rev. D 34, 969 (1986).
  * Kolb et al. (1986) E. W. Kolb, M. S. Turner, and T. P. Walker, Phys. Lett. B 175, 478 (1986).


# Sounding the metabolic orchestra:
A delay dynamical systems perspective on the glucose-insulin regulatory
response to on-off glucose infusion

Stefan Ruschel  Northumbria University, Department of Mathematics, Physics and
Electrical Engineering, Newcastle upon Tyne, NE1 8ST, United Kingdom  The
University of Auckland, Department of Mathematics, Auckland, 1142, New Zealand
<EMAIL_ADDRESS>Benoit Huard  Northumbria University,
Department of Mathematics, Physics and Electrical Engineering, Newcastle upon
Tyne, NE1 8ST, United Kingdom<EMAIL_ADDRESS>
###### Abstract

We investigate the consequences of periodic, on-off glucose infusion on the
glucose-insulin regulatory system on the basis of a system-level mathematical
model with two explicit time delays. Studying the effects of such infusion
protocols is mathematically challenging yet a promising direction for probing
the system response to infusion. We pay special attention to the interplay of
the infusion with intermediate-time-scale, ultradian oscillations that arise
as a results of the physiological response of glucose uptake and back-release
into the bloodstream. By using numerical solvers and numerical continuation
software, we investigate the response of the model to different infusion
patterns, and explore how these patterns affect the overall levels of glucose
and insulin, and can lead to entrainment. By doing so, we provide a road-map
of system responses that can potentially help identify new test strategies for
detecting abnormal responses to glucose uptake.

## Introduction

Cyclic rhythms are widely recognized for their significant role in regulating
the function of biological and physiological systems [1, 2]. Endogenic
oscillations are typically encountered in healthy individuals, while a
progressive lack of control of these rhythms is often associated with system
stress (e.g. sleep deprivation [3, 4]), and disease evolution in humans [5].

A prominent example of such endocrine oscillations in the human body is the
self-regulation of blood glucose levels [6]. When blood glucose levels
increase, insulin is released from the pancreas. Insulin then causes blood
glucose levels to decrease by stimulating body cells to absorb glucose from
the blood. Conversely, when blood glucose levels fall, pancreatic
$\beta$-cells release glucagon stimulating hepatic glycogenolysis and
neoglucogenesis. The level of blood glucose is then controlled by the rates of
insulin secretion (activation by glucose) and hepatic glucose production
(inhibition by insulin). Within the glucose-insulin regulatory system, both
rapid oscillations of insulin (period $\sim$ 6-15 minutes), and ultradian
oscillation of glucose and insulin (of similar period $\sim$ 80-180 minutes
[7]) have been observed during fasting, meal ingestion, continuous enteral and
intravenous nutrition [8].

The most important pathway to understand the underlying mechanisms of these
glucose-insulin oscillations is measuring the response to glucose infusions. A
large quantity of metrics and mathematical models have been devised for that
purpose. While the HBA1c metric remains an essential tool for the diagnostic,
prevention and control for T2 diabetes [9], clinical tests involving patterns
of glucose intake combined with mathematical models provide a mechanism for
evaluating the efficacy of internal regulation [10, 11, 12, 13]. The minimal
model devised by Bergman and Cobelli [14, 15] provides an effective method for
estimating insulin sensitivity from an intravenous or oral glucose tolerance
test, although it can lead to underestimation in individuals with a large
acute insulin response [16]. With the wider availability of continuous glucose
monitors and automated insulin pumps, the ability to detect diabetic
deficiencies relies on the capacity of models to reproduce more complex and
realistic dynamics under various routine life conditions such as, for example,
sleep deprivation [3].

The main goal of this article is to identify the types of behaviors in a
suitable mathematical model that can be expected as a response to periodic
glucose uptake, specifically periodic on-off glucose infusion, which can be
readily implemented in practice. We focus on the capacity of the system to
fall into lockstep with the frequency of the glucose stimulus (so-called
entrainment) which has been observed in numerous contexts at the ultradian and
circadian levels in endocrinology [17, 18, 13], but especially in models of
glucose-insulin oscillations with periodic infusion [19].

Many modeling efforts have been made to replicate the nonlinear response of
the glycolytic system; in particular, the mathematical modeling of the delayed
response of individual parts of the system by explicit time delays has proven
an effective means to explain the onset of self-sustained, ultradian
oscillations in the glycemic system [20, 21, 22]. A common approach to
modeling oscillatory behavior of complex biology is to consider time delays
[23]. In particular, models of endocrine regulation often incorporate explicit
delays to account for the time required for the synthesis, release, and action
of hormones or metabolites [17]. Various models of intrapancreatic rhythmic
activity have been proposed recently, see Ref. 13 for a review. For example,
it was shown that glucose oscillations can enhance the insulin secretory
response at the $\beta$-cell level when tweaked at a suitable amplitude and
frequency [24]. Negative delayed feedback has also been shown to provide a
suitable explanatory mechanism for the coordinated pancreatic islet activity
[25].

Figure 1:  Panel (a): Diagrammatic overview of the glucose-insulin regulatory
delayed-feedback model (1)–(2); see methods section for details. Panels
(b)–(e): Characteristic time series of system (1)–(2) (with positive initial
condition) for different patterns of glucose infusion with intervals of
fasting indicated by a white background and intervals of glucose infusion with
constant rate indicated by a light blue background. Units are [$G$] mg dl-1,
[$I$] mg dl-1, and [$t$] h. Infusion rates when not fasting are
$G_{\text{in}}=1.35$ mg dl-1min-1 in panels (c)–(d) and $G_{\text{in}}=24.3$
mg dl-1min-1 in panel (e); period of infusion is $T_{\text{in}}=1$ h in panel
(d) and $T_{\text{in}}=3$ h in panel (e); time of infusion is
$t_{\text{in}}=30$ min in panel (d) and $t_{\text{in}}=5$ min in panel (e).

In this article, we investigate a two-component, system-level mathematical
model (see Eqs. (1)–(2) in the methods section) for blood glucose level $G(t)$
and insulin levels $I(t)$ with two explicit time delays $\tau_{I}$ and
$\tau_{G}$ corresponding to pancreatic insulin and hepatic glucose production
pathways, see the methods section for details on the model. The model
incorporates the following physiological processes and factors that influence
glucose and insulin dynamics, see Fig. 1(a) for a schematic overview.

  * •

Glucose uptake: $G_{\text{in}}$ represents glucose uptake into the blood by
meal ingestion, continuous enteral or intravenous nutrition.

  * •

Insulin production: $f_{1}$ represents the production of insulin. It is
influenced by the concentration of glucose with a delay $\tau_{I}$ to account
for the time lag between high glucose levels triggering insulin production in
the pancreas and when it becomes available for reducing glucose in the
bloodstream.

  * •

Insulin-independent glucose utilization: $f_{2}$ describes the utilization of
glucose by tissues, mainly the brain, in an insulin-independent manner. It
does not rely on the presence of insulin.

  * •

Insulin-dependent glucose utilization: $f_{3}\cdot f_{4}$ represents the
utilization of glucose by muscle tissues in an insulin-dependent manner. It
reflects the capacity of tissues to utilize insulin for glucose uptake.

  * •

Glucose production by the liver: $f_{5}$ represents the production of glucose
by the liver. The delay $\tau_{G}$ represents the time between hepatic glucose
production and insulin stimulation.

  * •

Insulin degradation: The rate $d$ accounts for the degradation of insulin in
the body, primarily by the liver and kidneys. It combines both natural factors
(e.g., exercise) and artificial factors (e.g., medication) that influence the
rate of insulin degradation.

The nonlinear pathways $f_{1},f_{2},f_{3},f_{4},$ and $f_{5}$ are represented
using Hill functions, which are mathematical functions commonly used in
biological modeling. These functions introduce additional parameters that have
specific physiological interpretations and allow for a more accurate
representation of the underlying dynamics of the glucose-insulin system, see
methods section for details. The delays $\tau_{I}$ and $\tau_{G}$ are
important physiological parameters encapsulating the responsitivity of the
signaling and production pathways. They are assumed to be constant for the
purpose of this article, although in practice, they can vary between
individuals, as well as during the day and lifespan, and especially in the
presence of diabetes.

The model has been extensively analyzed by various authors in the case of
constant rates of glucose infusion [20, 26, 27, 28]. It originates from the
work of Sturis and collaborators who devised a model of glucose and insulin
ultradian oscillations which were observed experimentally under various
conditions [29]. We also remark here that the model belongs to a larger class
of models incorporating delays to capture secretion processes [11, 30]. We
extend these earlier efforts on the analysis of the model by studying its
response to periodic variations of the parameter $G_{\rm{in}}$, that is,
periodic variations of glucose uptake. In particular, we consider on-off
infusion, a form of periodic infusion that is comparatively easy to implement
in practice, where the rate of glucose infusion periodically switches between
a positive constant value and zero. Panels (b)–(e) of Fig. 1 show prototypical
examples for the response of system (1)–(2) for various types of glucose
uptake during fasting (b), glucose infusion with a (relatively high) constant
rate (c) and periodic on-off infusion (d)–(e). We first investigate the loss
of ultradian oscillations under sufficiently strong constant infusion, see
Fig. 1(b)–(c). We then aim to study the effects of different glucose infusion
patterns $G_{\rm{in}}$ on glucose homeostasis, in particular the transiton
from quasi-periodicity to entrainment, as shown in Figs. 1(d)–(e).

## Results

### Ultradian oscillations

It has been shown that for a fixed constant glucose infusion $G_{\rm{in}}$,
sufficiently large values of the response delays $\tau_{1}$ and $\tau_{2}$
lead to periodic oscillations in system (1)–(2) with periods closely
resembling the observed range for ultradian oscillations[20, 27, 28].
Mathematically speaking, the onset of oscillations is mediated by a
supercritical Hopf bifurcation that leads to a local topological change in the
solution space of system (1)–(2) from a stable equilibrium to a situation of
an unstable equilibrium surrounded by a small stable limit cycle close to the
bifurcation point[26]. For details on bifurcation theory and the Hopf
bifurcation, we refer the interested reader to Ref. 31. To witness the
bifurcation point, it necessary to vary at least one parameter of the system.
Here, we focus on the response delays $\tau_{I}$ and $\tau_{G}$. Allowing
these two parameters values to vary simultaneously, one obtains a one-
parameter curve $\mathbf{H(\omega)}=(\tau_{I}(\omega),\tau_{G}(\omega))$ of
Hopf bifurcation in the $(\tau_{I},\tau_{G})$-plane in terms of the Hopf
frequencies $\omega$, see Methods section for a detailed derivation. The curve
$\mathbf{H(\omega)}$ corresponds to the critical curve for oscillations in
system (1)–(2).

Figure 2:  Characterization of fasting oscillations with respect to response
delays. Panels show the period (a), maximum glucose value (b), and minimum
glucose value (c) as a function of response delays $\tau_{I}$ (min) and
$\tau_{G}$ (min). Shown are the critical curve for oscillations (black, Hopf
bifurcation), and iso-curves (blue) with constant period (a), glucose-maxima
(b) and minimum of $G$ (c). The light blue rectangle shows the physiological
range of delay values for comparison. See methods section for the model and
choice of parameters.

Figure 2(a) shows the curve $\mathbf{H}$ (black) during fasting, i.e. for
$G_{\text{in}}=0$, computed with the software package DDE-Biftool for Matlab
[32, 33]. It has been numerically verified that the curve $\mathbf{H}$ is
indeed supercritical for the range of parameter values considered. Figure 2(a)
can be interpreted as follows: First, for value pairs above the curve and for
$\tau_{I}\leq 20\leavevmode\nobreak\ \mathrm{min},$ $\tau_{G}\leq
60\leavevmode\nobreak\ \mathrm{min}$, any solution of the model starting in a
physiological range of glucose and insulin develops periodic oscillations, see
Fig. 1(b). Second, for value pairs $(\tau_{I},\tau_{G})$ below the curve
$\mathbf{H}$, oscillations in system (1)–(2) decay and approach the
equilibrium $(G^{\ast},I^{\ast}).$ This possibly reflects the situation where
an individual is administered a glucose dose that is too high to be managed in
an oscillatory manner within physiological glucose and insulin ranges, compare
Fig. 1(c).

Figure 2(a) also gives an overview of the resulting period of oscillation
above the critical curve $\mathbf{H}$ shown in the form of isocurves (blue) of
limit cycles with constant period. The physiological range of parameters
$(\tau_{I},\tau_{G})$ is highlighted by light blue square in the background
for convenience. The range of expected periods for ultradian oscillations as
predicted by the model thus ranges from $2.2$ to $4.2$ hours during fasting.
More generally, we observe the period of the limit cycle oscillation grow
approximately linear with the sum of the two delay values $\tau_{I}+\tau_{G}$.
We also observe that away from the curve $\mathbf{H}$, the limit cycle
oscillation becomes less and less sinusoidal, i.e. the nonlinearity of system
(1)–(2) has more and more of an effect on the limit cycle. Panels (b)–(c) of
Fig. 2 illustrate this effect by plotting isocurves of periodic orbit with
constant minimum and maximum glucose within one period of oscillation. We
observe that, whereas the glucose minimum decreases approximately linearly
with the sum of the delays $\tau_{I}+\tau_{G}$, the maximum $G$ remains almost
constant for the range of parameter values considered. Note that this
predicted effect of long response delays is potentially harmful and is
virtually undetectable by common testing methods.

On the other hand, we observe that, for fixed values of the delays, gradually
increasing the glucose infusion leads to a loss of oscillations. This
phenomenon has been observed before and can be interpreted as an insufficient
insulin secretion to accommodate the infusion, forcing the system to lower
glucose lower levels[26]. Figure 3 shows how the location of the curve
$\mathbf{H}$ changes for various levels of constant glucose infusion. We
observe a two different types of change for values in the approximate ranges
$0\leq G_{\text{in}}\leq 0.55$ mg dl-1 min-1 and $G_{\text{in}}>0.55$ mg dl-1
min-1 shown in panels (a) and (b) of Fig. 3, respectively. Figure 3(a)
suggests that low levels of $G_{\text{in}}$ promote oscillations in system
(1)–(2) as compared to the fasting case. This trend reverses at approximately
at $G_{\text{in}}=0.55$ mg dl-1 min-1, where the location of the curve
$\mathbf{H}$ starts moving to larger and larger values of $\tau_{I}$ and
$\tau_{G}$, see Fig. 3(b). Approximately at $G_{\text{in}}=1.2$ mg dl-1 min-1
the position of $\mathbf{H}$ is comparable with the starting location for
$G_{\text{in}}=0$. Further increasing $G_{\text{in}}$ moves $\mathbf{H}$
inside the physiological range of delay values (light blue) and finally beyond
causing all oscillations to cease in the physiological parameter regime.
Compare also Fig. 1(b)-(c) for an illustration of this transition and the loss
of oscillations for $(\tau_{I},\tau_{G})=(5,20)$.

Figure 3:  Position of the critical curve (curve of Hopf bifurcation) in the
$(\tau_{I},\tau_{G})$-plane for various values of $G_{\text{in}}$ ranging from
$0$ to $0.5$ and from $0.6$ to $1.6$ mg dl-1 min-1 (all black). The light blue
rectangle shows the physiological range of delay values for comparison.

### Entrainment and amplitude response to on-off glucose infusion

We now investigate the effect of periodic glucose infusion on baseline fasting
oscillations shown in Fig. 1(b), i.e. we fix $\tau_{I}=5$ min, $\tau_{G}=20$
min and periodically adjust the level of $G_{\text{in}}$ between $0$ and a
positive value $G_{0}$ to be specified. The natural frequency of ultradian
oscillation in this case is $T\approx 2.2$ h. We show that the resulting
glucose and insulin ranges depend sensitively on the period of the on-off
infusion. Figures 1(d)–(e) show two of the possible outcomes with different
maximal infusion strength $G_{\max},$ period of infusion $T_{\text{in}}$, and
infusion duration $t_{\text{in}}.$

Figure 4: Response of model (1)–(2) to glucose infusion protocol (13) with
maximum infusion rate $G_{\max}$ (mg/(dl min)) and length of infusion
$t_{\text{in}}=T_{\text{in}}/2$ (h). Shown is the maximum value of $G$ (mg/dl)
in colorcode (blue-white) obtained by integration for various $T_{\text{in}}$
and $G_{\max}$ over $100(\tau_{I}+\tau_{G}+T_{\text{in}})$ time units. The
maximum data is overlaid by curves of torus bifurcation (purple), curves of
fold bifurcation of periodic orbits (red) and curves of period doubling
bifurcation (magenta) bounding regions of locking to the infusion protocol.
Other parameters are $\tau_{I}=5$ min and $\tau_{G}=20$ min.

#### Long infusion time compared to period

Figure 1(d) shows the result of periodic infusion with $G_{\text{in}}=1.35$ mg
dl-1 min${}^{-1},$ for $t_{\text{in}}=30$ min every $T_{\text{in}}=60$ min,
resulting in so-called quasi-periodic oscillations. Indeed, quasi-periodic
oscillations are characterized by the presence of an oscillating envelope of
the oscillation that evolves on a much slower time-scale, compare Fig. 1(d).
This is in sharp contrast with panels (b) (no infusion) and (c) (constant
infusion with the same maximum rate) of Fig. 1, where we have either periodic
oscillations, or a decay of oscillations towards the equilibrium state. Quasi-
periodic oscillations can be expected to occur in oscillatory systems which
are externally driven by an input with non-commensurable period, here
$T_{\text{in}}/T_{0}=2.2$. In this case, the effect of infusion very much
depends on its current state: When insulin is low, glucose increases quickly;
when insulin is high, glucose cannot increase further and the infusion only
delays the expected decrease in glucose levels.

Periodicity of the oscillations can be restored by adjusting $G_{\text{in}}$
and $T_{\text{in}}$. Figure 4 summarizes the response of system (1)–(2) to
periodic forcing with different values of $G_{\text{in}}$ and $T_{\text{in}}$.
The locus in parameter space of the quasi-periodic oscillation shown in Fig.
1(d) is indicated by a green rectangle. Figure 4 shows the overall glucose
maximum (in color code) observed over a time span of
$100\cdot(T_{\text{in}}+\tau_{I}+\tau_{G})$ minutes. The various mechanisms
generating periodic rhythms can be understood from the numerically computed
bifurcation curves shown in 4. These correspond to curves of torus
bifurcations $\mathbf{T}$ (purple), curves $\mathbf{F}$ (red) of fold (or
saddle-node) bifurcations of periodic orbits, and curves $\mathbf{PD}$
(magenta) of period-doubling bifurcations of periodic orbits. These mark the
transition to periodic solutions and thus characterize the so-called
_entrainment_ of oscillations.

The curves $\mathbf{F}$ respectively enclose deltoid-like regions – called
resonance or locking tongues – extending from the line $G_{\text{in}}=0$,
inside of which we observe periodic oscillations. The curves $\mathbf{F}$
emerge pairwise from resonant points where the infusion period is a rational
multiple of the natural period of the system without infusion, i.e.
$pT_{\text{in}}=qT_{0}$ for integers $p,q$. Figure 4 shows the first three
principal resonances of system (1)–(2) where $p=1,2,3$ and $q=1$. It is
expected that such resonance tongues emanate from the line $G_{\max}$ at every
rational point $T_{0}$. These higher order resonances (exept $p=4$ and $q=1$
which is outside of the considered range of parameter values) has been
omitted/not computed as they are are typically very narrow and thus unlikely
to be physiologically relevant.

This behavior persists moving towards larger values of $G_{\max}$ into the
regions that are bounded approximately by the curves $T$, where the underlying
stable periodic orbit destabilizes and gives rise to a torus that corresponds
to quasi-periodic oscillations. We find numerical evidence that the direction
with which this torus emanates from the the curve $\mathbf{T}$ can change and
gives rise to the discontinuous transition between the observed maximum values
in Fig. 4. The curves $\mathbf{T}$ each emanate from either point of
intersection with a curve $\mathbf{F}$ or $\mathbf{PD}$. Intersections with
curves $\mathbf{PD}$ correspond to higher order locking between the ultradian
oscillations and the infusion. Overall, we observe that the strength and
period of the infusion have a crucial effect on the resulting amplitude of the
oscillations. For instance, forcing the system periodically with
$T_{0}=T_{\text{in}}$ and relative amplitude $G_{\text{in}}=1$ mg dl-1 min-1
leads to a $40\%$ increase of the overall amplitude of the oscillation (which
appears to be still in physiological range). In contrast, stimulating the
system with a gradually increasing $G_{\text{in}}$ in the 2:1 regime first
goes through phase during which glucose amplitudes remain relatively constant
before slowly increasing.

More generally, we observe that, for the assumed values of the response
delays, periodic infusion with $T_{\text{in}}=2t_{\text{in}}>T_{0}$ and
$G_{\text{in}}$ is sufficient for the resulting period of the resulting
glucose-insulin oscillation to be set by (locked to) the period of glucose
infusion.

#### Short infusion time compared to period

We note here that locking can be achieved when the same glucose dose is
delivered in a shorter period of time, resulting in a more concentrated and
intense infusion. To further explore this phenomenon, we conducted additional
experiments using an on-off glucose infusion protocol with a fixed infusion
period of $T_{\rm{in}}=180$ min. Figure 5 showcases the results obtained from
these experiments, where we varied both the infusion time $t_{\rm{in}}$ and
the average glucose dose per minute $\bar{G}$, represented by $G_{\max}\cdot
t_{\rm{in}}/T_{\rm{in}}$.

In this figure, we observe a locus in the parameter space that corresponds to
the quasi-periodic orbit illustrated in Fig. 1(e). This locus is denoted by a
distinctive yellow diamond marker, which highlights the specific combination
of infusion time and glucose dose that leads to the observed quasi-periodic
behavior. Additionally, we present a curve labeled as $\mathbf{T}$, which
represents a torus bifurcation curve. This curve serves as an indicator of the
critical transition point between entrainment and quasi-periodic oscillation
in response to the infusion protocol.

Figure 5:  Response of model (1)–(2) to glucose infusion protocol (13) with
average infusion rate $\bar{G}=G_{\max}\cdot t_{\text{in}}/T_{\text{in}}$ mg
dl-1 min-1 over the length of infusion $t_{\text{in}}$ (min) with constant
period $T_{\text{in}}=180$ (min). Shown is the maximum value of $G$ (mg/dl) in
colorcode (blue-white) obtained by integration for various $t_{\text{in}}$ and
$\bar{G}$. The maximum data is overlaid by curves of torus bifurcation
(purple), curves of fold bifurcation of periodic orbits (red) and curves of
period doubling bifurcation (purple) bounding regions of locking to the
infusion protocol. Other parameters are $\tau_{I}=5$ min and $\tau_{G}=20$
min.

## Discussion

It is well documented that glucose rhythms stimulate pulsatile pancreatic
insulin secretion at various timescales [29, 34]. For example, the 1:1
entrainment mode – namely one ultradian glucose oscillation per glucose
infusion cycle – was clinically shown to be present using a sinusoidal glucose
infusion in individuals without diabetes [8, 19]. Our analysis of periodically
driven ultradian oscillations highlights that a periodic on-off stimulus,
closer to normal daily conditions, also possesses the ability to entrain
glucose rhythms. Furthermore, the duration of each glucose input has a crucial
impact on the generation of periodic rhythms, as well as on attained glycemic
levels. This theoretically provides a method for delivering a fixed glucose
dose while minimising the amplitude of the resulting rhythm. This can be
achieved by either altering the period of the infusion, or the length of each
pulse. This is most observable in figure 5, where stretching the infusion
duration leads to lower glucose amplitudes. For example, consider a scenario
where glucose is infused every 180 minute over a 12-hour period. Infusing a
dose with $G_{\max}=2.4$ mg dl-1 min-1 over $t_{in}=30$ minutes leads to a
maximal glucose value around $150$ mg/dl. In contrast, a dose with
$G_{\max}=1.2$ mg dl-1 min-1 over $t_{in}=60$ minutes reduces the maximal
glucose level to around $125$ mg/dl. In both cases, the average dose per
minute is $\bar{G}=0.4$ mg dl-1 min-1, and a total dose of $288$ mg/dl is
infused over the 12 hour timespan.

Our study provides valuable insights into the system’s response to glucose
infusion patterns, providing multiple pathways for the production of stable
oscillatory rhythms and a similar entrainment structure is also expected for
simpler models of glucose-insulin regulation featuring delays, e.g. [35, 36,
30]. Nonetheless, there are several limitations that should be considered.
First, let us note that while the exact location of bifurcation curves would
depend on model parameters, the bifurcation types are likely to remain the
same for parameter ranges representing non-diabetic individuals. Our model
assumes fixed values for the delays in insulin and glucose production
pathways, represented by $\tau_{I}$ and $\tau_{G}$ respectively. In reality,
these delays can vary between individuals and change over short and long
timescales due to daily-life factors such as exercise, aging and the presence
of insulin resistance. Future research could incorporate individual-specific
delays to account for this variability and investigate their impact on the
system’s dynamics.

It is worth noting that our model relies solely on plasma glucose and insulin
measurements for prediction, which highlights the importance of accurate and
reliable measurements in clinical settings. The nonlinear structure of the
model allows for the description of nontrivial dynamics and enhances parameter
identifiability. This aspect is crucial for developing robust and accurate
models that can capture the complex dynamics of the glucose-insulin regulatory
system.

It is also worth noting that the timing of the glucose infusion does not
influence the bifurcation structure (4), nor the glucose-insulin ranges of the
periodic rhythms. In other words, the long term dynamics is not dependent on
the starting time of the periodic on-off glucose infusion. This does not mean
that the timing of glucose inputs bears little importance. While the
investigated infusion ranges ensured the positivity of glucose and insulin
values, values below or above healthy physiological ranges may appear in the
transient path to the limit cycle. In turn, additional dynamics may emerge
from interactions with other physiological feedback loops or subsystems, such
as the glucagon pathway or the hypothalamic-pituitary-adrenal axis, for which
the alignment with glucose regulation is essential for maintaining good health
[37]. The recent incorporation of glucagon [22] in models of the glucose-
insulin feedback system may help provide a more complete and quantitative
picture of dynamical interactions occurring within the pancreas [38, 39] which
can be used to improve quantitative tests for the detection and measurement of
insulin and glucagon resistance [40].

Another aspect to consider is the interaction between the glucose-insulin
regulatory system and other physiological processes. Our model focuses solely
on the glucose-insulin loop, but in reality, there are complex interactions
between various metabolic pathways, hormones, and organs. Integrating these
interactions into a comprehensive model could provide a more complete
understanding of the system’s behavior and its response to different stimuli.

## Conclusion

In this study, we employed a system-level mathematical model to investigate
the response of the glucose-insulin regulatory system to periodic glucose
infusion. By exploring different glucose infusion patterns and analyzing the
resulting dynamics, we gained insights into the system’s behavior and
identified key factors influencing its response.

Our findings demonstrate that the glucose-insulin regulatory system exhibits a
range of behaviors depending on the glucose infusion pattern. When a constant
glucose infusion is applied, the system shows ultradian oscillations
characterized by periodic variations in glucose and insulin levels. However,
as the glucose infusion rate exceeds a certain threshold, these oscillations
disappear, and the system focuses on reducing glucose levels without
exhibiting oscillatory behavior. This observation suggests a physiological
limit beyond which the system’s oscillatory capacity is overwhelmed.

We further investigated the effects of periodic on-off pulses, mimicking
repeated intravenous glucose tolerance tests. Our analysis revealed that the
period of the on-off pulses plays a crucial role in determining the system’s
dynamics and glucose-insulin ranges. Different patterns of oscillations,
including stable limit cycles and irregular oscillations, were observed for
varying infusion periods. This highlights the importance of considering the
frequency and duration of glucose stimuli in understanding the system’s
response.

The results of this study have important implications for understanding
glucose regulation in both normal and abnormal physiological conditions. By
elucidating the system’s response to different glucose infusion patterns, our
findings can inform the development of test strategies for evaluating the
system’s performance and identifying potential dysfunctions. Furthermore, they
provide insights into the underlying mechanisms governing glucose-insulin
dynamics, contributing to the broader understanding of metabolic regulation.

In conclusion, our study enhances our understanding of the glucose-insulin
regulatory system by investigating its response to periodic glucose infusion.
We identified the impact of different glucose infusion patterns on the
system’s dynamics and demonstrated the importance of various types of glucose
stimuli. These insights can aid in the development of diagnostic and
therapeutic strategies for glucose regulation and contribute to advancements
in the management of metabolic disorders. Future research should aim to
incorporate individual-specific delays and consider the broader physiological
context to further refine our understanding of glucose regulation and its
implications for human health.

## Methods

### The glucose-insulin regulatory delayed-feedback model

We consider the system-level mathematical model

$\displaystyle G^{\prime}(t)$
$\displaystyle=G_{\text{in}}(t)-f_{2}(G(t))-f_{3}(G(t))f_{4}(I(t))+f_{5}(I(t-\tau_{G}))$
(1) $\displaystyle I^{\prime}(t)$
$\displaystyle=I_{\text{in}}(t)+f_{1}(G(t-\tau_{I}))-dI(t)$ (2)

with variables $I(t)$ and $G(t)$ representing the concentrations and of
glucose (mg dl-1) and insulin (uU ml-1) in the plasma at time instant $t$.
System (1)–(2) explicitly depends on time delays $\tau_{I}$ and $\tau_{G}$
respectively representing the system’s response time to insulin production as
a result to glucose uptake, and the production of glucose by the liver as a
result of low insulin levels. Glucose intake and insulin infusion are modeled
by parameters, here called $G_{\rm{in}}$ and $I_{\rm{in}}$. The physiological
response of body is modeled by the nonlinearities

$\displaystyle f_{1}(G)$
$\displaystyle=\frac{R_{m}G^{h_{1}}}{G^{h_{1}}+(V_{g}k_{1})^{h_{1}}},$
$\displaystyle f_{2}(G)$
$\displaystyle=\frac{U_{b}G^{h_{2}}}{G^{h_{2}}+(V_{g}k_{2})^{h_{2}}},$
$\displaystyle f_{3}(G)$ $\displaystyle=\frac{G}{C_{3}V_{g}},$ $\displaystyle
f_{4}(I)$
$\displaystyle=U_{0}+\frac{(U_{m}-U_{0})I^{h_{4}}}{I^{h_{4}}+(1/V_{i}+1/(Et_{i}))^{-h_{4}}k_{4}^{h_{4}}},$
$\displaystyle f_{5}(I)$
$\displaystyle=\frac{R_{g}I^{h_{5}}}{I^{h_{5}}+(V_{p}k_{5})^{h_{5}}},$

where $R_{m}=210,$ $V_{i}=11$, $V_{g}=10$, $E=0.2$, $U_{b}=72$, $t_{i}=100$,
$C_{3}=1000$, $R_{g}=180$, $U_{0}=40$, $V_{p}=3$, $U_{m}=940$, $h_{1}=2$,
$k_{1}=6000$, $h_{2}=1.8$, $k_{2}=103.5$, $h_{4}=1.5$, $k_{4}=80$,
$h_{5}=-8.54$, and $k_{5}=26.7$ with corresponding units. Insulin degradation
is modeled by a constant rate $d$. Throughout the paper we fix $d=0.06$. The
model has been considered before and has been analyzed extensively by various
authors [20, 26, 27, 28]. In particular, it can be shown that, for the
parameter values considered and in the absence of infusion, there is a unique
equilibrium solution $(G^{\ast},I^{\ast})$; see for example [41]. The delay
parameters used for numerical simulation are $\tau_{I}=5$ and $\tau_{G}=20$ if
not stated otherwise. For the general theory of delay differential equations,
such as existence, uniqueness and the stability of solutions, we refer the
interested reader to classic textbooks on the topic[42, 43].

### Critical delay values for oscillatory behavior when infusion rate is
constant

The critical curve for oscillations in the ($\tau_{I},\tau_{G}$)-parameter
plane can be computed from the linearization of system (1)–(2) about the
equilibrium solution $(G^{\ast},I^{\ast})$ and imposing the condition
$\lambda=i\omega,\leavevmode\nobreak\ \omega>0$ (Hopf bifurcation) on
solutions of the corresponding characteristic equation

$0=\chi(\lambda):=\lambda^{2}+\alpha_{1}\lambda+\alpha_{0}+\beta_{1}e^{-\lambda\tau_{1}}+\beta_{2}e^{-\lambda\tau_{2}},$
(3)

where $\tau_{1}=\tau_{I},$ $\tau_{2}=\tau_{I}+\tau_{G}$ and
$\alpha_{1}=f_{2}^{\prime}(G^{\ast})+f^{\prime}_{3}(G^{\ast})f_{4}(I^{\ast})+d$,
$\alpha_{0}=d(f_{2}^{\prime}(G^{\ast})+f^{\prime}_{3}(G^{\ast})f_{4}(I^{\ast}))$,
$\beta_{1}=f^{\prime}_{1}(G^{\ast})f_{3}(G^{\ast})f^{\prime}_{4}(I^{\ast})$,
$\beta_{2}=-f^{\prime}_{1}(G^{\ast})f^{\prime}_{5}(I^{\ast})$. A detailed
derivation of Eq. (3) can be found in [28].
The equation $0=\chi(i\omega)$ can be solved parametrically for $\tau_{1}$ and
$\tau_{2}$ to give

$\displaystyle\tau_{1,2}(\omega)$
$\displaystyle=\frac{1}{\omega}\left(\arctan\left(\frac{\alpha_{1}\omega}{\omega^{2}-\alpha_{0}}\right)+\arccos\left(\frac{\beta_{2,1}^{2}-\beta_{1,2}^{2}-(\omega^{2}-\alpha_{0})^{2}-\alpha_{1}^{2}\omega^{2}}{2\beta_{1,2}\sqrt{(\omega^{2}-\alpha_{0})^{2}+\alpha_{1}^{2}\omega^{2}}}\right)\right),$
(4)

revealing the critical curve for oscillations
$\mathbf{H}\subset\mathbb{R}^{2}$ (curve of Hopf bifurcation)

$\displaystyle\mathbf{H}(\omega)$
$\displaystyle=(\tau_{I}(\omega),\tau_{G}(\omega))=(\tau_{1}(\omega),\tau_{2}(\omega)-\tau_{1}(\omega)).$
(5)

For the considered parameter values, we have that $\alpha_{1}>\alpha_{0}$ and
$\beta_{2}>\alpha_{0}$, ensuring the existence of $\mathbf{H}$. Indeed the
curve is a sharp threshold for oscillation, as it can been shown numerically
that for positive values ($\tau_{I},\tau_{G}$) below $\mathbf{H}$ the fixed
point $(G^{\ast},I^{\ast})$ is stable for any physiological range of starting
values $G$ and $I$. It is worth noting here that system (1)–(2) undergoes
further Hopf bifurcations, respectively at
$\tau_{I,k}(\omega)=\tau_{I}(\omega)+2\pi k/\omega$ and
$\tau_{G,l}(\omega)=\tau_{G}(\omega)+2\pi l/\omega$ with $k,l$ an integer;
however, for the parameter values considered, we can restrict ourselves to the
smallest positive such value pair to cover the physiological parameter range.
The range of relevant values of $\omega$ resulting in positive delays cannot
be computed explicitly, however, straightforward calculations show that the
boundaries $\omega_{I},\omega_{G}$ satisfying $\tau_{I}(\omega_{I})=0$ and
$\tau_{G}(\omega_{G})=0$ are given by

$\displaystyle\omega_{G}$
$\displaystyle=\sqrt{\alpha_{0}-\frac{\alpha_{1}^{2}}{2}+\sqrt{\left(\alpha_{0}-\frac{\alpha_{1}^{2}}{2}\right)^{2}+(\beta_{1}+\beta_{2})^{2}-\alpha_{0}^{2}}},$
(6) $\displaystyle\omega_{I}$
$\displaystyle=\sqrt{\alpha_{0}+\beta_{1}-\frac{\alpha_{1}^{2}}{2}+\sqrt{\left(\alpha_{0}+\beta_{1}-\frac{\alpha_{1}^{2}}{2}\right)^{2}+\beta_{2}^{2}-\alpha_{0}^{2}}},$
(7)

with the corresponding delay values

$\displaystyle\tau_{I}(\omega_{G})$
$\displaystyle=\frac{1}{\omega_{G}}\arctan\left(\frac{\alpha_{1}\omega_{G}}{\omega_{G}^{2}-\alpha_{0}}\right)+\frac{2\pi
k^{\ast}}{\omega_{G}},$ (8) $\displaystyle\tau_{G}(\omega_{I})$
$\displaystyle=\frac{1}{\omega_{I}}\arctan\left(\frac{\alpha_{1}\omega_{I}}{\omega_{I}^{2}-\alpha_{0}-\beta_{1}}\right)+\frac{2\pi
l^{\ast}}{\omega_{I}},$ (9)

where $k^{\ast},l^{\ast}$ are the smallest integers such that $\tau_{G}$ and
$\tau_{I}$ are positive.

We remark that the curve $\mathbf{H}$ vaguely resembles a straight line with
slope $-1$ in the $(\tau_{I},\tau_{G})$-plane. This can be understood by
exploiting the fact that parameter $\beta_{1}$ is small on the order of
$10^{-3}$. Imposing the regular perturbation ansatz
$\omega=\omega_{0}+\beta_{1}\omega_{1}+\mathcal{O}(\beta_{1}^{2})$ on the
imaginary part of (3) and comparing at zeroth and first order in $\beta_{1}$,
we formally obtain

$\displaystyle\omega_{0}$
$\displaystyle=\sqrt{\alpha_{0}-\frac{\alpha_{1}^{2}}{2}+\sqrt{\left(\alpha_{0}-\frac{\alpha_{1}^{2}}{2}\right)^{2}+\beta_{2}^{2}-\alpha_{0}^{2}}},$
(10) $\displaystyle\omega_{1}$
$\displaystyle=\frac{\omega_{0}}{\alpha_{1}-\alpha_{0}+\omega_{0}^{2}}\tau_{I}\leq\frac{1}{2\sqrt{\alpha_{1}-\alpha_{0}}}\tau_{I}.$
(11)

Thus, we can approximate $\mathbf{H}$ to first order in $\beta_{1}$

$\displaystyle\mathbf{H}(\omega_{0}+\beta_{1}\omega_{1}(\tau_{I}))$
$\displaystyle\approx\left(\tau_{I},\tau_{2}(\omega_{0}+\beta_{1}\omega_{1}(\tau_{I}))-\tau_{I}\right)$
(12)

by using the expression $\tau_{2}(\omega)$ in (4). As a result, $\mathbf{H}$
approaches the graph of the function
$\tau_{I}\mapsto\tau_{2}(\omega_{0})-\tau_{I}$ with slope $-1$ as
$\beta_{1}\to 0$, which can be considered as a zero order approximation of
$\mathbf{H}$.

### On-off periodic infusion

Periodic shot infusion is modeled by a smooth, periodic, quickly varying
function between $G_{\text{in}}(t)=I_{\text{in}}(t)=0$ in mg/(dl$\cdot$min)
(no infusion) and $G_{\text{in}}(t)=G_{\max}$ (glucose infusion),
respectively. We consider the specific form

$\begin{split}G_{\text{in}}(t)=G_{\max}\cdot
s(t-\sigma_{G}),\,s(t)=h(\sin(2\pi t/T_{\text{in}}))\cdot
h(\sin(2\pi(t-t_{\text{in}})/T_{\text{in}}-\pi)),\end{split}$ (13)

where the sigmoidal function $h(y)=(1+\exp(-ky))^{-1}$ can be considered a
smooth version of the Heaviside step function $H(y)=0$ if $y<0$, and $H(y)=1$
if $y\geq 0$ for sufficiently large $k$. The form of (13) was inspired by a
model for auditory perception[44]. $T_{\text{in}}$ is the time between
consecutive shots of glucose/insulin with duration $t_{\text{in}}$, and the
lag $\sigma_{G}$ can be used to specify the timing of the infusion with
respect to the underlying oscillation. The parameter $k$ models the initial
and terminal variations at the beginning and end of the shot application. A
detailed analysis of the influence of the $k$-parameter is beyond the scope of
this paper. We found the choice $k=100$ to be sufficient. Panels (a) and (b)
of Figure 6 show the shape of the infusion patterns used for numerical
computation of time series in Fig. 1.

Figure 6: Form of glucose infusion used to obtain Figs. 1(d)–(e).

### Numerical bifurcation analysis of time-delay systems with periodic
infusion

Numerical simulations have been obtained using pydelay [45]. Numerical
bifurcation analysis has been performed using the software package DDE-BIFTOOL
for Matlab/Octave[33]. For a general introduction to numerical continuation
methods available for delay differential equations and their application to
physiological systems see Refs. 46 and [32], respectively. Isocurves in Fig. 2
have been computed using numerical continuation of periodic orbits in two
parameters with the additional condition fixed period (a), and fixed maximum
value (b), fixed maximum value (c), where in cases (b) and (c) we also relaxed
the phase condition. For bifurcation analysis in the presence of periodic
infusion, we append the two-dimensional ordinary differential equation

$\displaystyle x^{\prime}(t)$ $\displaystyle=x-\omega
y(t)-x(t)(x(t)^{2}+y(t)^{2}),$ (14) $\displaystyle y^{\prime}(t)$
$\displaystyle=-\omega x(t)+y-y(t)(x(t)^{2}+y(t)^{2}),$ (15)

with known stable periodic solution $(x(t),y(t))=(\cos(\omega t),\sin(\omega
t))$ to system (1)–(2). The method has been employed in several other works,
see for example Ref. 47. We achieve the specific form of infusion (13) by
setting
$G_{\text{in}}(t)=G_{\max}h(y(t-\sigma_{G}))h(y(t-\sigma_{G}-t_{in}\omega-\pi))$,
where $\omega=2\pi/T_{\text{in}}$.

## References

  * [1] Goldbeter, A. & Yan, J.  Multi-synchronization and other patterns of multi-rhythmicity in oscillatory biological systems.  _Interface Focus_ 12, 20210089 (2022).
  * [2] Keener, J. & Sneyd, J.  _Mathematical physiology: II: Systems physiology_ (Springer, 2009).
  * [3] Sweeney, E. L., Jeromson, S., Hamilton, D. L., Brooks, N. E. & Walshe, I. H.  Skeletal muscle insulin signaling and whole-body glucose metabolism following acute sleep restriction in healthy males.  _Physiological reports_ 5, e13498 (2017).
  * [4] Sweeney, E. L., Peart, D. J., Ellis, J. G. & Walshe, I. H.  Impairments in glycaemic control do not increase linearly with repeated nights of sleep restriction in healthy adults: a randomised controlled trial.  _Applied Physiology, Nutrition, and Metabolism_ 46, 1091–1096 (2021).
  * [5] Spiga, F., Walker, J. J., Terry, J. R. & Lightman, S. L.  Hpa axis-rhythms.  _Comprehensive physiology_ 4, 1273–1298 (2011).
  * [6] Grant, A. D., Wilsterman, K., Smarr, B. L. & Kriegsfeld, L. J.  Evidence for a coupled oscillator model of endocrine ultradian rhythms.  _Journal of biological rhythms_ 33, 475–496 (2018).
  * [7] Scheen, A. J., Byrne, M. M., Plat, L., Leproult, R. & Van Cauter, E.  Relationships between sleep quality and glucose regulation in normal humans.  _American Journal of Physiology-Endocrinology And Metabolism_ 271, E261–E270 (1996).
  * [8] O’Meara, N. M., Sturis, J., Van Cauter, E., Polonsky, K. S. _et al._ Lack of control by glucose of ultradian insulin secretory oscillations in impaired glucose tolerance and in non-insulin-dependent diabetes mellitus.  _The Journal of clinical investigation_ 92, 262–271 (1993).
  * [9] Committee, A. D. A. P. P.  6\. Glycemic Targets: Standards of Medical Care in Diabetes—2022.  _Diabetes Care_ 45, S83–S96, DOI: 10.2337/dc22-S006 (2021).  https://diabetesjournals.org/care/article-pdf/45/Supplement_1/S83/637560/dc22s006.pdf.
  * [10] Ajmera, I., Swat, M., Laibe, C., Le Novere, N. & Chelliah, V.  The impact of mathematical modeling on the understanding of diabetes and related complications.  _CPT: pharmacometrics & systems pharmacology_ 2, 1–14 (2013).
  * [11] Makroglou, A., Li, J. & Kuang, Y.  Mathematical models and software tools for the glucose-insulin regulatory system and diabetes: an overview.  _Applied numerical mathematics_ 56, 559–573 (2006).
  * [12] Palumbo, P., Ditlevsen, S., Bertuzzi, A. & De Gaetano, A.  Mathematical modeling of the glucose–insulin system: A review.  _Mathematical biosciences_ 244, 69–81 (2013).
  * [13] Huard, B. & Kirkham, G.  Mathematical modelling of glucose dynamics.  _Current Opinion in Endocrine and Metabolic Research_ 100379 (2022).
  * [14] Bergman, R. N., Ider, Y. Z., Bowden, C. R. & Cobelli, C.  Quantitative estimation of insulin sensitivity.  _American Journal of Physiology-Endocrinology And Metabolism_ 236, E667 (1979).
  * [15] Bergman, R. N.  Origins and history of the minimal model of glucose regulation.  _Frontiers in Endocrinology_ 11, 583016 (2021).
  * [16] Ha, J., Muniyappa, R., Sherman, A. S. & Quon, M. J.  When minmod artifactually interprets strong insulin secretion as weak insulin action.  _Frontiers in physiology_ 508 (2021).
  * [17] Walker, J. J., Terry, J. R. & Lightman, S. L.  Origin of ultradian pulsatility in the hypothalamic–pituitary–adrenal axis.  _Proceedings of the Royal Society B: Biological Sciences_ 277, 1627–1633 (2010).
  * [18] Zavala, E. _et al._ Mathematical modelling of endocrine systems.  _Trends in Endocrinology & Metabolism_ 30, 244–257 (2019).
  * [19] Sturis, J. _et al._ Phase-locking regions in a forced model of slow insulin and glucose oscillations.  _Chaos: An Interdisciplinary Journal of Nonlinear Science_ 5, 193–199 (1995).
  * [20] Li, J., Kuang, Y. & Mason, C. C.  Modeling the glucose–insulin regulatory system and ultradian insulin secretory oscillations with two explicit time delays.  _Journal of theoretical biology_ 242, 722–735 (2006).
  * [21] Chen, C.-L., Tsai, H.-W. & Wong, S.-S.  Modeling the physiological glucose–insulin dynamic system on diabetics.  _Journal of theoretical biology_ 265, 314–322 (2010).
  * [22] Cohen, R. B. & Li, J.  A novel model and its analysis on the metabolic regulations of glucose, insulin, and glucagon.  _SIAM Journal on Applied Mathematics_ 81, 2684–2703 (2021).
  * [23] Glass, D. S., Jin, X. & Riedel-Kruse, I. H.  Nonlinear delay differential equations and their application to modeling biological network motifs.  _Nature communications_ 12, 1788 (2021).
  * [24] McKenna, J. P., Dhumpa, R., Mukhitov, N., Roper, M. G. & Bertram, R.  Glucose oscillations can activate an endogenous oscillator in pancreatic islets.  _PLoS computational biology_ 12, e1005143 (2016).
  * [25] Bruce, N. _et al._ Coordination of pancreatic islet rhythmic activity by delayed negative feedback.  _American Journal of Physiology-Endocrinology and Metabolism_ 323, E492–E502, DOI: 10.1152/ajpendo.00123.2022 (2022).
  * [26] Li, J. & Kuang, Y.  Analysis of a model of the glucose-insulin regulatory system with two delays.  _SIAM Journal on Applied Mathematics_ 67, 757–776 (2007).
  * [27] Huard, B., Easton, J. F. & Angelova, M.  Investigation of stability in a two-delay model of the ultradian oscillations in glucose–insulin regulation.  _Communications in Nonlinear Science and Numerical Simulation_ 26, 211–222 (2015).
  * [28] Huard, B., Bridgewater, A. & Angelova, M.  Mathematical investigation of diabetically impaired ultradian oscillations in the glucose–insulin regulation.  _Journal of theoretical biology_ 418, 66–76 (2017).
  * [29] Sturis, J., Polonsky, K., Mosekilde, E. & Van Cauter, E.  Computer model for mechanisms underlying ultradian oscillations of insulin and glucose.  _Am J Physiol_ 260, E801–E809 (1991).
  * [30] Shi, X., Kuang, Y., Makroglou, A., Mokshagundam, S. & Li, J.  Oscillatory dynamics of an intravenous glucose tolerance test model with delay interval.  _Chaos: An Interdisciplinary Journal of Nonlinear Science_ 27, 114324 (2017).
  * [31] Kuznetsov, Y. A.  _Elements of applied bifurcation theory_ , vol. 112 (Springer, 1998).
  * [32] Engelborghs, K., Luzyanina, T. & Roose, D.  Numerical bifurcation analysis of delay differential equations using dde-biftool.  _ACM Transactions on Mathematical Software (TOMS)_ 28, 1–21 (2002).
  * [33] Sieber, J., Engelborghs, K., Luzyanina, T., Samaey, G. & Roose, D.  Dde-biftool v. 3.0 manual—bifurcation analysis of delay differential equations.  _arXiv preprint arXiv:1406.7144_ (2014).
  * [34] Satin, L. S., Butler, P. C., Ha, J. & Sherman, A. S.  Pulsatile insulin secretion, impaired glucose tolerance and type 2 diabetes.  _Molecular aspects of medicine_ 42, 61–77 (2015).
  * [35] Panunzi, S., Palumbo, P. & De Gaetano, A.  A discrete single delay model for the intra-venous glucose tolerance test.  _Theoretical Biology and Medical Modelling_ 4, 35, DOI: 10.1186/1742-4682-4-35 (2007).
  * [36] Li, J., Wang, M., Gaetano, A. D., Palumbo, P. & Panunzi, S.  The range of time delay and the global stability of the equilibrium for an IVGTT model.  _Mathematical Biosciences_ 235, 128–137, DOI: 10.1016/j.mbs.2011.11.005 (2012).
  * [37] Zavala, E.  Misaligned hormonal rhythmicity: Mechanisms of origin and their clinical significance.  _Journal of Neuroendocrinology_ 34, e13144 (2022).
  * [38] Pedersen, M. G., Mosekilde, E., Polonsky, K. S. & Luciani, D. S.  Complex patterns of metabolic and ca2+ entrainment in pancreatic islets by oscillatory glucose.  _Biophysical Journal_ 105, 29–39, DOI: https://doi.org/10.1016/j.bpj.2013.05.036 (2013).
  * [39] Montefusco, F., Cortese, G. & Pedersen, M. G.  Heterogeneous $\alpha$-cell population modeling of glucose-induced inhibition of electrical activity.  _Journal of theoretical biology_ 485, 110036 (2020).
  * [40] Morettini, M. _et al._ Mathematical model of glucagon kinetics for the assessment of insulin-mediated glucagon inhibition during an oral glucose tolerance test.  _Frontiers in endocrinology_ 12, 611147 (2021).
  * [41] Bennett, D. & Gourley, S.  Global stability in a model of the glucose-insulin interaction with time delay.  _European Journal of Applied Mathematics_ 15, 203–221 (2004).
  * [42] Hale, J. K. & Lunel, S. M. V.  _Introduction to functional differential equations_ , vol. 99 (Springer Science & Business Media, 2013).
  * [43] Diekmann, O., Van Gils, S. A., Lunel, S. M. & Walther, H.-O.  _Delay equations: functional-, complex-, and nonlinear analysis_ , vol. 110 (Springer Science & Business Media, 2012).
  * [44] Ferrario, A. & Rankin, J.  Auditory streaming emerges from fast excitation and slow delayed inhibition.  _The Journal of Mathematical Neuroscience_ 11, 1–32 (2021).
  * [45] Flunkert, V. & Schoell, E.  Pydelay-a python tool for solving delay differential equations.  _arXiv preprint arXiv:0911.1633_ (2009).
  * [46] Krauskopf, B. & Sieber, J.  Bifurcation analysis of systems with delays: Methods and their use in applications.  In _Controlling delayed dynamics: advances in theory, methods and applications_ , 195–245 (Springer, 2022).
  * [47] Keane, A., Krauskopf, B. & Postlethwaite, C.  Delayed feedback versus seasonal forcing: Resonance phenomena in an el niño southern oscillation model.  _SIAM Journal on Applied Dynamical Systems_ 14, 1229–1257 (2015).

## Acknowledgements (not compulsory)

The authors thank Jan Sieber for helpful discussions on the implementation of
the numerical continuation methods in DDE-BIFTOOL.

## Author contributions statement

S.R. and B.H contributed equally to the study. Both authors conducted
numerical experiments, analyzed the results, and contributed to the draft of
the manuscript. Both authors reviewed the manuscript.

## Data available statement

Numerical procedures to generate figures are available from the corresponding
author on reasonable request.

## Additional information

The authors declare no competing interests.