Strong delayed negative feedback

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.


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 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 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 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: We first consider the case 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 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.

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

Hopf bifurcation analysis
From the linearized theory, we determine the first Hopf bifurcation located at We may construct a small amplitude periodic solution near τ = τ 0 by using the Lindstedt-Poincaré method (Erneux, 2009;Smith, 2011).We find 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 In other words, the domain where the amplitude of the periodic solution increases parabolically as

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.

Inner solution
We now examine Eqs 3, 4 near t = τ and x = τ.To this end, we introduce the variables s and X defined by We note from Figure 1 (square in the figure) that when t is close to τ. Eqs 3, 4 then implies that the leading order equation for X = X 0 is The solution of this equation needs to satisfy matching conditions as s → ±∞.They are obtained by first introducing ( 14) into (11).We find which implies the condition Second, by introducing ( 14) into (12), we obtain x which leads to the condition The solution of Eq. 16 is 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.
Both conditions requires that The solution Eq.21 now is given by 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 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 ( ) .( 28)

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) 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.(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 The characteristic equation for the growth rate λ is Inserting λ = iω into Eq.31, we obtain from the real part a simple expression for x p given by where z ≡ ωτ.From the imaginary part, we determine the following equation for τ as a function of z and b 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 The approximation of the upper bifurcation point is provided in parametric form (π/2 < z 0 < π is the parameter).

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 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τ ( ), (39) while the period is 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.
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.

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.
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.Using (Eq.48), we first determine the leading approximation of x p .We obtain 3 ) . (50) Second, we evaluate bx using (Eq.48) and (Eq.49).We find 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 The steady state solution u = u(b 1 ) in implicit form is and the conditions for a Hopf bifurcation are.
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.

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 τ.
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

FIGURE 2
FIGURE 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).

FIGURE 3
FIGURE 3 Hopf bifurcation line in the (τ, b) parameter space (p = 20).b a H1 and b a H2 denote the large p analytical approximations determined in the appendix.

Figure 8
Figure8compares 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.

FIGURE 7
FIGURE 7 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 n H2 0.048 is the Hopf bifurcation point obtained numerically from Eqs 3, 5 and b a H2 0.044 π/(2pτ) is its analytical approximation.The fixed parameters are p = 20 and τ = 1.8.

FIGURE 8
FIGURE 8 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 a H1 1.04 matches the bifurcation point determined numerically.Fixed parameters are p = 20 and τ = 1.8.

FIGURE 9
FIGURE 9Periodic 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.