Parameter Dependent Stability / Instability in a Human Respiratory Control System Model

In this paper a computational procedure is presented to study the development of stable/unstable patterns in a system of three nonlinear, parameter dependent delay differential equations with two transport delays representing a simplified model of human respiration. It is demonstrated using simulations how sequences of changes in internal and external parameter values can lead to complex dynamic behavior due to forced transitions between stable/unstable equilibrium positions determined by particular parameter combinations. Since changes in the transport delays only influence the stability/instability of an equilibrium position a stability chart is constructed in that case by finding the roots of the characteristic equation of the corresponding linear variational system. Illustrative examples are included.


1.
Introduction.Understanding unstable patterns in respiratory control has great importance when studying various forms of periodic breathing (PB)(including Cheyne-Stokes Respiration (CSR)).It is believed that PB is a feedback induced instability which can be influenced in a complicated way by many system and environmental parameters including the transport delays, the loop gains, and the levels of inspired O 2 and CO 2 see e.g.( [4], [5], [6], [7], [8], [9], [11]).A mathematical model of the respiratory control system is based on chemical balance equations (motivated by the work of Khoo et al.([6]), is described by a system of nonlinear, parameter dependent, delay differential equations.Internal system parameters can vary from person to person and for individuals they can change with age, heart condition, and activities such as sleep and exercise .Note that during sleep chemical control of respiration is dominant.There are two negative feedback loops systems in the central control loop with delay τ b and the peripheral control loop with delay τ p .Physiologically τ b is always greater than τ p due to the fact that the central and peripheral chemoreceptors are located in the brain and the carotid body, respectively.Investigating the regulation of the human respiratory process is very important in understanding symptoms including PB, CSR, Obstructive Sleep Apnea(OSA), Central Sleep Apnea (CSA), altitude sickness, sudden infant death syndrome (SIDS) which are all considered unstable patterns of respiration.Following the pioneering paper of Mackey and Glass [10], Cooke and Turi [2] introduced the necessary theoretical tools to study the local stability of the unique equilibrium in a two state model with a single transport delay.In this paper a three state model with two transport delays is considered and due to the increased complexity in the mathematical model the stability analysis is carried out computationally.(It is not clear at this time if the theoretical results in [2] can be extended to the system investigated here.)In particular, the following issues are addressed: 1.A physiologically correct model of respiration, where the transport delay from the lung to the brain compartment (τ b ) is always greater than the transport delay from the lung to the carotid body (τ p ) is considered.Simulations are performed corresponding to the time dependency of some of the internal and external system parameters to investigate the corresponding transition pattern in the system.For fixed parameter values (other than the transport delays, which do not play a role in the determination of the unique equilibrium) a linear variational system is derived and used to obtain the stability chart with respect to the transport delays.Responses to disturbance are computed at stable,unstable, and boundary points.2. Stability of the respiratory system is parameter dependent.When model parameters vary the unique equilibrium and also the corresponding stability chart are forced to change correspondingly.In this study special attention is paid to parameters such as central control gain (G C ), peripheral control gain (G P ), and the transport delays, (τ b ) and (τ p ) in order to understand the role of the central and peripheral control loops in the regulation of the respiratory system.The rest of the paper is organized as follows: In Section 2 the respiratory control system model is described and analyzed.In Section 3 the computational tools used in the paper are discussed.Case studies are presented in Section 4. Finally, Section 5 contains some discussion on our findings.
2. Respiratory Model.We are studying breathing patterns during sleep.Note that sleep is a dynamic activity on its own with several possible sleep stages.Some of the model parameters assume different values corresponding to different sleep stages.An overview of the respiratory control system is shown in Figure 1.It is based on two compartments (lung and brain) for blood gas exchange and storage.The model equations are balance equations for CO 2 and O 2 concentrations in the lung compartment and and for CO 2 concentration in the brain compartment.where the ventilation function, V (t) is given by Note that in the ventilation equation the first term, depending on delayed values of P AO2 and P ACO2 , is the contribution of the peripheral control loop; while the second term, depending on the delayed value of P BCO2 , is representing the contribution of the central control loop.The notation, x + , defined as x ≥ 0 0 otherwise, stands for the positive part of a quantity, x.In (1)-( 4) P X denotes the partial pressure of CO 2 or O 2 in the alveolar air (X = A), in the brain (X = B), in the venous blood (X = V ), and in the inspired air (X = I).Similarly, V X refers to the effective volume of CO 2 and O 2 in the lung (X = A), and brain (X = B) compartments, respectively.Furthermore, V is the total minute ventilation, E F the dead space factor, G C and G P are the central and peripheral controller gains, I C and I P are the apneic thresholds for the central and peripheral carbondioxide response, M R BCO 2 is the metabolic production rate for carbondioxide in the brain compartment, K CO2 and K BCO2 are the dissociation constants for carbondioxide at the artery and the brain, Q and QB are the cardiac output and the rate of blood flow to the brain compartment,τ p and τ b are the blood transport delays from lung compartment to the peripheral and central chemoreceptors, and M V ,M A , B V and B A are constants.To simplify the notation we rewrite equations ( 1) -( 4) as: where x 1 (t) = P ACO 2 (t), x 2 (t) = P AO 2 (t), and x 3 (t) = P BCO 2 (t) are the states, and the constants a ij with i, j = 1.2.3 are given by a 11 = 863 , and The following numerical values will be used for the biological parameters appearing in the a ij 'sthroughout the paper: QB = 0.75, Q = 6.0 in [l/min]; K CO2 = .0057,K BCO 2 = .0065;P VCO 2 = 46.0, It is easy to see that the system described by equations ( 5) -( 8) has a unique positive equilibrium, (x 1 , x2 , x3 ).
In particular, at equilibrium we have  9) and (10) we get x2 in terms of x1 Substituting the expressions for x2 and x3 into (9) we get a nonlinear equation for x1 which has a unique positive solution (constant steady state of system ( 5)-( 8)).
Note that the unique positive equilibrium of system ( 5)-( 8) depends on the parameters G P , G C , P ICO2 , and P IO2 , which means that a change in the values of those parameters forces the system to transit to the corresponding "new" equilibrium.During the transition process unstable respiratory patterns may arise which then can lead to further changes in parameter values (e.g., due to the resulting changes in sleep stages).The transport delays do not influence the location of the equilibrium, but can change its stability.
The local stability of the equilibrium of system ( 5)-( 8) can be investigated using the linear variational system obtained by linearizing system (5)-( 8) around the equilibrium.Let Then the linear variational system for ( 5)-( 8) is given by

PARAMETER DEPENDENT STABILITY/INSTABILITY IN A HUMAN RESPIRATORY 647
where The associated characteristic equation is For the computation of the roots of ( 14) the software TRACE-DDE is used providing stability/instability information on the equilibrium of the nonlinear system as long as the rightmost eigenvalues are not located on the imaginary axis.Note that simulation can always be performed corresponding to some perturbation of the constant initial function at the equilibrium point.
3. Numerical Approximation and Methodology.In this section following [3] we introduce a simple numerical scheme to approximate solutions of equations ( 1) -( 4).Let h be a fixed positive constant, and define the notation [t] h = t h h, where [•] is the greatest integer function.Note that [t] h as a function of t is piecewise constant, since [t] h = nh for t ∈ [nh, (n + 1)h).For a fixed h > 0 we associate the system ẋ1 to ( 1)-( 4).The solution x 1 h , x 2 h and x 3 h of equations ( 15) -( 18) is defined as continuous functions, which are differentiable and satisfy system (15) -(18) on each interval (nh, (n + 1)h) (n = 0, 1, 2, . ..).Since the right-hand-sides of (15)-( 18) are constant on each interval [nh, (n + 1)h), we get that x 1 h , x 2 h and x 3 h are piecewise linear continuous functions (linear spline functions).Therefore, they are determined by their values at the mesh points nh.Introduce the sequences  15) -( 17) from nh to t and taking the limit t → (n + 1)h−, we get by simple calculation that u n , v n and w n satisfy for n = 0, 1, 2, . ..,Where Note that for negative integer n the sequences u n and v n are defined by u n = x 1 (nh) and v n = x 2 (nh), the initial functions corresponding to the original system (1)-( 3).Therefore the sequences u n , v n and w n are well-defined and can be easily generated and the solutions of equations ( 15)-( 18) are uniquely determined.It follows from the result of [3] that x i h converge uniformly to x i (t) for i = 1, 2, 3 on any compact time interval.We shall use ( 19) -( 22) with changing parameter values and former equilibrium values as constant initial functions to run simulations on the respiratory system.
For the computation of the characteristic roots of ( 14) corresponding to given τ p and τ b at an equilibrium point identified by the values of G P , G C , P IO 2 and P ICO 2 we use the Matlab package TRACE-DDE (see [1]).We repeat the above calculation for various delay values and construct an approximate stability chart by identifying the behavior of the rightmost characteristic root of ( 14).If this root is close to the imaginary axis, we also run simulations on the nonlinear model in order to improve the accuracy of the approximate stability chart.and calculate the corresponding steady-state, X = (39.080,76.207, 47.695).Using the linear variational system for different τ p and τ b values we obtain an approximate stability chart by repeatedly computing the roots of the characteristic equation (see Figure 2.For illustration we select points inside (A), on the boundary (B), and outside (C) of the stability region and display the roots of (14) on Figures 3, 4, and 5, respectively.Starting from the constant initial function φ(t) = X and switching off the ventilation for 10 seconds at t= 0 (i.e., V (t) = 0 for t ∈ [0, 10]) we can observe numerically the responses of the system to this perturbation (see Figures. 6,7 and 8).115.17,0.22) to (109.935, 0.21) and (104.70,0.20), respectively.As the consequences of these changes the equilibrium positions of system (1) -(4) will move from X = (39.080,76.207, 47.695) to X = (38.885,72.987, 47.500) and X = (38.677,69.734, 47.292).We compute the approximate stability charts corresponding to those equilibrium points (see Fig 9).Example 3 : We fix P I O2 = 115.17, P I CO2 = 0.22 and change G P and G C from (35,1.2) to (45,1.2) and then to (45,1.5).The equilibrium position of system (1) -(4) will move from X = (39.296,75.029, 47.911) to X = (39.080,76.207, 47.695) and to X = (38.899,77.149, 47.514).We compute the approximate stability charts corresponding to those equilibrium points (see Figure 10).

Example 4
We fix P I O2 = 115.17,P I CO2 = 0.22, τ p = 9, τ b = 15 and vary (G P ) and (G C ). ( Of course, the equilibrium position will be changing).We find points with coordinates, (G P , G C ), for which the corresponding equilibrium is stable/unstable (see Figure 11).The arrows on Figure 11 indicate a sequence of G P and G C changes which moves the system from stable (A) to unstable (B) to stable (C) to unstable(D),... etc. Example 5 We fix P I O2 = 115.17,P I CO2 = 0.22, G P = 35 and G P = 1.5 and calculate X = (39.089,76.156, 47.705).We change G P and G C at t=0 from (35, 1.5) to (35, 1.2).Then at t = 20 and t = 40 we change them again from (35, 1.2) to (45, 1.2) and from (45, 1.2) to (45, 1.5), respectively and compute the corresponding trajectories (see Figure 12). 5. Discussion.We considered both peripheral and central controls in a three state model of the human respiratory control system.As new result we observed that the presence of the two control loops with different transport delays brought about the possibility of multiple stability switching as τ p and τ b varied.For example, for the transport delay from lung to carotid body (τ p ), Figure 2 indicated the possibility of moving from unstable to stable condition by either increasing or decreasing the lung to brain transport delay, (τ b ).Changing heart conditions lead to changes in both of the parameters τ p and τ b and the combination of those changes can lead to changes in stability/instability of the equilibrium in the respiratory system.Note that one can identify these parameters and their time dependency from measured data (see e.g., [4]) for details) and then use the method presented here to assess stability/instability.We demonstrated on Figures 9 and 10 how changes of individual parameter values influence the stability of the corresponding equilibrium point.Then on Figure 11 we investigated the simultaneous variations of G p and G C .Finally, on Figure 12 = a 31 + a 32 ( x1 − x3 ), (11) and V = G P exp(−.05x 2 )(x 1 − I P ) + + G C (x 3 − a 33 − I C ) + (12) where V is the equilibrium value of the ventilation.From equation (11) we obtain x3 in terms of x1 .Writing V = a 11 − a 12 x2 a 13 ( x1 − P ICO 2 ) = a 21 + a 22 x2 a 23 (P IO2 − x2 ) using ( ), and w n = x 3 h (nh) and let k = [ τp h ] and l = [ τ b h ].Then integrating equations (