Global Stability in Dynamical Systems with Multiple Feedback Mechanisms

A class of n-dimensional ODEs with up to n feedbacks from the n’th variable is analysed. The feedbacks are represented by non-specific, bounded, non-negative C1 functions. The main result is the formulation and proof of an easily applicable criterion for existence of a globally stable fixed point of the system. The proof relies on the contraction mapping theorem. Applications of this type of systems are numerous in biology, e.g., models of the hypothalamic-pituitary-adrenal axis and testosterone secretion. Some results important for modelling are: 1) Existence of an attractive trapping region. This is a bounded set with non-negative elements where solutions cannot escape. All solutions are shown to converge to a “minimal” trapping region. 2) At least one fixed point exists. 3) Sufficient criteria for a unique fixed point are formulated. One case where this is fulfilled is when the feedbacks are negative.


Outline
First, an n dimensional system with feedbacks from the n'th variable is introduced and some applications from bio-medicine and biochemistry are described.Then, analysis of a scaled version of the system is made including fixed point investigation.Finally, an easy applicable sufficient criterion for a unique, globally stable fixed point is formulated and proved.
Mathematically, the results in this paper follow from the dimensionless form of the equations stated in (6) of Section 2. But before turning to this form we motivate and discuss the dimensional form of the equations in Section 1 as we relate the system to applications and earlier results.

Introduction
Many applications of ODEs to physics, chemistry, biology, medicine, and life sciences give rise to non-linear non-negative compartment systems.These include metabolic pathways, membrane transports, pharmacodynamics, epidemiology, ecology, cellular control processes, enzyme synthesis, and control circuits in biochemical pathways [1]- [9].
This paper concerns the stability of the solutions of such models.More specifically, the paper presents criteria for both local and global stability of all systems of ODEs that can be presented as a compartment model with n compartments, on the form shown in Figure 1.Here the n'th variable may have a non-linear feedback on any of the variables.The main results of this paper are: • Existence of a "trapping region"-a compact set with non negative elements in which any solution will be trapped after finite time.• At least one fixed point exists and a real valued function of one variable and the system parameters determines the fixed point.• A unique, globally stable fixed point exists if the norm of a real valued function of one variable and the system parameters is less than 1.

Motivating Background
Figure 1 reflects typical hormone regulation.Since a hormone has to bind to a receptor to cause a feedback, a bounded number of receptors justify that the feedback functions i f are bounded.Examples of systems corresponding to Figure 1 with 3 n = are models of the hypothalamic-pituitary-adrenal axis (HPA axis) concerning the interplay of three hormones in the human body [1]- [5] [10] [11].Here cortisol exerts a feedback on two other hormones that are involved in the production of cortisol.The system is related to stress and depression.Also testosterone secretion has been modelled by a three dimensional compartment ODE-model including a single feedback [12] which is included in the system investigated here.Similar models exist of gonadotropin hormone secretion [13], for describing female fertility [14]- [16] and for cellular metabolism [17].
A two dimensional model of the HPA axis corresponding to Figure 1 is found in [18].Here the focus is on a sufficient criterion for a locally stable fixed point.However it is made clear that a global investigation is preferable.Criteria for global stability of solutions are rare.An example is through use of a Liapunov function [6] [12] that can be employed to some problems.Existence and construction of a Liapunov function are unfortunately not easily addressed in general, and Liapunov functions are not used in this article.
Some general and analytical considerations partly similar to our has been considered in previous papers [8] [19].However, [8] investigate only a feedback from compartment n to compartment 1.The approach of [20] proves the existence of periodic solutions but does not touch upon global stability.
The mathematical results derived in this article relate to the robustness of hormonal systems, cellular metabolism, etc.The existence of a trapping region ensures that non negative initial (hormone) values lead to (hormone) levels that stay non negative and bounded which is reasonable.Existence of locally stable fixed points may be interpreted as states where (hormone) levels may settle.Perturbing parameters such that a solution enters the basin of attraction to another fixed point may then be interpreted as a new (physiological) state (for a person).Or distinct stable fixed points may be interpreted as states for distinct groups (of people).In case of a unique, globally stable fixed point the long term behaviour is very robust to perturbations.

Mathematical Formulation
We consider an n dimensional system of differential equations with n non negative variables i x , x may exert a feedback on all the variables thus making the system non-linear.( ) with production rates 0 i k > and consumption rates 0 i w > .The feedback from n x on i x occurs through the function

( )
i n f x .The following demands are posed for the feedback functions: ( ) . The feedbacks are modelled to influence the positive stimulation of the variable in a compartment but with a saturation which justify why the feedbacks must be bounded functions.This means a feedback acts like an adjustable tap that affects the production of variable i x as a function of n x .When modelling many biochemical systems, such as hormone dynamics, saturation is present due to a finite number of binding sites, e.g., receptors.When all binding sites, or receptors, are occupied and work at maximum speed, then an increase in concentration has insignificant effect.The feedback functions must not attain negative values since this corresponds to reverting the flow.When the concentration of n x is zero the feedback functions must not cause the production rates to be zero.Therefore, ( ) In life sciences the consumption rates correspond to elimination rates in general and are therefore by and large constants.However, some results hold even if we allow the w i 's to be bounded non-negative functions of i x .The models outlined in [1] [3] [10] [18] [12] are covered by Equation ( 1).An example of a typical feedback function is the sigmoidal Hill-function ( ) and α being an integer.Such Hill-functions are often the result of underlying inter cellular enzymatic reactions regulating feedbacks in the quasi-steady-state approximation [21].In neural networks applications a utilized feedback function is the hyperbolic tangent [22].

Analysis
First a scaling is performed to facilitate the analysis.Defining dimensionless variables , i X τ by the equations

{ }
,for 1, , , where 0 t and i d are constants to be defined.

( ) (
) Choosing 0 d as a unit of inverse time we get [ ] and A scaling of Equation (1) thus leads to the dimensionless system ( ) ( ) ( ) . τ corresponds to a dimensionless time and the value of 0 d may be fixed arbitrarily.Differentiation with respect to τ will be noted by a dot such that d d

Existence and Uniqueness of Solutions
C and locally Lipschitz in ( ) local existence and uniqueness of solutions to Equation ( 6) are guaranteed given Since the right hand side of Equation ( 6) in addition fulfils, ( ) , , , , , , at least one global solution exists.Here we have made exclusive use of the fact that the i F 's are bounded.Combined with the aforementioned local uniqueness result a unique global solution exist.Alternatively one may combine the fact that Equation ( 6) is autonomous with theorem 3.22 of [23] to guarantee global existence and uniqueness of solutions to Equation (6).

Positivity of Solutions
Avoiding negative modelling hormone levels is necessary for a sound model and is proved in the following lemma.
Lemma 1.The non negative hypercube is an invariant solution set to Equation (6) Proof.Given a solution initially in the non negative hypercube we consider the behaviour at a boundary of the hypercube-a hyperplane defined by 0 . Considering Equation ( 6) and first considering 1 j = ( ) which is non negative for all non negative 1 , , n X X  . Then, considering which is a product of non negative factors for all non-negative 1 , , n X X  . This means a solution cannot pass a boundary given by the non negative hypercube due to the aforementioned (local) uniqueness property of solutions.

Existence of a Fixed Point
The fixed point condition of Equation ( 6) can be expressed This means that for each fixed point value nss X the fixed point values of the other variables are easily calculated.The equation ( ) may not be explicitly solvable for nss X .However, existence of a solution can be guaranteed and the solution can be numerically approximated.Define the functions ( ) Thus, finding fixed points of Equation ( 6) is equivalent to finding n X that fulfills ( ) ( ) Now choose .
Since L and R are continuous so is φ and notice that ( ) ( ) ( ) < and ( ) ( ) ( ) 0 Then there exists a ] [ . This means there exists at least one fixed point of the system.Notice that any fixed point is in the set

Sufficient Criteria for a Unique Fixed Point
We now discuss a sufficient criterion for existence of a unique fixed point of the system.Let nss X ′ denote the smallest existing fixed point of Equation (11).If ( ) , there can only be one fixed point.Since ( ) = a sufficient criteria for only one fixed point is ( ) If the feedback functions correspond to negative feedbacks or are independent of n X , the criteria is fulfilled .
Since i F only attains non negative values, none positive feedbacks guarantee that there exists exactly one fixed point.

Trapping Region
A trapping region is a set,

( )
W  , where a solution will never escape if it is once in there.It is a physiological desirable property of a model, since this guarantees, that reasonable initial values lead to reasonable levels of the variables for all future time.
Lemma 2. Let and define .
Then ( ) for all non negative values of the remaining variables.This means that  is a 'trapping region' for 1 X .Using this region for 1 X we can find a 'trapping region' for 2 X and so on by induction.Assume ( ) Then for ( ) This ensures that ( ) . This means there is a "hierarchy" when finding the trapping region i X has to be bounded before a bound on = is denoted the 'minimal' trapping region.Notice that any fixed point of Equation ( 1) is contained in U.

All Solutions Get Arbitrarily Close to U in Finite Time and Stay Close to U
For any 0 δ > we can choose 0 >  such that the distance between elements of ( ) . We will prove that for any 0 >  any solution enters ( ) W  in finite time (the time depends on the initial condition).Since ( ) W  is a trapping region this means that the solution stays less than δ away from U for all future time.

Lemma 3. Let
Proof.Follows by the comparison theorem for integrals.
> and since f is continuous there exists τ ′ ′ such that 0 τ τ ′′ ′ < < with ( ) X  is continuous then 1 X  has a maximum 1 0 m < on 1 K .Using lemma 3 and lemma 4 with and ( ) ( ) ( ) there exists a finite time 1 τ such that ( ) ( ) Hence by the proof of theorem 2 1 X stays in this region for all future time.If ( ) I  we will have to repeat the argument.In general assume ( ) ( ) ( ) then we are done.If ( ) ; .
Then by lemma 3 and 4 there exists trapped in this set for all future time.This argument ensures there exists n τ < ∞ such that ( ) ( ) W  is a trapping region, a solution once in ( ) W  for all future time.We emphasize that j X ,

{ }
2, , j n ∈  may be increasing for some time for some initial conditions outside U.
U is the 'minimal' trapping region.However, if ( ) n R X is strictly positive on ( ) then a smaller trapping region can be found using a lower bound on the derivatives which we will not pursue further here. .

Sufficient Criteria for a Globally Stable Fixed Point
This means H is the restriction of R to 0 D .H only attains non negative values since this is the case for R. To continue we assume H is positive and a contraction on 0 D which means we assume there exists ( ) H y H y p y y − ≤ − and ( ) , y y D ∀ ∈ .This ensures the existence of a unique fixed point of Equation (6).Moreover, any solution of Equation ( 6) in 0 D converge to the unique fixed point of the system which will be proven in this section.The approach relies on squeezing the solutions of the Equation ( 6) with solutions of linear systems.The contraction property then ensures the upper and lower bound converge towards the same limit why the solutions of Equation ( 6) must converge to that limit, the unique fixed point of Equation (6).For Thus, two linear systems of differential equations can be constructed with initial condition ( ) ( ) ( ) , for 2, , .
sums appearing in the solutions of the linear systems get arbitrarily small for increasing time.This means for any 1 0 there exists 1 0 ; , for for 1, , .
This means especially ( ) Choosing 1   sufficiently small makes since 0 1 ≤ .The argument may be repeated using the solutions of linear differential equations as bounds for the non-linear system but with a restricted domain for n X .Define ( ) From above there exists a finite time 1 τ such that ( ) where are well defined and finite.Then Since by assumption and ensures .  Due to the squeezing of the solutions using linear systems we have shown that if ( ) We may repeat the argument with bounding the solutions of the non-linear differential equations by solutions to linear systems of differential equations.This means there exists k τ < ∞ such that if ( ) > .We now want to prove that k D converges to { } nss X meaning that all points of k D converge to nss X .The idea of the proof is based on the convergence of ( ) by the Banach Fixed Point Theorem [24].However, there is also a large number of 'error terms' that we have to control.This is done by using the contraction property of H as well as a specific choice of the sequence k   which is decreasing and positive.Thus, all solutions of the non-linear differential equations converge to the unique fixed point of the system.We need the following two lemmas to prove this main result.
Lemma 7. Let p be the contraction constant for H. Then Proof.Follows from the contraction property and the triangle inequality. Similarly it follows.Lemma 8. Let p be the contraction constant for H. Then Lemma 7 and 8 means we can bound the maximum and minimum of H applied on a compact interval by the maximal distance between any two points in the interval and H evaluated at an end point of the interval.
As mentioned a specific choice of k   is chosen as a decaying sequence.Introducing a fixed 0 where  is given by Equation (34).Choose For simplicity we put ( ) To simplify notation further we introduce Thus, For later use we emphasize that { } Proof.The proof is by induction.Since 0 0 l l =  and 0 0 u u =  and 1 0 A =   a basis for the induction is justified.Let [ ] , .
Using the contraction property as shown in lemma 7 and lemma 8 , ; and .
From the definitions of , , , , Thus, we have upper and lower bounds for each of the sets ( ) ) 62)-(67) we get and applying Equation ( . Using Equation ( 52) which completes the proof. Lemma 10.Let H be defined as in Equation (23).If H is a contraction and positive on 0 D then a unique fixed point exists of Equation (6).All solutions in ( ) W  converge to the fixed point.Proof.Fix ˆ0 >  .We will first show that for any ( ) Then the convergence of the remaining i X follows easily.
( ) . From Equation (28) and Equation ( 29) this means that ( ) i h τ and ( ) i g τ converge towards the same limit as This means that all solutions with initial conditions in

( )
W  converge to the unique fixed point of Equation (6).Since all solutions outside ( ) W  in finite time then if H is a contraction and positive on 0 D all solutions converge to the fixed point as τ tends to infinity.This implies that no periodic solution exists.Assuming existence of a periodic solution there must be a positive distance between the fixed point and any periodic solution.Since we have just proved that any solution converge to the fixed point then, after some time we have a contradiction which means, there cannot exist any periodic solutions in the trapping region.

Sufficient Criteria for a Contraction
A sufficient, easily applicable criteria for H being a contraction can be formulated [24].
then all solutions of the non-linear system of differential Equation (6) converge to the unique fixed point.However since for this conclusion to hold.
With the results of Section 3 we now have established the main result of global stability of system (6).
∏  a unique, glo- bally stable fixed point exists of system (6).

Discussion
The general formulation and results in this paper guarantee that the hormone levels in the models [1]- [7] [10] [20] stay in a trapping region where non-negative concentrations are impossible which is a physiological necessity.A repeating pattern is often visible in hormone levels.However, for Equation (1) periodic solutions are impossible outside the "minimal" trapping region.This narrows the domain for interesting initial conditions.The one dimensional function contains a lot of relevant information about the system since it determines the fixed point(s) derivative gives a sufficient criterion for a globally stable fixed point.In [3]- [5] [7], the sufficient criteria for a globally stable fixed point are fulfilled for a subset of the physiologically relevant parameter space.In [1], the focus is on local stability of the fixed point.The investigation of global stability using the criteria found in this paper seems straight forward.Similarly for [18] when the external input to the system is independent of time.
A model of mRNA and Hes1 protein production fits Figure 1 [25] [26].However, a time delay in the feedback has to be included in order to obtain experimentally observed oscillations.A model of testorone dynamics including delay in the feedback is investigated in [27].Including time delays in models corresponding to Figure 1 has proved useful in the search for oscillatory behaviour [28]- [30].One may wonder whether the feedback itself can cause oscillations or if a time delay needs to be included.The contribution of this paper may help in quickly ruling out oscillatory behaviour in the case of no time delay.
Including time delay in the feedbacks, global stability criteria have been formulated for a subset of possible feedback functions in systems resembling 1 [31].This requires that all feedbacks are monotone functions.The approach is different from ours and relies on control theory.

Summary
A general formulation of an n-dimensional system of differential equations with up to n feedbacks from the n'th variable is formulated.The feedbacks may be non-linear but must be represented by bounded functions which are considered to be the case for some biological systems.Some relevant general results are shown.
• Existence and uniqueness of solutions are guaranteed.
• Non-negative initial conditions cause non-negative solutions for all future time.
• A trapping region, ( ) W  , with non-negative elements exists 0 ∀ ≥  . A "minimal" trapping region, ( ) , exists.The existence of a trapping region is a desirable property if e.g. the system is a model of hormone levels.Then moderate hormone levels are guaranteed for future time if the initial conditions are reasonable.
• All solutions of the system enter ( ) W  in finite time for 0 >  .Then any solution gets arbitrarily close to U in finite time.This eliminates the existence of possible limit cycles outside U.
• At least one fixed point exists and all fixed points are contained in U. Using ( ) ( ) = ∏  a sufficient criteria for uniqueness of the fixed point is ( ) If the feedback functions correspond to negative feedbacks or are independent of n X then a unique fixed point exists.

Figure 1 .
Figure 1.Compartment model of the system.

Lemma 5 ..
Consider Equation(6).For any 0 >  any initial condition leads to a solution in ( ) Assume we have an arbitrary non negative initial condition ( ) ( ) monomiums in τ determined from the initial conditions.If j c τ and ( ) ij d τ are constants.With ( ) ( ) ( )

D 1 D
Equation (36) is well defined and compact and The proof is done by induction.0 u and 0 l are given by the expressions is compact and H is continuous then 0 u and 0 l are well defined and finite.This guarantees that is well defined and compact.Since 0 c

ku
 and k l  are well defined since repeated use of a continuous function maps compact sets into compact sets.kl and k u are crucial for the of 1 k D + .We want to make bounds on k l and k u using k l  and k u  since we know the latter converges.In k D "error terms" ( k   ) are introduced at each step in the sequence.The following lemma helps bounding k D by a series in the 'error terms' and a sequence we can then estimate the two separately.Lemma 9.If H is a contraction and positive on 0 D then ; , .
domain of a function can only increase the minimum and decrease the maximum of the function values.Thus, from Equation

∏ 1 1
 a unique, globally stable fixed point exists.Then the fixed point can be found by choosing any for n → ∞ by Banach Fixed Point Theorem.