Nonlinear dynamics and entrainment in a continuously forced pulse-modulated model of testosterone regulation

Dynamical behaviors arising in a previously developed pulse-modulated mathematical model of non-basal testosterone regulation in the human male due to continuous exogenous signals are studied. In the context of endocrine regulation, exogenous signals represent, e.g., the influx of a hormone replacement therapy drug, the influence of the circadian rhythm, and interactions with other endocrine loops. This extends the scope of the autonomous pulse-modulated models of endocrine regulation to a broader class of problems, such as therapy optimization, but also puts it in the context of biological rhythms studied in chronobiology. The model dynamics are hybrid since the hormone metabolism is suitably captured by a continuous description and the control feedback is implemented in a discrete (i.e., event-based) manner by the hypothalamus of the brain. It is demonstrated that the endocrine loop with an exogenous signal entering the continuous part can be equivalently described by proper modifications in the pulse modulation functions of the autonomous model. The cases of a constant and a harmonic exogenous signal are treated in detail and illustrated by the results of bifurcation analysis. According to the model, adding a constant exogenous signal only reduces the mean value of testosterone, which result pertains to the effects of hormone replacement therapies under intact endocrine feedback regulation. Further, for the case of a single-tone harmonic positive exogenous signal, bistability and quasiperiodicity arise in the system. The convergence to either of the stationary solutions in a bistable regime is shown to be controlled by the phase of the exogenous signal thus relating this transition to the phenomenon of jet lag.

mone metabolism is suitably captured by a continuous description and the control feedback is implemented in a discrete (i.e., event-based) manner by the hypothalamus of the brain. It is demonstrated that the endocrine loop with an exogenous signal entering the continuous part can be equivalently described by proper modifications in the pulse modulation functions of the autonomous model. The cases of a constant and a harmonic exogenous signal are treated in detail and illustrated by the results of bifurcation analysis. According to the model, adding a constant exogenous signal only reduces the mean value of testosterone, which result pertains to the effects of hormone replacement therapies under intact endocrine feedback regulation. Further, for the case of a single-tone harmonic positive exogenous signal, bistability and quasiperiodicity arise in the system. The convergence to either of the stationary solutions in a bistable regime is shown to be controlled by the phase of the exogenous signal thus relating this transition to the phenomenon of jet lag.

Introduction
Oscillating nonlinear dynamical systems are standard mathematical models in life science capturing periodic patterns in living organisms [11]. Relevant examples are presented by models of biological clocks that are instrumental in timing of all basic biological processes, see, e.g., [29].
Periodicity arises due to natural (endogenous) phenomena within the system but is also affected by signals (cues) from the environment. When the exogenous signal is periodic, the so-called entrainment of the endogenous (unforced) dynamics can arise [33]. Entrainment, also known as a frequency (phase) locking, is a kind of synchronization which occurs in dynamical systems under external force (for which reason it is sometimes also referred to as a forced synchronization). Even though, in most cases, entrainment of periodical solutions is considered, there is also a more inclusive interpretation of the phenomenon, when nonperiodic endogenous solutions are entrained to the periodicity of the exogenous signal [33]. For instance, in the present paper, when an external periodic forcing is applied in the regime of periodic oscillations, the self-sustained oscillator displays regions of two-mode quasiperiodic dynamics interrupted by a dense set of resonance zones, where the internally generated periodic oscillations synchronize with the external forcing. For further details about entrainment and an insightful discussion of the terminology in synchronization theory, see [27].
A feedback mechanism is necessary for creating a self-sustained oscillation. An early and general mathematical construct to describe a simple periodical biological system is Goodwin's oscillator [12]. It was intended to portray the oscillations in a single gene that suppress themselves via the production of intermediate enzymes. From a control perspective, Goodwin's oscillator is just a third-order linear continuous-time system with a static nonlinear feedback parameterized by a Hill function. Already in this early model, two important properties shared by many mathematical models of biological oscillators have been heralded: One of them is the feedback nonlinearity exhibiting bilateral saturation and another one is the cascade (chain) structure of the linear part. The necessity of saturation is in the boundedness of the involved quantities while the cascade structure is ubiquitous in biochemistry and biology.
The original paradigm of Goodwin's oscillator fits well into the simplified structure of testosterone (Te) regulation in the male [15], where gonadotropinreleasing hormone (GnRH) produced in the hypothala-mus stimulates the production of luteinizing hormone (LH) in hypophysis, which hormone, in its turn, stimulates the production of Te in the testes. The concentration of Te exhorts negative feedback on the concentration of GnRH by inhibiting its release. Goodwin's oscillator is often called the Smith model [30] in the context of endocrine regulation.
Being a conceptual (phenomenological) model, Goodwin's oscillator, in its classical form, neither necessarily fits experimental data nor captures in detail the underlying biological mechanisms. In the endocrine regulation of Te, a significant modeling difficulty is presented by the fact that GnRH secretion by the hypothalamic neurons is not continuous but rather episodic. In fact, synchronized GnRH neurons collectively produce bursts of hormone concentration [16], whose amplitude and frequency are dependent on the concentration of Te. This pulse-modulated mechanism has been established experimentally [35] and implements a negative feedback as the amplitude and frequency of the GnRH pulses decrease with increasing Te levels; see [20] for experimental data.
To bring Goodwin's oscillator (the Smith model) in agreement with the compelling biological evidence, the original static nonlinear feedback of it is substituted with a frequency-amplitude pulse modulation mechanism in [3]. The resulting model is termed the impulsive Goodwin's oscillator. It possesses hybrid dynamics as the feedback is implemented by pulse modulation of the first kind [10] and thus introduces a first-order discrete subsystem into the closed loop of the oscillator.
Besides the hypothalamic-gonadal axes in the male and the female, endocrine loops with pulsatile secretion are found in, e.g., the hypothalamic-pituitary-adrenal axis [34] and growth hormone regulation [32]. Therefore, the dynamics of the latter can be mathematically described in a similar manner.
The most prominent property of the impulsive Goodwin's oscillator is the lack of equilibria that, together with boundedness of the solutions [3], agree well with the original biological function of producing oscillative temporal patterns. This is in contrast with what is experienced in the classical continuoustime version of the mathematical model. A diversity of signal shapes (hormone concentration profiles) is achieved through richness of the dynamics. Even for the impulsive Goodwin's oscillator without continuous time delay, high-periodic solutions, as well as deterministic chaos, are observed [37].
Hormone Te therapy is recommended for men who have both low levels of testosterone in the blood (less than 300 ng/dl) and show symptoms of low testosterone. Exogenous Te can be administered in several ways: injection, patch, transdermal gels, implantable deposits, buccal tablets, etc [24]. Mathematical analysis of what happens when exogenous Te interacts with the pulse-modulated feedback mechanism of the endogenous Te regulation has not been performed before. The different ways of drug delivery require distinct mathematical models. The focus of the present study is on continuous exogenous Te influx that can be achieved by, e.g., hormone patches. Injections are most properly modeled by impulses and thus contribute to the discrete (pulse-modulated) part of the model. This case is left to future research.
Concentration of Te in the male exhibits a circadian rhythm [13] with a period of approximately 24 h, typically modeled in chronobiology as a cosine wave, with a peak value between 7:00 and 7:30 am, [8]. Longer cycles of plasma testosterone levels with periods ranging between 8 and 30 days, with a cluster of periods around 20-22 days have been reported in [9]. How the circadian rhythm entrains endocrine regulation is not well understood so far. Mathematical constructs describing the effect of the circadian rhythm usually depict it as a periodical additive or multiplicative exogenous signal, e.g., [18]. Another approach is to implicitly induce a circadian rhythm in the model solutions by a certain choice of the parameters [31]. Somewhat related to the topic of the present paper, synchronization of (classical) Goodwin's oscillators has been recently treated in [28].
Actual oscillative biological data are never periodic in a mathematical sense. This allows for two interpretations of measured data with respect to the underlying dynamics: One is to see the data as a periodical solution that is disturbed by random impacting signals while another is to attribute it to a chaotic or quasiperiodic attractor. In both cases, entrainment phenomena are highly relevant and have not been studied previously in hybrid oscillators.
The main contribution of this work are as follows: -The equations of the impulsive Goodwin's oscillator are generalized to allow for an exogenous hormone influx governed by known linear dynamics. -Complex dynamical behaviors arising in two important special cases of the exogenous signal being constant and a positive sine wave are studied by means of bifurcation analysis. The former portrays a hormone replacement therapy with Te patches, while the latter describes the effects of circadian rhythm on Te regulation. -It is demonstrated that, for constant exogenous Te, no new dynamical model behaviors arise, compared to the autonomous case. Yet, due to a restricted depth of modulation, the diversity of behaviors is reduced and so are the frequency and amplitude of the GnRH pulses. -Entrainment of the autonomous periodic oscillations of the impulsive Goodwin's oscillator to a sine wave is observed for some combinations of the model parameters, while quasiperiodic and chaotic solutions appear often. -Bistability is discovered in the forced model in contrast to the autonomous one. When in a bistable mode, the convergence of solutions to either of the coexisting attractors can be controlled by the phase of the exogenous sine wave signal.
The obtained mathematical results have bearing on essential biological and medical aspects of testosterone regulation in the male. The model analysis shows that no increase in the cumulative Te level can be achieved through feeding a constant level of Te into the impulsive Goodwin's oscillator unless the pulse-modulated feedback of the model is saturated. The phase of the periodic exogenous signal in the bistable mode of the forced impulsive Goodwin's oscillator can be interpreted as the local time difference in a long-haul longitudinal flight thus relating this effect to jet lag [2]. Indeed, the coexisting periodic solutions have distinctly different mean values of Te and can provide an explanation to the observed endocrine symptoms.
The rest of the paper is organized as follows. First, the mathematical model of the impulsive Goodwin's oscillator is revisited and extended with the dynamics of an exogenous signal. Then a Poincaré map is given for the extended model. Finally, a detailed account of nonlinear dynamical phenomena in the model in hand is provided by bifurcation analysis.

Mathematical models
In this section, the equations governing the impulsive Goodwin's oscillator specialized to Te regulation in the male are summarized without derivation for further use in the bifurcation analysis part of the paper. The exposition in this section generally follows that of [21] but relies on a Sylvester matrix equation instead of convolution integrals in obtaning a pointwise map of the model dynamics.

Pulse-modulated autonomous model
The mathematical model of non-basal Te regulation introduced in [3] is formulated aṡ The continuous state variables of the model are the concentrations of the hormones involved in the regulation of Te. In particular, x 1 ∈ R + is the concentration of GnRH, x 2 ∈ R + is the concentration of LH, and x 3 ∈ R + is the concentration of Te. This implies positivity of the states, also reflected by the matrix A being Metzler. From the biochemistry of the system, the constants b i > 0, i = 1, 2, 3 reflect the half-life times of the involved hormones and g i , i = 1, 2 describe how the production of one hormone is stimulated by another one in the model.
At the discrete times t k ≥ 0, k = 1, 2, . . ., the state vector x(t) undergoes instantaneous jumps. The timing and weights of the impulses producing the jumps are defined by where the superscripts "−" and "+" denote the leftand right-side limits, respectively. In theory of pulsemodulated systems, see, e.g., [10], (·) is called the frequency modulation function, and F(·) is called the amplitude modulation function. The Te concentration is the argument of the modulation functions in (3), i.e., C = 0 0 1 . From the biology of the underlying system, the modulation functions are continuous, monotonous, strictly positive, and bounded from above and below. The function (·) is monotonically increasing, while F(·) is monotonically decreasing. This captures the experimentally observed fact that the hypothalamus responds with sparser GnRH impulses of lower amplitudes to elevated concentrations of Te.
The frequency modulation law in (3) introduces first-order discrete dynamics in the feedback of the impulsive Goodwin's oscillator. Due to the hybrid dynamics, a solution of closed-loop system (1), (2), (3) is defined by the initial conditions x(t 0 ), t 0 , for the continuous states and the discrete one, correspondingly.

Pulse-modulated continuously forced model
An extension to the model above with an exogenous input of Te is introduced in [21] and pursuits two practically motivated goals. One of them is incorporating the basal secretion of Te that is a slowly varying hormone level. It is not related to the pulses of LH but subject to circadian rhythm. The second goal is to describe Te replacement therapies that are administered through, e.g., hormone patches, gel, or implantable Te deposits. These treatments result in a continuous influx of exogenous Te in the closed-loop endocrine regulation system and can be captured by the following model where β(t) ∈ R + is the continuous exogenous input, while t k , as well as λ k , are still given by (3), and D = 0 0 1 . As shown in [21], any solution x(t) to (4)-(5) can be written as In the equations above, the firing times t k and the weights λ k are given by Without loss of generality, assume t 0 = 0.
With respect to the problem of incorporating basal secretion in model (1), (2), the term x p (t) describes the non-basal Te concentration component produced through the mechanism of the pulse-modulated feedback, whereas x f (t) is the Te concentration contributed by the basal level β(t). When a hormone therapy is modeled, x f (t) represents the concentration of Te introduced to the closed-loop system by the drug influx β(t).

Pointwise maps
In the interval t k t < t k+1 , the solution to the system (6), (7) is given by with Substituting (10) into (9) yields For t = t k+1 , the solution (11) has the form with In this way, the evolution of continuous-time system (6), (7) through the jump points t k is described by a discrete map [21] where Map (13) is therefore central in describing the neurally implemented control action of the hypothalamus on the continuous function of the endocrine axis. It is also instrumental in the analysis of nonlinear phenomena arising in the neuroendocrine control loop in an autonomous mode but as well when forced by an exogenous signal.
Constant exogenous Te: Consider first the simplest case, when the exogenous signal β(t) is constant, i.e., Since the solution to (8) with Hence, choosing the initial conditions at the stationary point to avoid (stable) transients as In this way, map (13) does not depend on t k and becomes time-invariant (or autonomous), i.e., with It has exactly the same form as the one in (13), with the only difference that the modulation functions have been shifted by a constant proportional to the exogenous signal β. This brings about an important insight regarding the modeling of basal Te level: It can be readily taken into account in autonomous model (1), (2) by properly modifying the modulation functions (·) and F(·) in (3). It also implies that no new complex nonlinear dynamics phenomena can arise due to constant exogenous signals in the impulsive Goodwin's oscillator as the modulation depth of F(·) and (·) (see [10]) is reduced. This is confirmed by the bifurcation analysis in Sect. 6.1.
Periodic exogenous Te: To capture the basal secretion variations in Te due to the circadian rhythm, consider a periodic exogenous signal of the form where N M > 0, implying that β(t) 0 for all t. In this case, the solution to (8) where the first term asymptotically vanishes. Hence, with the initial conditions becomes periodic The corresponding map then takes the form of (13) with: Note that the effect of a periodic β(t) on the pointwise map given by (13) is time-varying and depends on how the modulation instants t k are distributed on the continuous-time axis.
Similarly, approximative expressions of arbitrary accuracy for a general continuous periodic signal β(t) can be derived by retaining significant terms of the Fourier series.

Time-invariant pointwise map
As shown in [3], a solution x(t) to (1)-(2) satisfies at the jump times the following Poincaré map where The pointwise map in (13) is of the same order as that for autonomous model (18). However, a price to pay is that the discrete dynamics become time-varying. Therefore, the previously obtained results on the impulsive Goodwin's oscillator assuming a time-invariant pulse modulation algorithm cannot be directly applied to the case at hand. The distinction between timevarying and time-invariant modulation has also strong impact on observer design for the impulsive Goodwin's oscillator, as treated in [36].
An alternative approach preserving time invariance of the pointwise map is to augment the continuous state vector of the system thus also increasing the dimension of the map.
For instance, in the case of a shifted sinusoidal exogenous signal, one can augment the state vector with three auxiliary states as followṡ withx(t) = x(t) w(t) , the system model becomeṡ wherē Notice that the modulation variable is continuous despite the state augmentation asCB = 0. Hence, all the formulae developed in deriving the Poincaré maps for autonomous model (1), (2) still apply and the pointwise map is given bȳ As in the previous approach, the above reasoning could easily be generalized to more complicated periodic forms by extending the state vector with auxiliary states corresponding to significant terms in the Fourier series of the exogenous signal. Hence, in this way, a time-invariant map can be created even when the exogenous Te signal is time-varying, with the drawback that the state vector has to be enlarged. A complication arises though in stability analysis of (20) since the matrix A f describing sustained oscillations is not Hurwitz. As stability of the exogenous states w(t) is not of concern in the impulsive Goodwin's oscillator, this issue remains purely superficial.
By exploiting the block-triangular form of the matrix A, the Poincaré map for the augmented system can be represented in terms of the plant and exogenous signal state vectors.

Proposition 1 (Sylvester equation) For a given constant T > 0 and matrices F
is solved by X T = T 0 e F 1 τ Re F 2 τ dτ . The solution is unique when −F 1 and F 2 do not have common eigenvalues.
Proof The result is readily obtained by integrating both sides of The uniqueness condition follows by vectorization of (22) and the uniqueness of the solution to the corresponding linear algebraic system.
The utility of the proposition above is that it reduces the evaluation of a convolution of matrix exponentials to solving a linear matrix equation. (20), map (21) takes the form

Proposition 2 With respect to system model
where W T k is the unique solution of and Proof By direct calculation, one has Further, making use of Proposition 1, it follows that Since A is Hurwitz and A f has only eigenvalues on the imaginary axis in the complex plain as it describes sustained oscillations, W T is the unique solution to (24). It is also observed that the argument in (20) becomes while Cx (t) = (C x(t)). Substituting (25) and (26) into (20) concludes the proof.
It is worth noticing that the feedback firing intervals T k are not known beforehand unless the solutions is periodic. The Sylvester equation in (24) has to be generally solved at each iteration of the map. This also applies to the calculation of w(t − k ) = w(t + k ) = w(t k ). The latter equality is due to the continuity of w(t). However, since w is independent of x, the function w(t) can be calculated in advance for an arbitrarily long time interval and then sampled at the firing times t k obtained from the values of x and (·).

Periodic solutions
Endocrine regulation systems are generally believed to exhibit periodic [22], quasiperiodic [23] or chaotic behaviors [5]. All of those have been discovered in the impulsive model of Te regulation, see [37] and [38], by means of bifurcation analysis of pointwise map (19) and similar.
Define an m-iteration of the map Q(·) as A cycle of Q(·) is a (finite) set of points [17] x 0 , Q(x 0 ), Q (2) where m is the period of the cycle. Then, if, e.g., Q(x 0 ) = x 0 is fulfilled, the system exhibits a 1-cycle. Each point of this set is a fixed point of the m-th iterated map Consider an m-cycle of (13) characterized by the set of fixed points x p,0 , . . . , x p,m−1 with the initial condition on the exogenous system w(0) = w 0 so that Q (14) and (15), it follows that, The period of the m-cycle is then and The relationship between the period of β and the period of the forced m-cycle is provided by (27), (28). The period T of a periodic solution to the dynamical systems (4) and (6) is a multiple of the period of β(t).

Stability of periodic solutions
As argued in [3], Lyapunov stability is not relevant in impulsive systems (1)- (2). This obviously also holds for (4)- (5). The notion of input-to-state stability [14] is often suitable instead when the effects of an exoge-nous signal on the impulsive system behaviors are to be characterized. However, for the system in hand, the dynamics of the exogenous signal are assumed known and, therefore, an input-output framework is not plausible. Since (4)-(5) does not possess equilibria [3], only periodic, chaotic, and quasiperiodic solutions can arise. Notice that lack of equilibria is perfectly in line with the biology underlying the model because only oscillative behaviors of the latter can sustain the endocrine function. One can also argue that an (mathematically) periodic solution to the model in hand does not make more biological sense than, say, a quasiperiodic or chaotic one. The purpose of the pulsatile feedback in non-basal Te regulation captured by model (1), (2) is to produce a certain number of GnRH pulses over a given time interval, e.g., a day. Then the exact timing of the pulses actually does not matter that much and great variability in GnRH pulse frequency and amplitude is observed in biological data.
By making use of the results in, e.g., [1,4,17,26], it can be concluded that a fixed point defining a periodical solution is (locally) asymptotically stable if the Jacobian of the corresponding pointwise map is Schur stable, i.e., all the eigenvalues (multipliers) of it lie inside the unit circle. For instance, as considered in [3], an m-cycle of autonomous system (1)-(3) defined by the fixed point x 0 is asymptotically stable if the Jacobian of Q (m) (x 0 ) is Schur stable.
Since the map Q p (ξ, θ ) in (16) with constant exogenous Te and the mapQ(ξ ) in (21) have the same form as Q(ξ ), the analysis for constant and periodic exogenous Te can be handled in the same way as for the autonomous case. The Jacobian ofQ(ξ ) is readily evaluated as where the matrix exponential is given by (25). Notice that, due to the block-triangular structure ofĀ, the eigenvalue spectrum of J k also includes the eigenvalues of e A f T k that all lie on the unit circle and characterize the marginal stability of the linear periodic dynamics of w(t).
On the other hand, the stability condition can be obtained in terms of the Jacobian given by Unlike the previous case (29), this matrix is of low dimension (the same as A). Now, an m-cycle is asymp- i.e., is Schur stable.

Bifurcation analysis and simulation
To illustrate the possible dynamics arising in (4)-(5) due to exogenous positive harmonic signal (17) with the least period T * , assume M = N , i.e., where For simplicity of the index notation, rename the components of the continuous state vector and let Then the pointwise map Q(·) given by (13) can be rewritten as with Here Following [3], the modulation functions of the intrinsic pulsatile feedback are chosen as

Constant exogenous Te
In this section, the exogenous signal is assumed constant β(t) = const = M, 0 ≤ M ≤ 0.5 and the remaining parameters are 0.25 < b 1 < 0.65, b 2 = 0.014, b 3 = 0.15, g 1 = 0.6, g 2 = 1.5, k 1 = 50, k 2 = 100,  Fig. 1b, the system enters a 5-cycle window 5 through a saddlenode bifurcation at M 0 with a subsequent cascade of period-doubling bifurcations. Note that depending on parameters, the cascade of period-doubling bifurcations may be complete or not. Figure 1c illustrates an example of a finite sequence of period-doubling bifur-cations for p = 3 and Fig. 1d presents an example of a single period-doubling bifurcation. This transition takes place for k 2 = 220 and k 3 = 1.5. The curves shown in green in Fig. 1b-d represent saddle cycles.
Recall that the case of a constant exogenous signal is motivated by Te replacement therapies resulting in a steady influx of the hormone into the feedback endocrine system. The bifurcation diagram in Fig. 1b clearly shows that the Te concentration pulse ampli-tude monotonously decreases with the increasing magnitude of the exogenous signal. The exogenous signal also reduces the oscillativity of the solutions. Indeed, the character of the solutions changes from chaotic to a 1-cycle within the range of the bifurcation parameter. This brings about a somewhat unexpected conclusion that a Te replacement treatment can actually lead to lower cumulative levels of the hormone in the blood due to the action of the impulsive endocrine feedback. Actually, there is experimental biological evidence that collaborates the simulation result. According to [25] shorter treatment with anabolic steroids mixture (with similar effect to Te) decreased spontaneous physical activity in male mice.
Naturally, in the modeled scenario, the hypothalamic regulation of Te is assumed to be intact. The analyzed model does not capture the effect of Te on the organism and thus cannot be used for reasoning about the symptoms associated with Te deficiency.

Periodic exogenous Te
By contrast to the case of a constant exogenous signal, where the dynamics are defined only by features of the impulsive self-oscillatory system, a periodic (continuous) exogenous signal introduces an interplay between high-frequency (ultradian) self-oscillations of the autonomous pulse-modulated system and an exogenous periodic forcing signal that is assumed to be of low frequency (e.g., circadian). This interaction between the impulsive oscillator and the continuous one gives birth to a large variety of nonlinear dynamical phenomena, including quasiperiodicity, entrainment (modelocking), multistability, and chaotic dynamics.
In order to examine some of the basic aspects of these behaviors and transitions between them, the following parameter values are chosen: k 2 = 220, k 3 = 1.5, p = 4. Since the exogenous signal β represents the effect of circadian rhythm, the period T * = 1440 min = 24 h, and 0.2 < M < 2.0, (M = N ), 0.23 < b 1 < 0.69 are used. Other parameters are left as in the case of constant exogenous Te influx. Figure 2a provides an overview of the bifurcation structure that can be observed in the (b 1 , M)-parameter plane. Recall that M is the forcing amplitude in (31) while b 1 is the clearing rate constant of GnRH.
As one can see in Fig. 2a, the bifurcation diagram is characterized by a dense set of entrainment windows [27]. Between these windows, there are parameter combinations leading to quasiperiodic and chaotic dynamics (the corresponding regions are shown in white).
To understand the dynamics of map (32), consider the transitions that occur when the forcing amplitude M is changed along the scan B in Fig. 2a for b 1 =  0.42104. Figure 2b shows the dependence of the wave number W ϕ (see [6,7]) for the ϕ variable of map (32) on the parameter M for b 1 = 0.42104. The diagram displays the successively occurring regions of periodic behavior (entrainment regions) and aperiodic behavior.
Recall that the wave numbers, in the considered case, behave similar to the well-known rotation numbers. When the map exhibits an m-cycle and, during m iterations, the phase ϕ makes n rotations, the wave number of this cycle is W ϕ = n:m. Within an entrainment region, the wave number is constant. Several of the most prominent entrainment regions are marked in Fig. 2b with the corresponding wave numbers (1 : 14, 4 : 27, 3 : 20, 2 : 13, 5 : 32, 3 : 19, 1 : 6, and 4 : 25).
Depending on the value of M, the dynamics outside the entrainment region may be chaotic or quasiperiodic. As usual, information concerning the transitions from entrainment to quasiperiodicity or chaos and vice versa can be obtained by following the variation of the Lyapunov exponents. Figure 2c depicts the two largest Lyapunov exponents 1,2 as functions of the forcing amplitude M. The largest Lyapunov exponent 1 is positive in most of the considered interval 0.2 < M < 0.84 indicating chaotic dynamics. Inside the entrainment regions, 1 is negative. The values of 2 are negative everywhere, which signifies that no hyper-chaotic dynamics occur.
Numerical experiments show that, for small values of the forcing amplitude M, the dynamics of map (32) are either quasiperiodic with irrational wave numbers or periodic (entrained), when the wave number is rational. In both irrational and rational wave number cases, there exists a stable closed invariant curve.
When the wave number is irrational, the invariant curve is densely filled with points of quasiperiodic trajectories. When the wave number is rational, the closed invariant curve contains a pair of cycles, one of which is stable, while the other is a saddle. The attracting invariant curve is formed by the saddle-node or saddle-focus 14-cycle and includes the points of the saddle and stable 14-cycles (see Fig. 3c). The variation of the wave number for ϕ with respect to M is shown in Fig. 3d). When the parameter M increases or decreases, the stable and unstable fixed points collide and disappear in a fold (saddle-node) bifurcation at the points M R 0 and M L 0 , respectively (see Fig. 3a). Depending on the direction of the scan, two different bifurcation scenarios can be observed. When passing through the value M L 0 with decreasing M, the fold bifurcation leads to a transition from a stable 14-cycle to quasiperiodic dynamics. As illustrated in Fig. 3e, the largest Lyapunov exponent 1 becomes zero at the fold bifurcation point M L 0 . Outside the entrainment region to the left of the point M L 0 , the dynamics are quasiperiodic with 1 = 0. Another option is that the disappearance of the stable 14-cycle is followed by the appearance of a chaotic attractor. To the right of the point M R 0 , the largest Lyapunov exponent becomes positive, signaling the development of chaotic dynamics (Fig. 3e). The inset in Fig. 3d shows the enlargement of the devil's staircase structure for the wave number to the right of the point M R 0 . Figure 4a presents a magnified part of the diagram in Fig. 3e, that falls to the right of the point M R 0 . As illustrated in Fig. 4b, this region is characterized by a dense set of periodic windows (see also Fig. 3d). Figure 3f shows an example of the quasiperiodic attractor.
For large amplitudes M of the exogenous signal, there are other mechanisms of transition from/to entrainment, and they typically lead to multistability and chaotic dynamics. To illustrate the transition to entrainment in detail, Fig. 5a shows a magnified part of the bifurcation diagram that is outlined by the red rectangle in Fig. 5a. With an increase in M, the system enters a 1 : 6 entrainment region through a saddle-node bifurcation at the point M = M 1 . Note that the 1 : 6 entrainment region overlaps with the 3 : 19 window of periodicity, so that M 1 < M < M 0 is the region of bistability where the stable 19-cycle coexists with the stable 6-cycle (see also Fig. 2b). Here M 0 , M 1 are the saddle-node bifurcation points for 19-cycle and 6-cycle, respectively.
As M increases, the 6-cycle undergoes four perioddoubling bifurcations, two direct bifurcations, and then two reverse ones (Fig. 5b). The bifurcation diagram in Fig. 5c illustrates the direct and reverse perioddoubling bifurcations for 6-cycle and 12-cycle. Here M 1 is the saddle-node bifurcation point in which the stable node (branch 1) and saddle 6-cycles (branch 2) are born. M L 2 , M R 2 are the period-doubling bifurcation points for 6-cycle. M L 3 , M R 3 are the period-doubling bifurcation points for 12-cycle. Figure 5d shows a stable 6-cycle for M = 0.846 and b 1 = 0.42104. Multistability, i.e., coexistence of attractors in the phase space of a dynamical system, is a typical phenomenon in nonlinear dynamics. These attractors may arise through a saddle-node bifurcation and, with changing parameters, they can give rise to an infinite sequence of period-doubling bifurcations, leading to the transition to chaos (see Fig. 5). The latter results in parameter domains wherein, alongside with stable cycles, there are coexisting modes of chaotic oscillations. Under such conditions, an exogenous disturbance, even of low intensity, can cause a transition from one attractor to another. An important multistability property of the pulsemodulated Te regulation model is the following. Apparently, bistability in map (32) can be controlled through the phase of the periodic exogenous signal θ in (31). This is illustrated in Fig. 6 where the model is manipulated to first enter an 19-cycle and then, due to an instantaneous change in θ , transfer to a 6-cycle. As seen in the plot, a reverse transition is also possible, once again by means of controlling θ . Notably, the mean value of the Te concentration is higher for the 19-cycle. Thus, with only one change of the phase θ , the model can move to a stationary solution that corresponds to a lower hormonal activity and Te concentration. In the context of endocrine regulation, circadian disruptions [2] (day light shift) arise due to, e.g., long-distance longitudinal travel (jet lag) and have not been addressed before via mathematical modeling.
A similar problem, without any relation to the endocrine system, is studied in [19] with respect to the forced Kuramoto oscillator representing the neurons of the suprachiasmatic nucleus in the hypothala- Fig. 6 Bistability: Transition from 19-cycle to 6-cycle and back to 19-cycle controlled by the phase of the 24-h periodic exogenous excitation β(t) (lowest plot, in red). Time evolution of the continuous state variables x(t), y(t), z(t) is depicted mus implementing the circadian clock. The circadian system is responsible for tuning of physiological processes to the daily light cycle. The equilibrium points of the continuous model are analyzed to illustrate the different types of system dynamics. On the contrary, the model given by (4)-(5) does not have equilibria, thus exhibiting a completely different dynamical mechanism related to the phase of the exogenous signal.

Conclusions
Dynamical behaviors forced by a continuous exogenous signal in a previously developed pulse-modulated mathematical model of non-basal testosterone (Te) regulation are studied. The exogenous signal can represent, e.g., the influx of a drug used in a hormone replacement therapy, the dynamical effects due to circadian rhythm, or interactions with other endocrine loops of the organism.
Two equivalent ways of calculating the pointwise Poincaré maps that capture the evolution of the continuous states of the model from one firing of the pulsemodulated feedback to the next one are proposed: one making use of the linearity of the continuous part of the model and another via augmentation of the continuous state vector. Bifurcation analysis and simulations reveal intriguing properties of the model solutions that provide insights into experimentally observed biological phenomena. First, administering a constant influx of Te into the system consequently decreases the mean value of the Te concentration. This property arguably explains the adverse effects in Te replacement therapies. Second, a periodic exogenous signal entering the model in a periodical mode results in most cases in a non-periodic forced solution. Thus, even in this simple model, the circadian rhythm exhorts mainly quasiperiodicity. Finally, for periodic solutions of the model entrained by circadian rhythm, bistability is discovered. The convergence to either of the coexisting stable stationary solutions producing distinctively dissimilar signal shapes of the Te concentration can be controlled by the phase of the circadian rhythm model. The bistability phenomenon can offer a plausible explanation to endocrine disorders due to jet lag and shift work.