Abstract

The paper deals with the model-based estimation of hormone concentrations that are inaccessible for direct measurement in the blood stream. Previous research demonstrated that the dynamics of nonbasal endocrine regulation can be closely captured by linear continuous models with time delays under a pulse-modulated feedback. The presence of continuous time delays is inevitable in such a model due to transport phenomena and the time necessary for an endocrine gland to produce a certain hormone quantity. Yet, thanks to the finite-dimensional reducibility of the linear time-delay part of the system, a finite-dimensional model can be used to reconstruct both the continuous and discrete states of the hybrid time-delay plant. A hybrid observer exploiting this possibility is suggested and analyzed by means of a discrete impulse-to-impulse mapping.

1. Introduction

Hormones mediate communication between organs and tissues through the bloodstream carrying chemical messages that regulate many aspects in the human body, that is, metabolism, growth as well as the sexual function and the reproductive processes. Hormones are secreted by endocrine glands directly into the bloodstream in continuous (basal) or pulsatile (nonbasal) manner. Endocrine glands, interacting via hormone concentrations in blood, build up dynamical control loops characterized by self-sustained oscillations of the involved physiological quantities [1].

The endocrine system of testosterone regulation in the male essentially consists of three hormones, namely, gonadotropin-releasing hormone (GnRH), luteinizing hormone (LH), and testosterone (Te). GnRH is produced in the hypothalamus of the brain and released in short pulses. Reaching the pituitary gland, GnRH stimulates production of LH, which in turn stimulates production of Te in the testes. Finally, both the GnRH outflow and the LH secretion are subject to feedback inhibition by Te [2]. However, the inhibition of LH has a relatively small effect on the dynamics of the closed-loop system and therefore not considered in this paper.

An impulsive mathematical model of testosterone regulation was proposed in [3] and is shown to comport with experimental data in [4]. It constitutes an impulsive version of Goodwin oscillator, a mathematical model that is well known in mathematical biology (see, e.g., [59]). The impulsive Goodwin oscillator consists of a continuous and an impulsive part [10], thus possessing hybrid dynamics and presenting a special version of an impulsive differential system [1014]. It mathematically portrays the concept of pulsatile hormone regulation described in medical literature (see, e.g., [15]).

More recently, the impulsive Goodwin oscillator was augmented with a time delay in the continuous part of the system [16, 17], making it more aligned with the biological nature, as transport phenomena and biosynthesis are omnipresent in endocrine and metabolic systems [1825]. With the time delay taken into account, the pulse-modulated model of endocrine regulation acquires an infinite-dimensional continuous part. The closed-loop dynamics become therefore both hybrid and infinite-dimensional, and this combination is mathematically challenging and so far rarely treated. However, the cascade structure of the continuous part, together with the impulsive feedback, allow application of the concept of finite-dimensional reducibility (FD-reducibility), [16, 17]. In particular, it was shown [17] that the dynamics of an impulsive time-delay system with an FD-reducible continuous part coincide on certain time intervals with the dynamics of a delay-free impulsive system. This idea plays a key role in the present study.

Concentrations of the hormones secreted in human hypothalamus that is located in the lower central part of the brain are not available for direct measurement due to ethical reasons and need to be estimated. It poses an unusual observation problem. A considerable number of papers is devoted to the observability of hybrid systems, for example, [2628]. The discrete states of a system are usually assumed known, while observers for hybrid systems that are able to reconstruct discrete states from only continuous measurements are not so well covered in the literature.

In endocrine systems with pulsatile secretion, the highest degree of uncertainty is associated with the discrete (impulsive) part whose states have to be reconstructed from hormone concentration measurements. Two model-based estimation approaches are currently known. The first one is based on batch deconvolution techniques (blind system identification) [29, 30], while the relatively recent second one employs a state observer, whose estimates are corrected by output estimation error feedback [31, 32]. An extension of the observer scheme proposed in [31] to impulsive systems with time delay in continuous part was considered in [33]. Unlike the case treated in [33], the observer proposed here does not explicitly involve a delay but is rather based on a finite-dimensional plant model. Hence, the main contribution of the paper is in the novel structure and subsequent analysis of a hybrid observer exploiting a finite-dimensional model to reconstruct the states of the time-delay system.

Notice that impulsive feedback in the observer treated below is not contributed by design to achieve a performance objective but rather constitutes an integral and unmeasurable part of the plant model. On the contrary, in the impulsive observers for state estimation of linear and nonlinear continuous systems proposed in [3437], the observer state is updated in an impulsive fashion in order to achieve, for example, faster convergence. This distinction results in a major complication in observer design for plants with intrinsic impulsive feedback as the timing and weights of the impulses are unknown and have to be estimated by the observer.

A preliminary version of the present material without proofs of the main statements was presented in [38].

The paper is organized as follows. First, an impulsive time-delay model is summarized and reduced to an equivalent delay-free one. Then, making use of the reduced model, a hybrid observer is proposed and a pointwise (impulse-to-impulse) mapping describing its dynamics is derived. Further, the properties of the mapping pertaining to the observer performance are investigated. Then the impulsive time-delay model of testosterone regulation is described. Numerical simulations and calculations illustrating the observer design performance are also provided.

2. System Equations

Consider an impulsive time-delay model [16] given by the equationswhere , , , , are constant matrices and .

In (1), is the scalar controlled output, is the measurable output vector, is the state vector, and is a constant time delay. The amplitude modulation function and frequency modulation function are continuous and bounded: is nonincreasing and is nondecreasing.

System (1) is considered for subject to the initial condition , , where is a continuous initial vector function. The state vector of system (1) experiences jumps at the times , . The condition ensures that the modulating signal is continuous in time.

Only the time-delay values that are less than the minimal distance between two consecutive impulses are considered:so that for all . This condition implies that only one firing of the pulse-modulated feedback in (1) is possible within a time interval whose length is equal to the time-delay value.

Suppose that the linear part of the system possesses the property of finite dimensional (FD) reducibility [16, 17], implying thatThe notion of FD-reducibility is a formalization of the so-called “linear chain trick” originating from [39, 40] for the system in question.

3. Reduction to a Delay-Free Impulsive System

Define the matrices , . Introduce a delay-free impulsive system:

The following lemma obtained in [17] reveals the relationship between the solutions of system (1) and those of system (4).

Lemma 1. Consider solutions , of systems (1), (4), respectively. Assume that and . Then it holds that , and for all . Moreover,At the same time, generally, the solutions do not coincide entirely

The result above will be exploited further in the paper to design a finite-dimensional observer for the infinite-dimensional hybrid system in (1). Note that the value of the time delay in the delay-free impulsive system still influences the system dynamics as affects the matrix coefficients , of (4).

4. A Hybrid Observer

The purpose of state observation in hybrid closed-loop system (1) is to produce estimates of the impulse parameters . Notice that, unlike in the conventionally treated hybrid state estimation problem formulations, the jump times are considered to be unmeasurable in (1) and require estimation. In fact, the problem solved by the proposed observer is synchronization of the firings in the feedback of the plant representing its discrete state and those of the observer.

From Lemma 1, it follows that one can produce estimates of the impulse parameters of (1) by exploiting the delay-free model in (4). To evaluate , an estimate of the continuous state vector of (4), that is, is produced by the hybrid observer:

Notice that , are generally discontinuous in time.

The switched feedback gain is zero in the time intervals where the solutions of system (1) and those of system (2) do not coincide, while the static feedback gain is chosen to satisfy the stability conditions derived in Section 8.

5. Synchronous Mode

Keeping in mind that the purpose of the hybrid estimation here is essentially synchronization, and following [31], introduce the notion of a synchronous mode for the plant-observer system (4), (7). Let be a solution of plant equations (4) with the parameters , , and . Suppose that the plant is already running at the moment when the observer is initiated, that is, , for some integer .

Consider the solution of observer equations (7) subject to the initial conditions , , that yields , , , , and for . Such a solution will be called a synchronous mode with respect to . Thus, a synchronous mode is a null solution of the state estimation error dynamics of hybrid observer (7) on .

Following [31], a synchronous mode will be called locally asymptotically stable if, for any solution of (7) such that the initial estimation errors and are sufficiently small, it follows that and as . The latter implies as . In the definition above, the operator stands for any vector norm.

To ensure practical usefulness of the observer, stability properties of the synchronous mode have to be investigated. By choice of , the synchronous mode has to be rendered asymptotically stable with a suitable convergence rate and domain of attraction.

6. Pointwise Mapping and Its Properties

Consider the pointwise mapping describing the evolution of the observer hybrid state from one firing of the impulsive part in (7) to the next one:

For any integer numbers and , , define the sets

Hence, each point of the observer hybrid state belongs to one of the sets , that is, to each one can uniquely match two points and of the observed system (if , these points coincide) such that , .

Introducewhere . Note that the functions are piecewise continuous due to the definition of .

Define at withFor brevity sake, denote

Theorem 2. Pointwise mapping (8) is given by the equations

Proof. See Appendix A.

Theorem 3. The mapping is continuous.

Proof. See Appendix B.

It will be shown in the next section that the mapping is not continuously differentiable in the whole state space. However, due to its local differentiability, local stability properties of mapping (12) characterizing the dynamics of the observer state can be investigated via linearization.

7. Linearization of the Discrete-Time Map

The behaviors of pointwise mapping (12) in vicinity of the points will be studied with respect to local stability of a synchronous mode.

To show the smoothness of the mapping introduced below at the points , divide each set for all into two subsets and by the hyperplanes , that is, , where

Note that the set () is not divided into subsets and due to the assumption on delay in (2) and thus implying and .

Consider a point for some . It can be seen that the closures of the four sets , , , intersect at the point . Moreover, . For a sufficiently small neighborhood of the point , the mapping can only take the values , where and is one of the four sets , , , ; see Figure 1.

Denote for brevity ,

Theorem 4. If and the scalar functions have continuous derivatives, then the partial derivatives of are continuous at the points , and given by the following expression:

Proof. See Appendix C.

Introduce additional notation referring to mapping (8) and augmenting the continuous state vector with the discrete state. Define the function

Set for . Then by the definition of in Section 6 one has , and from Theorem 2 it follows that , where .

Propagation of the observer hybrid dynamics through the firing times is described by iterations of the operator . The th iteration is defined asTheorem 4 implies that the operator can be linearized in a vicinity of .

From Section 5, it follows that the synchronous mode with respect to is completely characterized by the vector sequencewhere is a number from the definition of a synchronous mode.

Then the Jacobian of at the point is calculated aswhere

By the chain rule, it follows that, for any , the Jacobian of the th iteration of the mapping is given by the expression

8. Local Stability of a Synchronous Mode with Respect to an -Cycle

A solution of (4) is called -cycle if it is periodic with exactly pulse modulation instants in the least period. The existence conditions of an -cycle in pulse-modulated time-delay system (4) were studied in [16, 17].

Let be an -cycle of plant (4), where is some integer, . The existence conditions of an -cycle of pulse-modulated system with time delay are readily derived in [17]. Then , , . Consider a synchronous mode of observer (7) with respect to and let be the corresponding vector sequence as in (17), such that is satisfied.

Consider previously defined matrices . Since , the sequence contains at most distinct matrices, namely, . The theorem below provides a simple tool for checking local stability of observer (7).

Theorem 5. Let the matrix product be Schur stable; that is, all the eigenvalues of this matrix lie strictly inside the unit circle. Then the synchronous mode with respect to is locally asymptotically stable.

Proof. The result can be proved along the lines of Theorem  3 in [31].

Theorem 5 formulates a stability condition guiding the choice of the observer gain that appears in the matrix . As pointed out above, the condition is local and depends not only on the coefficients of the system, but also on the parameters of the observed periodic mode. In particular, the multiplicity of the periodical solution in the plant has to be known. The spectral radius of Jacobian (18) reflects the local convergence rate of the linearized observer dynamics. To optimize the observer performance, the static gain can be chosen numerically to fulfill the conditions of Theorem 5 while minimizing the spectral radius of the Jacobian.

9. Mathematical Model of Testosterone Regulation

To model testosterone regulation in the human male [17], the case of a third-order system (1) with the matricesis considered. Here , , , , are given positive parameters reflecting the kinetics of the involved hormones. From the biology of the system, one has for . The elements of correspond to the concentrations of GnRH (), LH (), and Te ().

The presence of a constant time delay in closed loop relates to a delay in the hormone transport in the blood stream and a delay occurring in hormone synthesis prior to secretion [41]. The contribution of the transport delays is relatively smaller than that due to synthesis of testosterone. In the simulations, the delay value is selected so that the minimal distance between two consecutive impulses does not exceed the sum of the testosterone synthesis and hormone transport delays, that is, . This is in line with the data provided in [18, 20].

The concentrations of Te and LH can be measured in the blood, while the concentration of GnRH is typically not available in humans and has to be estimated. Nonetheless, the level of testosterone is usually overly more noisy than the level of LH (see, e.g., [4]), and it is difficult to distinguish between the basal and pulsatile components. Thus the structure of the output row vector is chosen so that only the measurement of LH concentration is taken into account.

Within a feedback construct, pulsatile secretion of a hormone gives rise to a dynamic system where amplitude and frequency modulation are employed to control the concentrations of other hormones, ostensibly in order to induce sustained oscillations in the closed-loop system.

As the amplitude modulation function and frequency modulation function Hill functions with the following (continuous) parameterizations are chosenwhere , , , , , are positive parameters and is integer. It is easy to check that the functions and are smooth, strictly monotonic and bounded.

Since , system (1) with matrix coeffecients (21) is FD-reducible, and, hence, can be represented in the form of delay-free system (4) with , .

The matrix exponentials are given bywhereIntroduce the numbersThen it can be easily seen that

The state vector of system (1) with matrix coeffecients (21) experiences jumps at the times , portraying nonbasal (episodic) release of GnRH. However, because of the matrix relationship , the assumptions of Theorem 4 are valid and the impulse-to-impulse mapping is smooth.

Below the observer design and performance are exemplified by two cases of periodical solutions in the plant arising for different values of the time delay within the considered interval. Notice that, for the numerical values in question, the multiplicity of the periodical solutions in the plant decreases with increasing delay.

10. Numerical Examples

Assume the following values in model (1): , , , , , andSince , then the time-delay value is within , to make the analysis of this paper applicable.

10.1. Observation of a 4-Cycle

Let . Then, the plant has a stable 4-cycle withwhere denotes transpose.

Choose the observer feedback gain in the formwith , , . ThenHence, the characteristic polynomial of is independent of and equal toSince , , , are positive and , are nonnegative, is Hurwitz stable.

To ensure the (locally) fastest convergence rate, find , for which the synchronous mode is locally asymptotically stable and the spectral radius of is minimal. By inspection of Figures 2, 3, and 4, such values of , are

Hybrid observer performance can be measured in numerous ways. The convergence to a synchronous mode is characterized here by the first time instant when comes into -neighborhood of and never leaves it:This criterion somehow captures the most demanding state estimation error in the hybrid observer since all the information regarding the discrete state in (1) comes from the continuous measurements. The relationship between the value of the threshold in (33) and is depicted in Figure 5.

10.2. Observation of a 2-Cycle

Increasing the time delay to yields a stable 2-cycle with

A search for , that render a locally asymptotically stable synchronous mode and minimal spectral radius of gives (see Figures 6, 7, and 8):

The relationship between the value of the threshold in (33) and for a certain stabilizing observer gain is depicted in Figure 9.

11. Conclusions

A state estimation problem motivated by unmeasured hormone concentrations in the system of nonbasal testosterone regulation in the human male is considered. The system dynamics are modeled by a linear continuous time-delay system under intrinsic pulse-modulated feedback. The continuous part of the model is known to possess the property of finite-dimensional reducibility that opens up for the use of a finite-dimensional (delay-free) model for the reconstruction of the discrete and continuous states of the process model. A hybrid observer exploiting this possibility is introduced and analyzed by means of a discrete impulse-to-impulse mapping.

Appendices

A. Proof of Theorem 2

Consider the difference in the interval and suppose thatfor some and , such that (Figure 10). Obviously, satisfies a differential equationwhereat all the points , where has no jumps.

Derive explicit formulas for the map (8). Introduce a number such that .

For (i.e., ) the vector function has no jumps in the interval . Hence,implies (12) for .

For the function has jumps at the points , . The assumption guarantees that for , while the point may lie either in the interval or in the interval , so these two cases should be considered separately.

Suppose that (i.e., ). Then .

Find .(i)If , one has(ii)If , one has Hence .

Thus,implying (12) for .

Suppose that . Prove that , wherefor some , such that .(i)If , then(ii)If , then

For one has

Thus one concludes thatwhere

Then (A.12) can be rewritten as

Sinceequality (A.14) implies (12).

B. Proof of Theorem 3

Sinceit is straightforward to see thatfor , . Since the functions are continuous for all , then the function can have gaps only on the surfaces in space , where either or for some .

Yet, from (B.2), (B.3) and because and , it follows thatand the function is continuous everywhere.

C. Proof of Theorem 4

Sincehas a gap on the surface , then is continuous everywhere on , except for (if then is continuous everywhere on ).

It is easy to see that for any there exists a sufficiently small neighborhood of a point such that . Hence, the partial derivatives of can have gaps either on the surface or on the surface .

Let . From (B.2) it follows that and . Hence and have no gaps on .

Let . From (B.3) it follows thatHence and have no gaps on .

Formula (14) follows by direct calculations along the lines of Theorem  2 in [31].

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

A. Medvedev and D. Yamalova were in part financed by the European Research Council, Advanced Grant 247035 (SysTEAM) and Grant 2012-3153 from the Swedish Research Council. A. Churilov was partly financed by the Russian Foundation for Basic Research, Grant 14-01-00107-a. D. Yamalova and A. Churilov acknowledge Saint Petersburg State University for research grant 6.38.230.2015.