Front. Netw. Physiol. Frontiers in Network Physiology Front. Netw. Physiol. 2674-0109 Frontiers Media S.A. 1399272 10.3389/fnetp.2024.1399272 Network Physiology Original Research Strong delayed negative feedback Erneux 10.3389/fnetp.2024.1399272 Erneux Thomas * Université Libre de Bruxelles, Optique Nonlinéaire Théorique, Bruxelles, Belgium

Edited by: Eckehard Schöll, Technical University of Berlin, Germany

Reviewed by: Cristina Masoller, Universitat Politecnica de Catalunya, Spain

Kathy Lüdge, Technische Universität Ilmenau, Germany

 *Correspondence: Thomas Erneux, thomas.erneux@ulb.be 05 06 2024 2024 4 1399272 11 03 2024 17 04 2024 Copyright © 2024 Erneux. 2024 Erneux

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

In this paper, we analyze the strong feedback limit of two negative feedback schemes which have proven to be efficient for many biological processes (protein synthesis, immune responses, breathing disorders). In this limit, the nonlinear delayed feedback function can be reduced to a function with a threshold nonlinearity. This will considerably help analytical and numerical studies of networks exhibiting different topologies. Mathematically, we compare the bifurcation diagrams for both the delayed and non-delayed feedback functions and show that Hopf classical theory needs to be revisited in the strong feedback limit.

 network physiology delayed negative feedback Mackey-Glass equation delay differential equation hopf bifurcation time periodic oscillations singular perturbation theory section-at-acceptance Networks of Dynamical Systems 1 Introduction

The new multi-disciplinary field of Network Physiology concentrates on coordinated network interactions among distinct organs in the human body (Ivanov et al., 2016; Ivanov, 2021; Schöll et al., 2022). These coordinated network interactions are essential to generating distinct physiological states such as wake, sleep and sleep stages, rest and exercise, stress and anxiety, cognition, consciousness and unconsciousness. Disrupting organ communications can lead to dysfunction of individual systems or trigger a cascade of failures leading to a breakdown and collapse of the entire organism, such as sepsis, coma and multiple organ failure. In Refs. (Bashan et al., 2012; Ivanov et al., 2014), the authors considered a dynamical network consisting of ten nodes representing six physiological systems: brain activity (five EEG waves), cardiac, chin muscle tone, leg and eye movements. They observed changes in network topology during different sleep stages (deep, light, and wake). In addition, they recorded time delays between fluctuations in the output signals of one physiological system, such as cardiovascular, and the emergence of corresponding modulations in another, such as the respiratory. According to the authors, the longer the period during which this delay is constant the stronger the coupling between the two systems.

To develop adequate tools for network physiology, recent efforts focused on understanding the network dynamics of coupled excitable or oscillatory units. Traxl et al. (Traxl et al., 2014) study the effects of noise and global coupling strength on coupled oscillators with different network topologies and different node dynamics. They report a general scaling law for the synchronization of such networks. The inclusion of time delays between interacting nodes has a clear impact on the stability of the network. Inspired by leaky integrate-and-fire models for neuronal networks (Politi and Luccioli, 2010), Mafahim et al. (Mafahim et al., 2015) investigate the dynamics of interacting neurons described by d x i d t = S − γ x i + k ∑ j ≠ i L i j f t − τ j where k is the control parameter and L ij describes the coupling between neurons. Note that i = 1, .N where N is the total number of neurons (nodes). The function f(t) is a Dirac delta-function. Each neuron moves along the x − axis starting at the rest state x = 0 and fires when it reaches the threshold x = 1. When the neuron fires it forces all the neurons linked to it to make a step ahead or backward by the quantity k according to whether L ij = 1 (excitatory) or L ij = −1 (inhibitory). The authors highlight the role of inhibitory links in controlling global network dynamics. While considering a simple delayed coupling mechanism between neurons is reasonable for populations of active neurons, delayed nonlinear feedbacks could be more appropriate as communication mechanisms between distinct organs in the body. The mathematical problem then takes the form d x i d t = g x i + ∑ j ≠ i A i j f x j t − τ j where g(x i ) describes the dynamics of x i in the absence of coupling. The A ij measures the (small or strong) coupling strengths between nodes. The nonlinear function f (x j (t − τ j )) models the delayed feedback of node j with respect to node i. The complexity of the dynamical problem when N > 2 have motivated simplifications which have been explored both analytically and numerically. Networks of delayed coupled Kuramoto oscillators are popular dynamical problems because the state of an oscillator is described by a single angular variable (Laing, 2016; Bick et al., 2020). Another simplification is to consider ring geometries of (unidirectional or bidirectional) coupled nodes (Yuan and Campbell, 2004; Bungay and Campbell, 2007; Ibrahim et al., 2021; Bukh et al., 2023). But the main difficulty remains the fact that we are dealing with coupled delay differential equations (DDEs). If the feedback is strong, however, the feedback function may approach a function exhibiting a threshold nonlinearity which will considerably simplify Eq. 2. In this paper, we consider two delayed negative feedback functions of biological interest and analyze the strong feedback limit. This analysis has never been done and, as we shall demonstrate, Hopf bifurcation theory needs to be revisited.

Negative feedback is one of fundamental mechanisms in cellular networks (Tyson et al., 2003; Tsai et al., 2008; Alon, 2019), which fulfils a variety of functions such as mediating adaptation (Yi et al., 2000; Ma et al., 2009; Ni et al., 2009), stabilizing the abundance of biochemical components (Hasty et al., 2002; Tyson et al., 2003; Alon, 2019), inducing oscillations (Tsai et al., 2008; Elowitz and Leibler, 2000; Kholodenko, 2000; Novak et al., 2007) and decoupling signal and response time (Tyson et al., 2003). Negative feedbacks are shown to be present in many biochemical systems including bacterial adaptation (Yi et al., 2000; Kollmann et al., 2005), mammalian cell cycle (Novak et al., 2010; Ferrell et al., 2011), stress response in yeast (Klipp et al., 2005; Schaber et al., 2012).

A negative feedback control slows of stops a reaction. It may involve a time delay which is needed for signal transduction and transcription, translation and formation of biochemical species (Hoffmann et al., 2002; Börsch and Schaber, 2016). If the delay is too large, however, the control loop loose its landmarks (it does not remember its state so long ago) and exhibit oscillations. The simplest model problem is described by the first order DDE d x d t = f x t − τ − b x where prime means differentiation with respect to time t, x(t) is the state variable, and x (t − τ) is its value at time t − τ. τ > 0 is the delay and b > 0 is a constant that measures the rate to equilibrium in the absence of feedback. The nonlinear function f(x) corresponds to a negative feedback loop: f = 1 if x (t − τ) is small (production is activated) and f = 0 if x (t − τ) is large (production stops).

We first consider the case f x = − tanh κ x  and  b = 0 , and analyze the limit κ → ∞. implying the limit f(x) = ∓1 as x → ±∞. Eqs 3, 4 appear in the modeling of delayed coupled cells (Yuan and Campbell, 2004; Bungay and Campbell, 2007) and for a minimal description of ENSO oscillations (Ghil et al., 2008; Keane et al., 2017). Compared to a purely cubic nonlinearity, the negative feedback function 4) saturates as x increases and is a more realistic feedback function.

We next consider the bifurcation diagram of Eq. 3 with the Hill function f x = 1 1 + x p  and  b > 0 , and analyze the limit p → ∞ implying the limit f(x) → 1 − H (x − 1) where H(y) is the Heaviside function. Originally, Eqs 3, 5 were modeling the control of hematopoiesis (production of blood cells). Proliferation and maturation of blood cells takes time, so there is a delay, τ, between the detection of a deficiency in a circulating population, x, and the appearance in the bloodstream of cells to replenish this population (Mackey and Glass, 1977; Glass and Mackey, 1979). Today, it is known as the Mackey-Glass equation and is considered as a reference DDE for any biological process involving a delayed negative feedback [(Fall et al., 2002) p249, (Beuter et al., 2003) p263, (Milton and Ohira, 2014) p236].

The plan of the paper is as follows. In Section 2, we consider the delayed sigmoidal feedback function 4) and determine the bifurcation diagram of the time-periodic solutions. The diagram shows two distinct domains, namely, one close to the Hopf bifurcation point where the amplitude grows parabolically and a larger domain where the amplitude increases linearly. In Section 3, we analyze the delayed Hill feedback function 5). The bifurcation diagram again exhibits two domains with different oscillatory waveforms. Close to the bifurcation point, the small amplitude oscillations quickly change from harmonic to pulsating oscillations. It motivates the analysis of two singular Hopf bifurcations detailed in Section 3.2. In the last section, we emphasize the role of a delayed exponential function appearing in several negative feedback problems and discuss the limit of large delays as another singular limit of physical interest.

 2 Sigmoidal feedback function

In this section, we analyze Eqs 3, 4 using τ as our bifurcation parameter.

2.1 Hopf bifurcation analysis

From the linearized theory, we determine the first Hopf bifurcation located at τ = τ 0 ≡ π 2 κ . We may construct a small amplitude periodic solution near τ = τ 0 by using the Lindstedt-Poincaré method (Erneux, 2009; Smith, 2011). We find x = 2 τ − τ 0 τ 0 1 / 2 ⁡ cos s + 1 6 τ − τ 0 τ 0 3 / 2 ⁡ cos 3 s + O τ − τ 0 τ 0 5 / 2 . By comparing the first two terms in (Eq. 7) in the limit κ large, we note that this expansion becomes non uniform if (τ − τ 0)/τ 0 = O (1), or equivalently, if τ − τ 0 = O κ − 1 . In other words, the domain where the amplitude of the periodic solution increases parabolically as x = ± 2 τ − τ 0 τ 0 1 / 2 is only valid if τ − τ 0 ≪ κ −1.

 2.2 Sawtooth oscillations

By contrast to our Hopf bifurcation analysis where we were looking for a small amplitude solution and then investigated its behavior for large κ, we now seek a periodic solution of arbitrary amplitude but take advantage of the large value of κ. A typical numerical solution for κ = 10 and τ = 0.4 > τ 0 = 0.157 is shown in Figure 1. This solution consists of a succession of straight lines connected at extrema located at t = (1 + 2n)τ (n = 0, 1, … ). It motivates to construct an analytical solution by using the method of matched asymptotic expansions (Kevorkian and Cole, 1996; Bender and Orszag, 1999; O’Malley, 2014). The method considers two distinct approximations valid for different intervals of time. The outer approximation, valid for a large subdomain, is obtained by treating the problem as a regular perturbation problem. The inner approximation solves a separate perturbation problem valid in a small subdomain where the outer solution is inaccurate. This area is often referred to as a transition layer. Outer and inner solutions are then combined through a process called “matching” in such a way that a solution for the whole domain is obtained.

Periodic solution of Eqs 3, 4 for κ = 10 and τ = 0.4. The figure shows x(t) (black), x (t − τ) (red), and the leading asymptotic approximation provided by (Eqs 26, 27) (grey). The square shows x (t − τ) ≃ t − τ when t is close to τ.

 2.3 Outer solution

Noting that tanh (κx) = 1 if κx ≫ 1 and tanh (κx) = −1 if κx ≪ 1, the leading approximation of Eqs 3, 4 satisfies d x 0 d t = − 1  if  x 0 t − τ > 0 1  if  x 0 t − τ < 0 . Consequently, x 0(t) is alternatively increasing and decreasing as x 0 = t 0 < t < τ , x 0 = τ − t − τ τ < t < 3 τ , x 0 = − τ + t − 3 τ 3 τ < t < 4 τ , and so on. Eq. 10 has been studied by Fridman et al. (Fridman et al., 2002) who showed that only the 4τ-periodic solution is stable, whereas the 4τ/(4n + 1) -periodic oscillations (n = 1, 2, … ) are unstable.

 2.4 Inner solution

We now examine Eqs 3, 4 near t = τ and x = τ. To this end, we introduce the variables s and X defined by t = τ + κ − 1 s ,   x = τ + κ − 1 X We note from Figure 1 (square in the figure) that x t − τ = t − τ = κ − 1 s when t is close to τ. Eqs 3, 4 then implies that the leading order equation for X = X 0 is d X 0 d s = − tanh s . The solution of this equation needs to satisfy matching conditions as s → ±∞. They are obtained by first introducing (14) into (11). We find x = τ + κ − 1 X 0 = t = τ + κ − 1 s which implies the condition X 0 = s   s → − ∞ . Second, by introducing (14) into (12), we obtain x 0 = τ + κ − 1 X 0 = τ − t − τ = τ − κ − 1 s which leads to the condition X 0 = − s   s → ∞ . The solution of Eq. 16 is X 0 = − ln cosh s + C where C is a constant of integration. We examine the limits s → ±∞ of (Eq. 21) which need to match (18) and (20). We find the conditions. X 0 s → − ∞ → − ln exp − s / 2 + C = s + ln 2 + C = s X 0 s → ∞ → − ln exp s / 2 + C = − s + ln 2 + C = − s . Both conditions requires that C = − ln 2 . The solution Eq. 21 now is given by X 0 = − ln 2 ⁡ cosh s .

Figure 2 represents the numerical bifurcation diagram of the periodic solutions for κ = 10 together with Hopf local approximation 9) and the large κ approximation given by (Eq. 28). Similar inner solutions may be constructed for the other extrema. An uniform solution combining outer and inner solutions leads to. x = τ − κ − 1 ⁡ ln ( 2 ⁡ cosh κ t − τ   0 < t < 2 τ , x = − τ + κ − 1 ⁡ ln ( 2 ⁡ cosh κ t − 3 τ   2 τ < t < 4 τ and so on. These approximations are compared to the numerical solution (grey line in Figure 1). The reason for such good agreement comes from the fact that the first correction to the leading outer approximation x = x 0(t) is not O (κ −1) but much smaller like O (exp (−κ)). This is because the expansion of tanh (κx) as κx → ±∞ is tan(κx) = ±1–2 exp (∓κx) + as κx → ±∞. The extrema of the oscillations are obtained from (Eq. 26) and (Eq. 27) at t = τ and t = 3τ, respectively: x = ± τ − κ − 1 ⁡ ln 2 .

Bifurcation diagram of the periodic solutions of Eqs. 3, 4 for κ = 10. The numerical bifurcation diagram of the extrema (black) is compared to Hopf local approximation (9) (blue). The straight lines (red) correspond to the approximation (28).

 3 Hill feedback function

By the end of the seventies two independent papers devoted to the development of red blood cells generated considerable mathematical interest. The paper by Wazewska-Czyzewska and Lasota (Wazewska-Czyzewska and Lasota, 1976) and the one by Mackey and Glass (Mackey and Glass, 1977) appeared in 1976 and 1977, respectively. Without knowing each other at that time, these authors published almost simultaneously two models very similar in several points. The one from Wazewska and Lasota is given by (Wazewska-Czyzewska and Lasota, 1976) d x d t = a ⁡ exp − c x t − τ − b x where a, b, and c are all positives. The other, today known as one of the two Mackey-Glass equations, is given by Eqs. 3, 5 (Mackey and Glass, 1977; Glass and Mackey, 1979) where p > 0 and b > 0. The Wazewska-Lasota Eq. 29 was derived from an age structured partial differential equation, and delay was a consequence of its integration. On the other hand, the Mackey-Glass equation Eqs. 3, 5 had been set up directly into a delay differential equation. The nonlinear function 5) is Hill function which is based on the law of mass action for the binding of molecules (Milton and Ohira, 2014). Eqs. 3, 5 has been the source of many numerical and analytical studies. In particular, the limit of a strong feedback (p → ∞) allows to simplify 5) and obtain an analytical approximation. Our objective is to compare its bifurcation diagram with the one obtained numerically from the original DDE with a fixed value of p. As we shall demonstrate, the agreement between the two diagrams is excellent except near the Hopf bifurcation points.

Eqs. 3, 5 admit a unique steady state which is unstable if b H1 < b < b H2. The critical points b = b H1 and b = b H2 are Hopf bifurcation points. Their analytical determination is documented at several places (Fall et al., 2002) p249, (Milton and Ohira, 2014), p243 and we briefly detail their conditions. From the steady state equation, we first determine b as a function of x b = 1 x 1 + x p . The characteristic equation for the growth rate λ is λ = − p x p − 1 1 + x p 2 exp − λ τ − b . Inserting λ = iω into Eq. 31, we obtain from the real part a simple expression for x p given by x p = − 1 p ⁡ cos z + 1 > 0 where z ≡ ωτ. From the imaginary part, we determine the following equation for τ as a function of z and b τ = − z b ⁡ tan z . Eqs 30, 32, 33 are the equations defining the Hopf bifurcation in parameter space. Using Eq. 32 for x p and determining x = ( x p ) 1 / p for x, we obtain b = b(z) from Eq. 30. The expression for b is then introduced in Eq. 33 allowing us to determine τ = τ(z). By continuously increasing z (π/2 < z < π), we determine the Hopf bifurcation line in the (τ, b) parameter space. See Figure 3. The lines denotes by b H 1 a and b H 2 a are the large p approximations of the upper and lower parts of the Hopf bifurcation line. They are determined in the appendix and their expressions will be useful in the next sections. The lowest Hopf bifurcation point admits the approximation b H 2 a = π 2 p τ . The approximation of the upper bifurcation point is provided in parametric form (π/2 < z 0 < π is the parameter). τ = − z 0 tan z 0 , b H 1 a = 1 + ln p p − 1 p x 1 + exp x 1 , x 1 = − ln − cos z 0       π / 2 < z 0 < π .

Hopf bifurcation line in the (τ, b) parameter space (p = 20). b H 1 a and b H 2 a denote the large p analytical approximations determined in the appendix.

 3.1 Bifurcation diagrams

By the end of the seventies and early eighties, Mathematicians discovered that piecewise linear (Glass and Mackey, 1979; Mackey and an der Heiden, 1984) or piecewise constant functions (An der Heiden and Walther, 1983) as nonlinearities can make dynamics generated by a scalar delay differential equation accessible, and that one can compute periodic solutions explicitly. In the large p limit, the nonlinear function 5) approaches the function 1 − H (x − 1) where H(y) is the Heaviside step function. Consequently, Eqs 3, 5 simplify as d x d t = − b x + 0 if  x t − τ > 1 1  if  x t − τ < 1 . A typical periodic solution of Eq. 38 is shown in Figure 4. Eq. 38 consists of a pair of ordinary differential equations which can be solved by the method of steps (An der Heiden and Mackey, 1982). The application of the method is well documented in (Mackey and Glass, 1977; Mackey et al., 1996). The method is also used for a delayed negative feedback problem (Milton, 2003) modeling changes in pupil size. The periodic solution consists of increasing and decreasing exponentials. The extrema of the oscillations are given by. x min = exp − b τ , x max = 1 − b − 1 exp − b τ + b − 1 while the period is P = − b − 1 ⁡ ln 1 − b x max 1 − b x min . x min x max .

Periodic solution of Eq. 38. b = 0.4 and τ = 1.8.

The expressions (Eq. 39) and (Eq. 40) for the extrema and the Period (Eq. 41) are compared to the numerical bifurcation diagrams obtained from Eqs 3, 5 with p = 20. See Figure 5.

Bifurcation diagram. The fixed parameters are τ = 1.8 and p = 20. The black lines show the extrema and the period of the limit-cycle oscillations of Eqs 3, 5. The red lines are the approximations of the extrema and period provided by Eqs 39–41.

The same construction of the solution is proposed in Ref. (Coombes and Laing, 2009). but with τ as the bifurcation parameter instead of b. The amplitude of the oscillations x max − x min increases like τ and saturates at a fixed value as τ → ∞.

The analytical approximations obtained in the limit p → ∞ correctly match the bifurcation branches obtained numerically from Eqs 3, 5 except near the two Hopf bifurcation points where the period becomes infinite. According to Hopf bifurcation theory, the oscillations near the bifurcation point should be nearly sinusoidal and exhibit a fixed period. So how may we understand the radical change of the oscillations from harmonic to pulsating in the vicinity of the two Hopf bifurcation points ? To resolve this problem, we need to take into account the large value of p in the construction of a small amplitude solution near each bifurcation points. To this end, we plan to scale the deviation b − b H with respect to p −1 and then reexamine the large p limit.

 3.2 Singular hopf bifurcations

We note from Eq. 39 and Eq. 40 that if b → 0+, x min → 1, x max → 1 + τ, and P → b −1 ln (1 + τ) → ∞. On the other hand, if b → 1−, x min → exp (−τ), x max → 1+, and P → − ln (1 − b) → ∞. Eq. 38 fails to provide the solution of Eqs 3, 5 near b = 0 and b = 1 because the period P become infinite at these points. We also need to realize that our analytical construction of the limit-cycle assumed that x(t) is sequentially larger and less that 1. This is not the case near the two Hopf bifurcation points where the oscillations remains either above or below 1. Figure 6 shows the limit-cycle oscillations obtained numerically from Eqs 3, 5 for b slightly above b H2 ≃ 0.048. The oscillations are sinusoidal for b = 0.05 and are clearly above x = 1 while the oscillations for b = 0.1 have their minima close to x = 1.

Limit-cycle oscillations close to the Hopf bifurcation point b = b H2 = 0.048. p = 20 and τ = 1.8. The value of b is indicated in the figure. The oscillations for b = 0.05 are sinusoidal with a period P = 3.89 close to the Hopf bifurcation period P H2 ≃ 2π/(π/2 + p −1) = 3.88. The oscillations for b = 0.1 exhibit minima slightly below x = 1 and the waveform approaches two successive exponentials. The period has increased and equals P = 4.88.

Our asymptotic theory based on the large value of p needs to be revised near the two Hopf bifurcation points. We first consider the lower Hopf bifurcation point b = b H2 ∼ 0 for which the analysis is simpler than the case b = b H1 ∼ 1.

3.2.1 <italic>b</italic> = <italic>b</italic> <sub> <italic>H</italic>2</sub> ∼ 0

The analysis of the Hopf bifurcation point detailed in the appendix suggests that x p = O(p) and b = O(p −1). We introduce the new bifurcation parameter b 1 = O(1) defined by b = p − 1 b 1 and take into account that x p (t − τ) is an O(p) large quantity. Eqs 3, 5 then simplifies as x ′ = x − p t − τ 1 + O p − 1 − p − 1 b 1 x . We next introduce the new dependent variable u defined by x = 1 + p − 1 ⁡ ln p + p − 1 u where the p −1 ln(p) term is motivated by the expansion of x at the Hopf bifurcation b = b H2 (see Appendix). We determine x −p (t − τ) and obtain 1 x − p t − τ = 1 p exp − u t − τ . From Eq. 43, we then find that the leading order problem is O (p −1) and is given by u ′ = exp − u t − τ − b 1 .

Eq. 46 belongs to the family of Wright’s equation (Wright’s equation is Eq. 46 with b 1 = 1). It admits a Hopf bifurcation at b 1 = π/(2τ). The bifurcation diagram of Eq. 46 is shown in terms of the extrema of x in Figure 7 2 . The agreement between the minima of the oscillations is excellent but the maxima diverges as soon as b > 0.06.

Bifurcation diagrams near b = b H2. The black lines correspond to the bifurcation diagram of Eqs 3, 5. The red dots mark the bifurcation diagram of Eq. 46. b H 2 n = 0.048 is the Hopf bifurcation point obtained numerically from Eqs 3, 5 and b H 2 a = 0.044 = π / ( 2 p τ ) is its analytical approximation. The fixed parameters are p = 20 and τ = 1.8.

 3.2.2 <italic>b</italic> = <italic>b</italic> <sub> <italic>H</italic>1</sub> ∼ 1

The analysis of the Hopf bifurcation point detailed in the appendix indicates that x p = O (p −1) for the upper Hopf bifurcation branch. Eqs 3, 5 then simplifies as x ′ = 1 − x p t − τ + O p − 2 − b x We introduce the new dependent variable u and new control parameter b 1 = O(1) as. x = 1 − ln p p + 1 p u , b = 1 + ln p p + 1 p b 1 where the ln(p)/p correction term is motivated by the asymptotic expressions of x and b at b = b H1 (see Supplementary Appendix). Using (Eq. 48), we first determine the leading approximation of x p . We obtain 3 x p t − τ = 1 p exp u t − τ . Second, we evaluate bx using (Eq. 48) and (Eq. 49). We find b x = 1 + 1 p b 1 + u . Inserting (Eq. 48), (Eq. 50), and (Eq. 51) into Eq. 47, we find that the leading problem for u is O(p −1) and is given by u ′ = − exp u t − τ − u + b 1 . The steady state solution u = u(b 1) in implicit form is b 1 = − u + exp u and the conditions for a Hopf bifurcation are. cos z = − exp − u , τ = − z / tan z .

The expression (Eq. 49) with (Eq. 53) and x 1 replacing u is identical to Eq. 64 in the appendix. Eqs 54, 55 are identical to (66) and (63) in the appendix with x 1 replacing u and z 0 replacing z.

Figure 8 compares the bifurcation diagram of the original equations (Eqs 3, 5) and the bifurcation diagram obtained using the reduced Eq. 52. The agreement between the maxima is excellent but the minima quickly diverges as we deviate from the Hopf bifurcation point.

Bifurcation diagrams near b = b H1. The black lines correspond to the bifurcation diagram of Eqs 3, 5. The red dots mark the bifurcation diagram of Eq. 52. The analytical Hopf bifurcation point at b H 1 a = 1.04 matches the bifurcation point determined numerically. Fixed parameters are p = 20 and τ = 1.8.

 4 Discussion

The new field of network physiology is based on the observation that a healthy body requires good synchronization between different organs. When perturbing elements disturb this equilibrium, many physiological processes are changing from metabolism, immune function to cardiovascular regulation. An example of a simple and well studied network is the circadian network. Circadian rhythms are generated by the autonomous circadian clock, the suprachiasmatic nucleus (SCN), and clock genes that are present in all tissues (Buijs et al., 2016). The SCN times these peripheral clocks, as well as behavioral and physiological processes. Recent studies have shown that frequent violations of conditions set by our biological clock, such as shift work, jet lag, sleep deprivation, or simply eating at the wrong time of the day, may have deleterious effects on health. On the long run, these perturbations are desynchronizing the circadian network.

In this paper, we hypothesize that strong delayed negative feedback loops between elements of the network are essential for a good synchronization. This idea is motivated by the importance of negative feedback in cellular processes. We have considered two delayed negative feedback which have proven to be useful for combined analytical and numerical studies. The limit of strong feedback allows to reduce the delayed function to a function exhibiting a threshold nonlinearity. We have shown that Hopf bifurcation theory needs to be revisited in the case of a strong negative feedback. By treating the Hopf problem as a singular perturbation problem, we determine small amplitude solutions which are quickly changing waveforms as we deviate from the bifurcation point. Like Wazewska and Lasota Eq. 29, the reduced problems for the two Hopf bifurcations of Mackey-Glass equation (Eqs 3, 5) exhibit a delayed exponential nonlinearity. The latter also appeared in a minimal model for periodic or episodic star formation (Alice et al., 2008).

The singularity of the Hopf bifurcation caused by the strong feedback limit is not the only one of physical interest. The limit of large delay is another case where harmonic oscillations quickly become 2τ − periodic square-waves as we deviate from the Hopf bifurcation point (Erneux et al., 2004). From the analytical solution of Eq. 38, we note that the square-wave is switching from 0 to b −1 through fast transition layers consisting of decaying exponentials. Figure 9 shows the time-periodic solution of Mackey-Glass equations Eqs 3, 5 for a large value of the delay τ.

Periodic solution obtained numerically from Eqs 3, 5. The parameters are τ = 100, b = 0.1, and p = 20. The red curves are the large p approximations obtained from solving Eq. 38. Fast transition layers appears at t = 0, t = t m and t = P.

 Author contributions

TE: Writing–original draft.

 Funding

The author(s) declare that no financial support was received for the research, authorship, and/or publication of this article.

The author acknowledges useful discussions with MC Mackey.

 Conflict of interest

The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

 Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

 Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnetp.2024.1399272/full#supplementary-material

ln ( x − p ) = − p ⁡ ln ( x ) = − p ⁡ ln ( 1 + p − 1 ⁡ ln ( p ) + p − 1 u ) = − p ( p − 1 ⁡ ln ( p ) + p − 1 u ) and thus: x −p = p −1 exp (−u)

In our simulations, we consider the logistic equation equivalent to Eq. 46 after the change of variables u = − ln(v). It is given by v ′ = v b 1 − v t − τ . and x = x(v) then is x = 1 + ln p / v 1 / p .

ln ( x p ) = p ⁡ ln ( x ) = p ⁡ ln ( 1 − ln ( p ) p + 1 p u ) = p − ln ( p ) p + 1 p u = − ln(p) + u and thus: x p = p −1 exp(u)

ln(x p ) = p ln(x)

= p ⁡ ln 1 + p − 1 ⁡ ln ( p ) + p − 1 x 1 + ⋯

= ln(p) + x 1 + ⋯ as p → ∞. Thus: x p = exp ln ( p ) + x 1 + ⋯ = p ⁡ exp ( x 1 ) + ⋯

ln(x p ) = p ln(x)

= p ⁡ ln 1 − p − 1 ⁡ ln ( p ) + p − 1 x 1 + ⋯

= − ln(p) + x 1 + ⋯ as p → ∞. Thus: x p = exp ln ( p − 1 ) + x 1 + ⋯ = p − 1 ⁡ exp ( x 1 ) + ⋯

References Alice C. Quillen A. C. Bland-Hawthorn J. (2008). When is star formation episodic? A delay differential equation ‘negative feedback’ model. Mon. Not. R. Astron. Soc. 386, 2227–2234. 10.1111/j.1365-2966.2008.13193.x Alon U. (2019) An introduction to systems biology: design principles of biological circuits. Boca Raton: Chapman and Hall/CRC. An der Heiden U. Mackey M. C. (1982). The dynamics of production and destruction: analytic insight into complex behavior. J. Math. Biol. 16, 75–101. 10.1007/bf00275162 An der Heiden U. Walther H. O. (1983). Existence of chaos in control systems with delayed feedback. J. Diff. Equations 47, 273–295. 10.1016/0022-0396(83)90037-2 Bashan A. Bartsch R. P. Kantelhardt J. W. Havlin S. Ivanov P. C. (2012). Network physiology reveals relations between network topology and physiological function. Nat. Commun. 3, 702. 10.1038/ncomms1705 Bender C. M. Orszag S. A. (1999) Advanced mathematical methods for scientists and engineers. New York: Springer. Beuter A. Glass L. Mackey M. C. Titcombe M. (2003) Nonlinear dynamics in physiology and medicine. New York: Springer. Bick C. Goodfellow M. Laing C. R. Martens E. A. (2020). Understanding the dynamics of biological and neural oscillator networks through exact mean-field reductions: a review. J. Math. Neurosci. 10 (1), 9. 10.1186/s13408-020-00086-9 Börsch A. B. Schaber J. (2016). How time delay and network design shape response patterns in biochemical negative feedback systems. BMC Syst. Biol. 10, 82. 10.1186/s12918-016-0325-9 Buijs F. N. León-Mercado L. Guzmán-Ruiz M. Guerrero-Vargas N. N. Romo-Nava F. Buijs R. M. (2016). The circadian system: a regulatory feedback network of periphery and brain. Physiol. (Bethesda) 31 (3), 170–181. 10.1152/physiol.00037.2015 Bukh A. V. Shepelev I. A. Elizarov E. M. Muni S. S. Schöll E. Strelkova G. I. (2023). Role of coupling delay in oscillatory activity in autonomous networks of excitable neurons with dissipation. Chaos 33, 073114. 10.1063/5.0147883 Bungay S. Campbell S. A. (2007). Patterns of oscillation in a ring of identical cells with delayed coupling. Int. J. Bifurcation Chaos 17 (9), 3109–3125. 10.1142/s0218127407018907 Coombes S. Laing C. R. (2009). Instabilities in threshold-diffusion equations with delay. Phys. D. 238, 264–272. 10.1016/j.physd.2008.10.014 Elowitz M. B. Leibler S. (2000). A synthetic oscillatory network of transcriptional regulators. Nature 403 (6767), 335–338. 10.1038/35002125 Erneux T. (2009) Applied delay differential equations. New York: Springer. Erneux T. Larger L. Won Lee M. Goedgebuer J. P. (2004). Ikeda Hopf bifurcation revisited. Phys. D. 194, 49–64. 10.1016/j.physd.2004.01.038 Fall C. P. Marland E. S. Wagner J. M. Tyson J. J. (2002) Computational cell biology. New York: Springer. Ferrell J. E. Tsai T. Y.-C. Yang Q. (2011). Modeling the cell cycle: why do certain circuits oscillate? Cell 144 (6), 874–885. 10.1016/j.cell.2011.03.006 Fridman L. Fridman E. Shustin E. (2002). “Steady modes and sliding modes in relay control systems with delay,” in Sliding mode control in engineering. Editors Barbot J. P. Perruguetti W. (New York: Marcel Dekker), 264–295. Ghil M. Zaliapin I. Thompson S. (2008). A delay differential model of ENSO variability: parametric instability and the distribution of extremes. Nonlin. Process. Geophys. 15, 417–433. 10.5194/npg-15-417-2008 Glass L. Mackey M. C. (1979). Pathological conditions resulting from instabilities in physiological control systems. Ann. N. Y. Acad. Sci. 316, 214–235. 10.1111/j.1749-6632.1979.tb29471.x Hasty J. Dolnik M. Rottschäfer V. Collins J. J. (2002). Synthetic gene network for entraining and amplifying cellular oscillations. Phys. Rev. Lett. 88 (14), 148101. 10.1103/PhysRevLett.88.148101 Hoffmann A. Levchenko A. Scott M. L. Baltimore D. (2002). The IkappaB-NF-kappaB signaling module: temporal control and selective gene activation. Science 298 (5596), 1241–1245. 10.1126/science.1071914 Ibrahim M. M. Kamran M. A. Mannan M. M. N. Jung I. H. Kim S. (2021). Lag synchronization of coupled time-delayed FitzHugh–Nagumo neural networks via feedback control. Sci. Rep. 11 (1), 3884. 10.1038/s41598-021-82886-x Ivanov P. C. (2021). The new field of network physiology: building the human physiolome. Front. Netw. Physiol. 1, 711778. 10.3389/fnetp.2021.711778 Ivanov P. C. Bartsch R. P. (2014). “Network physiology: mapping interactions between networks of physiologic networks,” in Networks of networks: the last frontier of complexity. Editors D’Agostino G. Scala A. (Berlin, Germany: Springer International Publishing Switzerland), 203–222. Ivanov P. C. Liu K. K. L. Bartsch R. P. (2016). Focus on the emerging new fields of network physiology and network medicine. New J. Phys. 18, 100201. 10.1088/1367-2630/18/10/100201 Keane A. Krauskopf B. Postlethwaite C. M. (2017). Climate models with delay differential equations. Chaos 27 (11), 114309. 10.1063/1.5006923 Kevorkian J. Cole J. D. (1996) Multiple scale and singular perturbation methods. Berlin: Springer-Verlag. Kholodenko B. N. (2000). Negative feedback and ultrasensitivity can bring about oscillations in the mitogen-activated protein kinase cascades. Eur. J. Biochem. 267 (6), 1583–1588. 10.1046/j.1432-1327.2000.01197.x Klipp E. Nordlander B. Krüger R. Gennemark P. Hohmann S. (2005). Integrative model of the response of yeast to osmotic shock. Nat. Biotechnol. 23 (8), 975–982. 10.1038/nbt1114 Kollmann M. Løvdok L. Bartholomé K. Timmer J. Sourjik V. (2005). Design principles of a bacterial signalling network. Nature 438 (7067), 504–507. 10.1038/nature04228 Laing C. R. (2016). Travelling waves in arrays of delay-coupled phase oscillators. Chaos 26, 094802. 10.1063/1.4953663 Ma W. Trusina A. El-Samad H. Lim W. A. Tang C. (2009). Defining network topologies that can achieve biochemical adaptation. Cell 138 (4), 760–773. 10.1016/j.cell.2009.06.013 Mackey M. C. (1996). “Mathematical models of hematopoietic cell replication and control,” in The art of mathematical modelling: case studies in ecology, physiology and biofluids. Editors Othmer H. G. Adler F. R. Lewis M. A. Dallon J. C. (Hoboken, New Jersey, USA: Prentice Hall), 149–178. Mackey M. C. an der Heiden U. (1984). The dynamics of recurrent inhibition. J. Math. Biol. 19, 211–225. 10.1007/BF00277747 Mackey M. C. Glass L. (1977). Oscillation and chaos in physiological control systems. Science 197, 287–289. 10.1126/science.267326 Mafahim J. U. Lambert D. Zare M. Grigolini P. (2015). Complexity matching in neural networks. New J. Phys. 17, 015003. 10.1088/1367-2630/17/1/015003 Milton J. (2003). “Pupil light reflex: delays and oscillations,” in Nonlinear dynamics in physiology and medicine. Interdisciplinary applied mathematics. Editors Beuter A. Glass L. Mackey M. C. Titcombe M. S. Berlin: Springer. Milton J. Ohira T. (2014) Mathematics as a laboratory tool. New York: Springer. Ni X. Y. Drengstig T. Ruoff P. (2009). The control of the controller: molecular mechanisms for robust perfect adaptation and temperature compensation. Biophys. J. 97 (5), 1244–1253. 10.1016/j.bpj.2009.06.030 Novak B. Kapuy O. Domingo-Sananes M. R. Tyson J. J. (2010). Regulated protein kinases and phosphatases in cell cycle decisions. Curr. Opin. Cell Biol. 22 (6), 801–808. 10.1016/j.ceb.2010.07.001 Novak B. Tyson J. J. Gyorffy B. Csikasz-Nagy A. (2007). Irreversible cell-cycle transitions are due to systems-level feedback. Nat. Cell Biol. 9 (7), 724–728. 10.1038/ncb0707-724 O’Malley R. E. (2014) Historical developments in singular perturbations. New York: Springer. Politi A. Luccioli S. (2010). “Dynamics of networks of leaky-integrate-and-fire neurons,” in Network science: complexity in nature and technology (Berlin, Germany: Springer), 217–242. Schaber J. Baltanas R. Bush A. Klipp E. Colman-Lerner A. (2012). Modelling reveals novel roles of two parallel signalling pathways and homeostatic feedbacks in yeast. Mol. Syst. Biol. 8, 622. 10.1038/msb.2012.53 Schöll E. Sawicki J. Berner R. Ivanov P. C. (2022). Editorial: adaptive networks in functional modeling of physiological systems. Front. Netw. Physiol. 2, 996784. 10.3389/fnetp.2022.996784 Smith H. (2011) An introduction to delay differential equations with applications to the life Sciences. New York: Springer. Traxl D. Boers N. Kurths J. (2014). General scaling of maximum degree of synchronization in noisy complex networks. New J. Phys. 16, 115009. 10.1088/1367-2630/16/11/115009 Tsai T. Y.-C. Choi Y. S. Ma W. Pomerening J. R. Tang C. Ferrell J. E. (2008). Robust, tunable biological oscillations from interlinked positive and negative feedback loops. Science 321 (5885), 126–129. 10.1126/science.1156951 Tyson J. J. Chen K. C. Novak B. (2003). Sniffers, buzzers, toggles and blinkers: dynamics of regulatory and signaling pathways in the cell. Curr. Opin. Cell Biol. 15 (2), 221–231. 10.1016/s0955-0674(03)00017-6 Wazewska-Czyzewska M. Lasota A. (1976). Matematyczne problemy dynamiki ukladu krwinek czerwonych (Mathematical problems of the dynamics of red blood cell population). Available at: https://wydawnictwa.ptm.org.pl/index.php/matematyka-stosowana/article/view/1173. Yi T. M. Huang Y. Simon M. I. Doyle J. (2000). Robust perfect adaptation in bacterial chemotaxis through integral feedback control. Proc. Natl. Acad. Sci. 97 (9), 4649–4653. 10.1073/pnas.97.9.4649 Yuan Y. Campbell S. A. (2004). Stability and synchronization of a ring of identical cells with delayed coupling. J. Dyn. Diff. Equat. 16, 709–744. 10.1007/s10884-004-6114-y Appendix The large <italic>p</italic> limit of the Hopf bifurcation points of Eqs. <xref ref-type="disp-formula" rid="e3">(3)</xref> and <xref ref-type="disp-formula" rid="e5">(5)</xref>

In this appendix, we determine the large p limit of the upper and lower parts of the Hopf bifurcation line b = b(τ) shown in Figure 3. The conditions for the Hopf bifurcation are provided by the steady state equation (30) and Eqs. (32), (33).

 The lower Hopf bifurcation <italic>b</italic> = <italic>b</italic> <sub> <italic>H</italic>2</sub> ≪ 1

Numerical simulations suggest that x p = O(p) and x ∼ 1 for the lower Hopf bifurcation branch. It motivates to seek a solution for x of the form x = 1 + p − 1 ⁡ ln p + p − 1 x 1 + ⋯ where the p −1 ln(p) correction term is needed when we determine x p and x 1 = O(1). We find 4 x p = p ⁡ exp x 1 as the leading approximation. From Eqs. (30) and (32), we then determine b and z ≡ ωτ as b = 1 p exp − x 1 , z = π 2 + p − 1 + p − 2 ⁡ exp − x 1 . Last, we evaluate τ from (33) and obtain τ = π 2 p b , or equivalently, b = b H 2 a ≡ π 2 p τ .

5.2 The upper Hopf bifurcation <italic>b</italic> = <italic>b</italic> <sub> <italic>H</italic>1</sub> ∼ 1

Numerical simulations now suggest that x p = O(p −1) and x ∼ 1 for the upper Hopf bifurcation branch. We seek a solution for x of the form x = 1 − p − 1 ⁡ ln p + p − 1 x 1 + ⋯ where the −p −1 ln(p) correction term is needed when we determine x p and x 1 = O(1). We obtain 5 x p = p − 1 ⁡ exp x 1 . From Eq. (30), we then determine b as b = b H 1 a ≡ 1 + p − 1 ⁡ ln p − p − 1 x 1 + exp x 1 . Inserting z = z 0 + p − 1 z 1 + ⋯ into Eq. (32), we find that the leading equation is cos(z 0) = − exp(−x 1). It then provides an expression for x 1 = x 1(z 0) given by x 1 = − ln − cos z 0      π / 2 < z 0 < π . Last, we evaluate τ from Eq. (33) and find τ = − z 0 tan z 0 . In summary, the large p limit of the upper Hopf bifurcation branch shown in Figure 3 is provided in parametric form by Eq. (64) with x 1 = x 1(z 0) determined from (66), and by Eq. (67) (z 0 is the parameter).