\documentclass[reqno]{amsart} \usepackage{graphicx} \AtBeginDocument{{\noindent\small 2004 Conference on Diff. Eqns. and Appl. in Math. Biology, Nanaimo, BC, Canada.\newline {\em Electronic Journal of Differential Equations}, Conference 12, 2005, pp. 65-78.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu \newline ftp ejde.math.txstate.edu (login: ftp)} \thanks{\copyright 2005 Texas State University - San Marcos.} \vspace{9mm}} \setcounter{page}{65} \begin{document} \title[\hfilneg EJDE/Conf/12 \hfil Numerical stability] {Numerical stability analysis in respiratory control system models} \author[L. E. Koll\'ar, J. Turi \hfil EJDE/Conf/12 \hfilneg] {L\'aszl\'o E. Koll\'ar, J\'anos Turi} % in alphabetical order \address{L\'aszl\'o E. Koll\'ar \hfill\break Department of Applied Sciences, The University of Quebec at Chicoutimi, 555 Boul.~de l'Universite, Chicoutimi, Quebec G7H 2B1, Canada} \email{laszlo\_kollar@uqac.ca} \address{J\'anos Turi \hfill\break Department of Mathematical Sciences, The University of Texas at Dallas, P. O. Box 830688, MS EC 35, Richardson, Texas 75083-0688, USA} \email{turi@utdallas.edu} \date{} \thanks{Published April 20, 2005.} \subjclass[2000]{92C30, 93C23} \keywords{Delay differential equations; human respiratory system; \hfill\break\indent transport delay; numerical analysis} \begin{abstract} Stability of the unique equilibrium in two mathematical models (based on chemical balance dynamics) of human respiration is examined using numerical methods. Due to the transport delays in the respiratory control system these models are governed by delay differential equations. First, a simplified two-state model with one delay is considered, then a five-state model with four delays (where the application of numerical methods is essential) is investigated. In particular, software is developed to perform linearized stability analysis and simulations of the model equations. Furthermore, the Matlab package \textsc{dde-biftool} v.~2.00 is employed to carry out numerical bifurcation analysis. Our main goal is to study the effects of transport delays on the stability of the model equations. Critical values of the transport delays (i.e., where Hopf bifurcations occur) are determined, and stable periodic solutions are found as the delays pass their critical values. The numerical findings are in good agreement with analytic results obtained earlier for the two-state model. \end{abstract} \maketitle \numberwithin{equation}{section} \section{Introduction} In the present work we examine stability/instability of breathing patterns in two models of the human respiratory system described by nonlinear parameter dependent delay differential equations with discrete circulatory transport delays. Due to the complexity of the delay systems involved, the only feasible way to carry out such study is by computational means. The main goal of this work is to develop numerical tools in order to study the stability of the human respiratory system with the additional benefit that these tools may also be applied to examine other delay differential equations of the same type. In Section 2 these numerical tools are described briefly. We developed a Matlab code based on \cite{Manitius} to compute the roots of characteristic functions associated with linear differential-difference equations. This code is also applicable to find the critical value of the delay where stability changes. An approximation technique for the simulation of nonlinear delay equations is also described which is widely applicable to simulate systems even with time- and state-dependent delays \cite{Hartung}. Moreover, the bases of a Matlab toolbox for numerical bifurcation analysis are sketched \cite{Engelborghs}. This toolbox is appropriate to compute, continue and analyze stability of steady state and periodic solutions of nonlinear systems even with state-dependent delays, to compute and continue steady state fold and Hopf bifurcations and to switch, from the latter, to an emanating branch of periodic solutions. We use these tools to study nonlinear retarded differential-difference equations with constant delays, but some of them are applicable even for more general classes of delay equations, therefore some of the results also cover more general cases. The purpose of the human respiratory system is to exchange the unwanted gas byproducts of metabolism, such as CO$_2$ for O$_2$ which is necessary for metabolism. The CO$_2$ is exchanged for O$_2$ by passive diffusion. The primary determiners of diffusion are the partial pressure gradients across the blood/gas barrier between capillaries and alveoli, where gas transfer occurs. The respiratory control system varies the ventilation rate in response to the levels of \( \mathrm{CO}_2 \) and \( \mathrm{O}_2 \) in the body. Delay is introduced into the control system due to the physical distance which CO$_2$ and O$_2$ levels must be transported to the sensory sites before the ventilatory response can be adjusted. Detailed discussion of the human respiratory system and its control mechanism can be found in \cite{Khoo} and \cite{West}. Models of the respiratory control system date back to the early 1900's to study a wide range of features of this complex system. The phenomena collectively referred to as periodic breathing have important medical implications. Physiological studies led to the hypothesis that periodic breathing is the result of delays in the feedback signals to the respiratory control system. A number of dynamic models of \( \mathrm{CO}_2 \) regulation have been proposed since the 1950's and transport delays have been introduced. A five-state model involving three compartments and two control loops with multiple delays is studied in \cite{Khoo} and a two-state system with one delay modeling partial pressures of CO$_2$ and O$_2$ in the lung and the peripheral controller is considered in \cite{CookeTuri}. In \cite{CookeTuri}, system parameters are kept constant except the delay and analysis is done to illuminate the effect of the delay on the stability of the unique steady state. In Section 3 we study the two-state model presented in \cite{CookeTuri} and the five-state model of \cite{Khoo} numerically by the tools described in Section 2. We linearize the systems, then find the characteristic roots with largest real parts to determine stability/instability of the equilibria. Critical values of transport delays (where Hopf bifurcations occur) can also be obtained numerically by using our Matlab code. In case of the two-state system, the critical delay can also be determined analytically according to the computation proposed in \cite{CookeTuri}. Simulations predict that after the Hopf bifurcation point the equilibrium is unstable but a stable periodic solution appears. Numerical bifurcation analysis gives the branch of the periodic solution including the amplitude and the frequency of oscillations which show good agreement with the simulation results. The obtained periodic solution represents a medically important phenomenon referred to as periodic breathing. Simulations of the five-state system are carried out for three different cases: (i) all the parameters are fixed and delays take the same value, (ii) some of the system parameters are state dependent and delays take the same value, and (iii) some of the system parameters are state dependent and delays may take different values. The results obtained make it possible to find approximately the critical values of the delays. As a particular consequence, our tools are appropriate to determine existence and stability of steady state and periodic solutions of the respiratory system and to find the critical value of the transport delay which has medical importance. \section{Numerical Stability Analysis of Delay Differential Equations} In this section, we describe three numerical tools used in our computations to study the human respiratory system that is modelled by nonlinear delay differential equations of retarded type. The stability of steady state solutions of such systems can be determined by roots of the characteristic equation. The characteristic equations of such systems have infinitely many characteristic roots, but only finitely many of them have real parts greater than any fixed real number. It follows that the stability can be determined by finding finitely many roots numerically. An algorithm to find characteristic roots is exhibited in Subsection \ref{rootfind}. The linearized model provides information about the stability of the steady state solution and it is applicable to determine a critical value of the parameter of a parameter dependent system where the equilibrium loses its stability. However, other methods are needed to study the nonlinear system after the loss of stability. An approximation scheme is presented in Subsection \ref{appr} which is the basis of simulations. Simulations of our models predict that there is a Hopf bifurcation at the parameter value where the equilibrium loses its stability and stable periodic solutions exist if the value of parameter increases. Subsection \ref{bif} gives a brief review of numerical bifurcation analysis techniques which provide suitable tools to find bifurcation points, branches of periodic solutions, and determine stability, amplitude and frequency of periodic solutions. \subsection{Computation of Eigenvalues Associated with Linear Delay Differential Equations}\label{rootfind} This section describes a method to compute the eigenvalues associated with systems of linear retarded differential-difference equations. The method proposed in \cite{Manitius} finds the roots of a given characteristic equation. The eigenvalues contained in a bounded region around the origin are approximately computed by a combinatorial algorithm suggested in \cite{Kuhn} and referred to as Kuhn's method. The eigenvalues of large modulus are computed using asymptotic formulae obtained directly from the coefficients of the characteristic equation \cite{Bellman}. The well known Newton's method is used for refinement of the approximations, and to verify that all the eigenvalues have been found in a finite domain of the complex plane, a procedure proposed by \cite{Carpentier} is applied. We have developed a Matlab implementation of the algorithm based on the above discussion \cite{thesis}. The user enters the characteristic equation together with its derivative (because Newton's method requires that), a rectangle region and specifies the mesh size on it for Kuhn's method, the matrix \( P \) for computing the roots of large modulus, and the radius of the disc where the test on the number of eigenvalues is desired. The program computes roots and also plots them. Our code is applicable to find characteristic roots of linear systems as the delay changes and to determine its critical value where the equilibrium loses its stability. We need further tools to examine what happens in the nonlinear system after the loss of stability. Two suitable numerical methods are discussed in the following two subsections. \subsection{Approximation Schemes for Delay Differential Equations}\label{appr} One way to examine nonlinear delay differential equations numerically is to discretize the initial function and the equations with respect to time and approximate the solution by solving the discretized system. A numerical approximation technique is presented in \cite{Hartung}, while the convergence of this method is proved in \cite{Hartungconv}. The computational scheme is based on approximation of delay differential equations by equations with piecewise constant arguments. The approximating equations, generated in the above process, lead naturally to discrete difference equations, well suited for computational purposes, and thus provide an approximation framework for simulation studies. This method is applicable for any delay differential equation with several discrete delays even if delays are time- and state-dependent. \subsection{Numerical Methods for Bifurcation Analysis}\label{bif} An effective way to study nonlinear delay differential equations is bifurcation analysis. We use \textsc{dde-biftool} v.~2.00 which is a collection of Matlab routines for numerical bifurcation analysis of systems of delay differential equations with several constant and state-dependent delays \cite{Engelborghs}. Models of the respiratory system that we consider contain one constant delay which is chosen to be the bifurcation parameter. This Matlab package finds the equilibrium of our system, and determines its stability by linearization and finding characteristic roots with greater real parts than a prescribed value. Changing the bifurcation parameter, a critical value is obtained where a complex conjugate characteristic root pair crosses the imaginary axis and the equilibrium loses its stability. According to the Hopf bifurcation theorem \cite{Stech85} and \cite{Stech87}, there is a Hopf bifurcation at this value of the bifurcation parameter in the nonlinear system. After the Hopf bifurcation point, the steady state solution is unstable but a stable periodic solution appears. \textsc{dde-biftool} finds that periodic solution in the nonlinear system by continuation. The amplitude and the frequency of the periodic motion can be computed for every value of the bifurcation parameter. At last, the stability of the periodic solution is determined by finding the Floquet multipliers. In the \textsc{dde-biftool}, the user must provide the governing equations, parameters (including delays), and system derivatives with respect to state variables and parameters up to second order. Then using an execution file, the steady state solutions with their stability, the Hopf bifurcation point and periodic solutions with their stability, as well as amplitude and frequency can be obtained and the results such as location of characteristic roots, bifurcation diagrams, a period of periodic solutions and location of Floquet multipliers can be plotted. \section{Mathematical Models of the Respiratory System} We study a simplified two-state model and a more complete, five-state model of the respiratory system. The five-state model is constructed in \cite{Khoo} where lung compartment, brain compartment and general tissue compartment are described. The equations for the model arise from straightforward development of mass balance equations utilizing Fick's law, Boyle's law and variations of Henry's law relating the concentration of the gas in the solution to the partial pressure of the gas interfacing with the solution. The symbol sets used in the model equations are provided in Appendix 5.1. In \cite{Khoo} the state equations for the system are described in terms of concentrations. Here, we employ partial pressures instead using the relations: \begin{equation} \label{dissoc} \begin{aligned} C_{a_{CO_2}} &= K_{CO_2} P_{a_{CO_2}} + K_1, \\ C_{V_{CO_2}} &= K_{CO_2} P_{V_{CO_2}} + K_1, \\ C_{B_{CO_2}} &= K_{B_{CO_2}} P_{B_{CO_2}} + K_1, \\ C_{a_{O_2}} &= m_a P_{a_{O_2}} + B_a, \\ C_{V_{O_2}} &= m_v P_{V_{O_2}} + B_v. \end{aligned} \end{equation} The dissociation curves for O$_2$ are represented approximately by piecewise linear functions. The dependence of parameters \( m_a, B_a \) and \( m_v, B_v \) on \( P_{a_{O_2}} \) and \( P_{V_{O_2}} \), respectively, are provided in Table 2 in Appendix \ref{parameters}. Straightforward manipulations of the state equations yield: \begin{equation}\label{model5eq} \begin{aligned} \frac{\mathrm{d}P_{a_{CO_2}}(t)}{\mathrm{d}t} &= \frac{863 Q K_{CO_2} ( P_{V_{CO_2}}(t-\tau_V)-P_{a_{CO_2}}(t) )}{M_{L_{CO_2}}} \\ &\quad+\frac{E_F V_I (t) ( P_{I_{CO_2}}-P_{a_{CO_2}}(t) )}{M_{L_{CO_2}}}\,, \\ \frac{\mathrm{d}P_{a_{O_2}}(t)}{\mathrm{d}t} &= \frac{863 Q ( m_v P_{V_{O_2}}(t-\tau_V)-m_a P_{a_{O_2}}(t)+B_v-B_a )} {M_{L_{O_2}}} \\ &\quad+\frac{E_F V_I (t) ( P_{I_{O_2}}-P_{a_{O_2}}(t) )}{M_{L_{O_2}}}\,, \\ \frac{\mathrm{d}P_{B_{CO_2}}(t)}{\mathrm{d}t} &= \frac{MR_{B_{CO_2}}}{M_{B_{CO_2}}K_{B_{CO_2}}} + \frac{Q_B (P_{a_{CO_2}}(t-\tau_B)-P_{B_{CO_2}}(t))}{M_{B_{CO_2}}}\,, \\ \frac{\mathrm{d}P_{V_{CO_2}}(t)}{\mathrm{d}t} &= \frac{MR_{T_{CO_2}}}{M_{T_{CO_2}}K_{CO_2}} + \frac{Q_T (P_{a_{CO_2}}(t-\tau_T)-P_{V_{CO_2}}(t))}{M_{T_{CO_2}}}\,, \\ \frac{\mathrm{d}P_{V_{O_2}}(t)}{\mathrm{d}t} &= \frac{Q_T (m_a P_{a_{O_2}}(t-\tau_T)-m_v P_{V_{O_2}}(t)+B_a-B_v) -MR_{T_{O_2}}}{M_{T_{O_2}}m_v}\,, \end{aligned} \end{equation} where \( V_I \) is the ventilation rate or ventilation function which depends on the signals sent from the peripheral sensors and includes transport delay to the peripheral controller. The ventilation function is assumed in the form \begin{equation}\label{model5vf} \begin{aligned} V_I(t)&= G_P \mathrm{e}^{-0.05 P_{a_{O_2}}(t-\tau_a)} \big( P_{a_{CO_2}} (t-\tau_a)-I_P \big) \\ &\quad+ G_C \big( P_{B_{CO_2}}(t) -\frac{MR_{B_{CO_2}}}{K_{CO_2}Q_B}-I_C \,. \end{aligned} \end{equation} The two-state model considers lung compartment only, and contains two equations describing CO$_2$ and O$_2$ arterial partial pressures and a control equation responsive to these arterial partial pressures with one transport delay to the peripheral controller. The state equations are the first two equations of system (\ref{model5eq}) without delay, while the ventilation function \( V_I \) is assumed in the form \begin{equation}\label{model2vf} V_I(t)=G_P \mathrm{e}^{-0.05 P_{a_{O_2}}(t-\tau)} ( P_{a_{CO_2}} (t-\tau)-I_P )\,, \end{equation} where notations are as before, and \( \tau \) is the only delay. The following simplified form is derived in \cite{CookeTuri} \begin{equation} \begin{aligned} \frac{\mathrm{d} \tilde x(t)}{\mathrm{d}t}&= p-\alpha W(t) ( \tilde x(t)-x_I ), \\ \frac{\mathrm{d} \tilde y(t)}{\mathrm{d}t}&= -\sigma+\beta W(t) ( y_I-\tilde y(t) ), \end{aligned}\label{model1eq} \end{equation} where the notation is provided in Appendix \ref{symbols}. System (\ref{model1eq}) can be transformed to a more convenient form by introducing the new state variables \( x(t)=a ( \tilde x(t)-x_I ) \) and \( y(t)=b ( y_I-\tilde y(t) ) \). Setting \( a=1/p \) and \( b=1/\sigma \), we obtain the equations \begin{equation} \begin{aligned} \frac{\mathrm{d}x}{\mathrm{d}t}&= 1-\alpha V(t) x(t), \\ \frac{\mathrm{d}y}{\mathrm{d}t}&= 1-\beta V(t) y(t), \end{aligned} \label{mod1} \end{equation} with the ventilation function \( V(t)=V ( x(t-\tau), y(t-\tau) )=W ( \tilde x (t-\tau), \tilde y(t-\tau) ) \). Note that state variables are concentrations in system (\ref{mod1}) and partial pressures in (\ref{model5eq}). However, the relationship between concentration and partial pressure is assumed to be linear, therefore system (\ref{mod1}) can easily be transformed into a system that contains partial pressures. \subsection{Numerical Results on the Two-State Model}\label{model1} In this subsection, we consider system (\ref{mod1}). First, we summarize theoretical results, then we examine this system for a certain parameter setting. Parameters \( \alpha \) and \( \beta \) are fixed, the ventilation function is given and we let \( \tau \) be a free parameter. The equilibrium is determined which is the same for all values of the delay and we study its stability as the delay changes. The linear stability analysis is carried out by using the rootfinder discussed in Subsection \ref{rootfind}. As a critical value of the delay is found where the equilibrium loses its stability, following examination of the nonlinear system is needed. We present simulation results and the results of numerical bifurcation analysis. The following assumptions are considered in our model. The arterial CO$_2$ concentration is always larger than its inspired value and the arterial O$_2$ concentration never exceeds its inspired value. It means that \( \tilde x>x_I \) and \( \tilde y0 \) and \( y>0 \). It also appears to be biologically realistic to assume that the ventilation function has the following properties: (i) \( V(x,y) \) is increasing in both \( x \) and \( y \), (ii) \( V(x,y) \geq 0 \) and \( V(0,0)=0 \), (iii) \( V(x,y) \) is differentiable, and (iv) \( V_x=\partial V(x,y)/\partial x>0, \quad V_y=\partial V(x,y)/\partial y>0 \). It is proved in \cite{CookeTuri} that if the above assumptions are held then system (\ref{mod1}) has a unique positive equilibrium. First, we recall the conditions for asymptotic stability of the equilibrium. Let \( (\bar x,\bar y) \) denote the equilibrium then letting \( \xi(t)=x(t)-\bar x, \eta(t)=y(t)-\bar y \) in system (\ref{mod1}) and removing the nonlinear terms, we obtain the linear variational system \begin{equation}\label{mod2lin} \frac{\mathrm{d}}{\mathrm{d}t} \begin{pmatrix} \xi(t) \\ \eta(t) \end{pmatrix} +A \begin{pmatrix} \xi(t) \\ \eta(t) \end{pmatrix} +B \begin{pmatrix} \xi(t-\tau) \\ \eta(t-\tau) \end{pmatrix}= \begin{pmatrix} 0 \\ 0 \end{pmatrix}\,, \end{equation} where $$ A=\begin{pmatrix} \alpha \bar V & 0 \\ 0 & \beta \bar V \end{pmatrix} , \quad B= \begin{pmatrix} \alpha \bar x \bar V_x & \alpha \bar x \bar V_y \\ \beta \bar y \bar V_x & \beta \bar y \bar V_y \end{pmatrix} , $$ \( \bar V=V(\bar x, \bar y), \bar V_x=V_x (\bar x, \bar y) \) and \( \bar V_y=V_y (\bar x, \bar y) \). The following statements are proved in \cite{CookeTuri} \begin{enumerate} \item If \( \bar V \geq \bar x \bar V_x+\bar y \bar V_y \), then the equilibrium \( (\bar x,\bar y) \) is asymptotically stable for all delay \( \tau \geq 0 \). \item If \( \bar V < \bar x \bar V_x+\bar y \bar V_y \), then there exists \( \tau_0 > 0 \) such that the equilibrium \( (\bar x,\bar y) \) is asymptotically stable if \( 0 \leq \tau < \tau_0 \) and unstable if \( \tau > \tau_0 \). \end{enumerate} Obviously, the latter case is interesting and subject of further study. If \( \tau \) is regarded as a parameter, then as \( \tau \) passes through its critical value \( \tau_0 \) the equilibrium \( (\bar x, \bar y) \) loses its stability. Then it can be shown that there is a Hopf bifurcation with emergence of a periodic solution. The critical delay \( \tau_0 \) is also computed in \cite{CookeTuri}. In the following, we examine system (\ref{mod1}) by setting \( \alpha= 0.5, \beta=0.8 \) and with the ventilation function \( V(t)=0.14 \mathrm{e}^{-0.05 ( 100-y(t-\tau) )} x(t-\tau) \). The state variables of this model are concentrations, while the five-state model considers partial pressures, therefore the form of the ventilation function slightly differs from Equation (\ref{model2vf}), and so do the values of the parameters from those given in Table 1 in Appendix \ref{parameters} for the five-state model. The equilibrium \( (\bar x, \bar y) \) can be found by using the Matlab code based on the discussion of Subsection \ref{appr}. We can easily check that \( V < x V_x + y V_y \) for any \( x > 0 \) and \( y > 0 \), so there exists a critical delay \( \tau_0 \) where the equilibrium loses its stability and it can be computed by following the above discussion. For the given \( \alpha \) and \( \beta \), the equilibrium is \( (\bar x, \bar y)=(29.1842, 18.2401) \) and the critical delay is \( \tau_0=30.8017 \) [s] (if time is given in seconds). The characteristic root with largest real part can be obtained for any value of the parameter, i.e., the transport delay, by using the rootfinder Matlab code discussed in Subsection \ref{rootfind}. The real part of the characteristic root with largest real part is computed for several values of \( \tau \) approaching \( \tau_0 \) in each step. Figure \ref{mod1tau0}(a) shows the computed values and the approximate \( \tau_0 \) is the arithmetic mean of the two delays when the computed real parts are closest to zero. If the desired accuracy is 0.01 then the approximate critical delay is \( \tau_{cr}=30.8076 \) [s]. This coincides with the analytically computed value, since \( |\tau_0-\tau_{cr}|<0.01 \). In Figure \ref{mod1tau0}(b), the characteristic roots of smallest modulus are given for \( \tau= 30.8 \) [s]. It can be seen that a complex conjugate characteristic root pair passes the imaginary axis which predicts that there is a Hopf bifurcation at the critical value of the transport delay. It is concluded that the equilibrium is stable if \( \tau<\tau_0 \) and unstable if \( \tau>\tau_0 \). The naturally arising question is what happens in the nonlinear system if we pass the critical delay. \begin{figure}[t] \begin{center} \includegraphics[width=0.45\textwidth]{fig41a} \quad \includegraphics[width=0.45\textwidth]{fig41b} \end{center} \caption{(a) The real parts of the eigenvalue with largest real part \enskip (b) Eigenvalues of smallest modulus for $\tau=30.8$ [s]} \label{mod1tau0} \end{figure} Simulation results are drawn in Figure \ref{mod1sim}. The CO$_2$ and O$_2$ concentration, i.e., \( x(t) \) and \( y(t) \) are given in time for \( \tau=15 \) [s] in Figures \ref{mod1sim}(a) and \ref{mod1sim}(c), respectively. The equilibrium is stable, \( x \) and \( y \) approach \( \bar x=29.1842 \) and \( \bar y=18.2401 \), respectively. The same relationships are shown in Figures \ref{mod1sim}(b) and \ref{mod1sim}(d) for \( \tau=40 \) [s]. According to the stability analysis, the equilibrium is unstable and these figures also show that \( x \) and \( y \) do not converge to \( \bar x \) and \( \bar y \), respectively. The solutions appear periodic. \begin{figure}[t] \begin{center} \includegraphics[width=0.45\textwidth]{fig42a}\quad \includegraphics[width=0.45\textwidth]{fig42b}\\[5pt] \includegraphics[width=0.45\textwidth]{fig42c}\quad \includegraphics[width=0.45\textwidth]{fig42d} \end{center} \caption{Simulation results (a) \( x(t) \) for \( \tau=15 \) [s] \enskip (b) \( x(t) \) for \( \tau=40 \) [s] \enskip (c) \( y(t) \) for \( \tau=15 \) [s] \enskip (d) \( y(t) \) for \( \tau=40 \) [s]} \label{mod1sim} \end{figure} The Matlab package \textsc{dde-biftool} mentioned in Subsection \ref{bif} can also be used to find the equilibrium by entering an initial approximation. (Recall that the equilibrium is stable if \( \tau<\tau_0 \) and unstable if \( \tau>\tau_0 \).) The \textsc{dde-biftool} located a Hopf bifurcation point at \( \tau_0 \). This is a supercritical Hopf bifurcation where a branch of periodic solutions emerges. This branch is shown in Figure \ref{mod1bif}(a) where the bifurcation parameter is the delay \( \tau \) and the amplitude of \( x \) can be seen as the delay changes. The Floquet multipliers can be determined for any value of the delay. Since all the multipliers are inside the unit circle of the complex plane, except the trivial one which is exactly one, periodic solutions are stable. The amplitude and also the frequency of periodic solutions can be computed for any value of the delay. As an example, if \( \tau=40 \) [s] then the amplitude of \( x \) is 39.8 as it can also be checked approximately in Figure \ref{mod1bif}(a), while the period is 111 [s], so the frequency is 0.009 [1/s]. Compare these results with those obtained by the simulations, i.e., see Figure \ref{mod1sim}(b). The maximum of the periodic solution is 47.4, while the minimum is 7.7, so the amplitude is 39.7. There are 8 periods between the maxima at \( t=97 \) [s] and \( t=985 \) [s], so the period is (985-97)/8=111 [s]. We can conclude that the simulation results and the results of numerical bifurcation analysis coincide very well. An interesting periodic solution is shown in Figure \ref{mod1bif}(b). The transport delay is \( \tau=65 \) [s] and the amplitude of \( x \) is 69.3. Since \( \bar x=29.1842 \) in the equilibrium, that would imply a negative value for the minimum of \( x \). However, in reality, the CO$_2$ concentration is never less than the inspired value, i.e., \( x>0 \) (and similarly, \( y>0 \)). Our model considers these conditions and therefore the peaks of periodic solutions do not go under zero. \begin{figure}[htbp] \begin{center} \includegraphics[width=0.45\textwidth]{fig43a} \quad \includegraphics[width=0.45\textwidth]{fig43b} \end{center} \caption{(a) Amplitude of \( x \) as \( \tau \) changes \enskip (b) 1 period of \( x \) and \( y \) for \( \tau=65 \) [s]} \label{mod1bif} \end{figure} \subsection{Numerical Results on the Five-State Model} The five-state model (\ref{model5eq}) with the ventilation function (\ref{model5vf}) is examined numerically in this section. Due to the complicated form of this system, theoretical study similar to that of the preceding subsection is not carried out here. This emphasizes the importance of the numerical tools that may be used to investigate delay differential equations. We approach the final and most complete model of the present study in three steps. First, all the parameters of the system are fixed and it is assumed that all four delays take the same value. The equilibrium of this system is computed by the Matlab code based on the approximation scheme presented in Subsection \ref{appr}. A critical value of the delay where the equilibrium becomes unstable and stable periodic solutions occur is determined by applying the rootfinder Matlab code described in Subsection \ref{rootfind}. Simulations are also carried out for both cases, i.e., when the equilibrium is stable and when the periodic solutions are stable. In the second step, we assume that some parameters are not constant, but they take different values in different regions of some of the state variables. The critical delay is determined in the region that describes the system in the vicinity of the equilibrium and the system is simulated. Finally, the different values of time delays are also taken into account, and the obtained most complete system is simulated. The parameter values used in the first step are given in Table 1 in Appendix \ref{parameters}. The equilibrium of the system is \( (\bar P_{a_{CO_2}}, \bar P_{a_{O_2}}, \bar P_{B_{CO_2}}, \bar P_{V_{CO_2}}, \bar P_{V_{O_2}}) \) = (38.40, 94.25, 47.63, 47.44, 38.97), and the critical value of the time delay is \( \tau_{cr} = 0.9535 \) [min]. It is calculated by the rootfinder Matlab code in the same fashion as it is described in Subsection \ref{model1}. The characteristic roots of smallest modulus for \( \tau = 1.1 \) [min] are given in Figure \ref{mod5fig}(a). It can be seen that the pair of characteristic roots with the largest real part is already greater than 0, therefore the equilibrium is unstable and the periodic solution is stable. Simulation results are shown in Figures \ref{mod5fig}(b) and \ref{mod5fig}(c) for \( \tau = 0.5 \) [min] and for \( \tau = 1.1 \) [min], respectively. The system was simulated for further values of the delay, and we obtained stable equilibrium for \( \tau = 0.84 \) [min], but periodic solution appeared for \( \tau = 0.85 \) [min]. Hence, the critical value of the delay can be found between these values which means that there is an approximately 10 \% difference between the results obtained by the rootfinder as well as by simulation. \begin{figure}[htbp] \begin{center} \includegraphics[width=0.45\textwidth]{fig44a}\quad \includegraphics[width=0.45\textwidth]{fig44b}\\[5pt] \includegraphics[width=0.45\textwidth]{fig44c}\quad \includegraphics[width=0.45\textwidth]{fig44d}\\[5pt] \includegraphics[width=0.45\textwidth]{fig44e}\quad \includegraphics[width=0.45\textwidth]{fig44f} \end{center} \caption{(a) characteristic roots for \( \tau = 1.1 \) [min], (b) simulation for \( \tau = 0.5 \) [min], (c) simulation for \( \tau = 1.1 \) [min], (d) simulation with state-dependent parameters for \( \tau = 0.5 \) [min], (e) for \( \tau_B = 0.7098 \) [min], \( \tau_T = 0.9102 \) [min], \( \tau_V = 0.9102 \) [min] and \( \tau_a = 0.5100 \) [min], \enskip (f) for \( \tau_B = 1.1830 \) [min], \( \tau_T = 1.5170 \) [min], \( \tau_V = 1.5170 \) [min] and \( \tau_a = 0.8500 \) [min]} \label{mod5fig} \end{figure} From now on four parameters, namely \( m_v, B_v, m_a \) and \( B_a \), are considered as state-dependent parameters. More precisely, \( m_a \) and \( B_a \) depend on \( P_{a_{O_2}} \), while \( m_v \) and \( B_v \) depend on \( P_{V_{O_2}} \). The dependence is shown in Table 2 in Appendix \ref{parameters}. In the equilibrium \( P_{a_{O_2}} = 94.25 \) and \( P_{V_{O_2}} = 38.97 \), thus, the conditions \( P_{a_{O_2}} \ge 70 \) and \( P_{V_{O_2}} < 55 \) are satisfied. It gives exactly the system with the parameters given in Table 1 that was discussed in the previous paragraph. The critical time delay was \( \tau_{cr} = 0.9535 \) [min] according to the rootfinder, and \( \tau_{cr} = 0.85 \) [min] according to simulations. If the different regions defined by the state-dependent parameters are considered, then the critical value is the same with accuracy of 0.01 [min] according to simulations. A simulation result is drawn in Figure \ref{mod5fig}(d) where \( \tau = 0.5 \) [min]. In the last step, we consider that the delays take different values. Measured values of the delays are published in \cite{Khoo} that are obtained after carrying out experiments on healthy people: lung to brain delay \( \tau_B = 0.1183 \) [min], lung to tissue transport delay \( \tau_T = 0.1517 \) [min], venous side transport delay from tissue to lung \( \tau_V = 0.1517 \) [min], and lung to carotid artery delay \( \tau_a = 0.0850 \) [min]. In what follows, the delays are increased, in particular they are multiplied by 6 and 10, and simulation results are presented in Figures \ref{mod5fig}(e) and \ref{mod5fig}(f), respectively. It can be observed that the equilibrium loses its stability when values of delays are between the two sets corresponding to these figures. It means that critical values of the delays are in the same range as in the case of the system with one delay, the critical value of the lung to carotid artery delay, \( \tau_a \), the lung to brain delay, \( \tau_B \), the lung to tissue delay, \( \tau_T \), and the venous side tissue to lung delay, \( \tau_V \), are around 0.8 [min], 1.0 [min], 1.3 [min] and 1.3 [min], respectively. The comparison of the main results obtained for the two-state and the five-state systems lead to the following conclusion. Both systems show the same behavior, i.e.~increasing transport delay causes loss of stability of the equilibria and occurrence of periodic solutions. The five-state system, however, predicts larger value for the critical transport delay than the two-state system. \section{Conclusions} We have developed computational tools and demonstrated that they are effectively applicable to study stability and bifurcations in two sets of delay equations which have been proposed as possible models of the human respiratory system. The equilibria of the examined systems and their stability were determined. The transport delay was considered as parameter in the two-state model and its critical value was found when the equilibrium loses its stability. This is a supercritical Hopf bifurcation point where a stable periodic solution emerges. The branch of periodic solutions was found (amplitude and frequency obtained for any value of the transport delay). This periodic solution describes the medically important periodic breathing. Periodic breathing occurs when the transport delay exceeds its critical value and the control system reacts to information which no longer describes the state of the system. The five-state model with four delays was examined in three steps. First, the four delays were assumed to take the same value, and the critical delay was almost double of that obtained in the two-state model. Then, the piecewise linear dependence of the parameters on the state variables in the dissociation curves are considered, but simulation results were only slightly different from that obtained for constant parameters. Finally, different values of the delays were taken into account, and simulation results are presented for a set of delays based on transport processes in healthy human bodies. Due to the complexity of the five-state model the only feasible way to conduct a nonlinear stability analysis is by computational means, which shows the significance of the availability of the numerical tools used in this paper. \section{Appendix} \subsection{Nomenclature}\label{symbols} %\begin{itemize} \begin{description} \item[$a, b$] Reciprocal of $p$ and $\sigma$ \item[$B_a$, $m_a$] Constants referring to the relationship between arterial concentration and partial pressure of CO$_2$ \item[$B_v$, $m_v$] Constants referring to the relationship between venous concentration and partial pressure of CO$_2$ \item[$C_{a_{CO_2}}$, $C_{a_{O_2}}$] Arterial CO$_2$ and O$_2$ concentrations \item[$C_{B_{CO_2}}$] CO$_2$ concentration in brain \item[$C_{V_{CO_2}}$, $C_{V_{O_2}}$] Venous CO$_2$ and O$_2$ concentrations, \item[$E_F$] Proportionally constant reflecting the reduction in the alveolar ventilation, \item[$G_C$, $G_P$] Control gain and peripheral control gain \item[$I_C$, $I_P$] Apneic threshold and peripheral cutoff threshold \item[$K_1$, $K_{CO_2}$, $K_{B_{CO_2}}$] Constants referring to the relationship between arterial concentration and partial pressure of CO$_2$ \item[$M_{B_{CO_2}}$, $M_{L_{CO_2}}$, $M_{T_{CO_2}}$] Effective CO$_2$ volume in brain, lung and tissue compartment \item[$M_{L_{O_2}}$, $M_{T_{O_2}}$] Effective O$_2$ volume in lung and tissue compartment \item[$MR_{B_{CO_2}}$, $MR_{T_{CO_2}}$] Metabolic rate for CO$_2$ in brain and tissue compartment \item[$MR_{T_{O_2}}$] Metabolic rate for O$_2$ in tissue compartment \item[$Q$] Volume of blood per unit time \item[$Q_B$, $Q_T$] Blood flow to the brain and tissue \item[$p$] CO$_2$ production rate \item[$P_{a_{CO_2}}$, $P_{a_{O_2}}$] Arterial partial pressures of CO$_2$ and O$_2$ \item[$P_{B_{CO_2}}$] Partial pressure of CO$_2$ in brain \item[$P_{I_{CO_2}}$, $P_{I_{O_2}}$] Inspired partial pressures of CO$_2$ and O$_2$, \item[$P_{V_{CO_2}}$, $P_{V_{O_2}}$] Venous partial pressures of CO$_2$ and O$_2$ \item[$t$ ] Time \item[$V$] Transformed ventilation function of the two-state model \item[$V_I$] Ventilation function of the five-state model \item[$W$] Ventilation function of the two-state model \item[$x$, $y$] Transformed state variables of the two-state model \item[$\tilde x$, $\tilde y$] Arterial CO$_2$ and O$_2$ concentrations, \item[$x_I$, $y_I$] Inspired CO$_2$ and O$_2$ concentrations \item[$\alpha$, $\beta$] Positive constants referring to the diffusibility of CO$_2$ and O$_2$, \item[$\xi$, $\eta$] Perturbation of $x$ and $y$ around their equilibrium \item[$\sigma$] O$_2$ consumption rate \item[$\tau$] Transport delay \item[$\tau_0$ ] Analytically computed critical delay \item[$\tau_{cr}$] Approximate critical delay \item[$\tau_a$] Lung to carotid artery delay \item[$\tau_B$] Lung to brain delay, \item[$\tau_T$] Lung to tissue transport delay \item[$\tau_V$] Venous side transport delay from tissue to lung. \item[$\bar{\quad}$] The superscript $\bar{\quad}$ indicates quantities in equilibrium. %\end{itemize} \end{description} %\subsection*{Appendix: Parameter Values} \begin{table}[htb] \caption{Values of the parameters of the five-state system \label{parameters}} \begin{center} \begin{tabular}{|l|l|l|} \hline Quantity & Unit & Value \\ \hline $P_{I_{CO_2}}$ & mmHg & 0.24 \\ $P_{I_{O_2}} $ & mmHg & 125.64 \\ $M_{L_{CO_2}}$ & l & 3.2 \\ $M_{L_{O_2}}$ & l & 2.5 \\ $M_{B_{CO_2}}$ & l & 0.9 \\ $M_{T_{CO_2}}$ & l & 15.0 \\ $M_{T_{O_2}}$ & l & 6.0 \\ $MR_{B_{CO_2}}$ & l/min & 0.042 \\ $MR_{T_{CO_2}}$ & l/min & 0.235 \\ $MR_{T_{O_2}}$ & l/min & 0.29 \\ $Q$ & l/min & 6 \\ $Q_B$ & mmHg/min & 0.7 \\ $Q_T$ & mmHg/min & 4.0 \\ $K_{CO_2}$ & $l_{STPD}$/(l\ mmHg) & 0.0065 \\ $K_{B_{CO_2}}$ & $l_{STPD}$/(l \ mmHg) & 0.0065 \\ $E_F$ & 1 & 0.8 \\ $m_v$ & $l_{STPD}$/(l \ mmHg) & 0.0021 \\ $B_v$ & $l_{STPD}$/l & 0.0662 \\ $m_a$ & $l_{STPD}$/(l \ mmHg) & 0.00025 \\ $B_a$ & $l_{STPD}$/l & 0.1728 \\ $G_P$ & 1/(min \ mmHg) & 26.5 \\ $I_P$ & mmHg & 35.5 \\ $G_C$ & 1/(min \ mmHg) & 3.2 \\ $I_C$ & mmHg & 35.5 \\ \hline \end{tabular} \end{center} \end{table} \begin{table}[htb] \caption{Dependence of $m_a, B_a, m_v, B_v$ on $P_{a_{O_2}},P_{V_{O_2}}$} \begin{center} \begin{tabular}{|c|c|c|c|c|} \hline & $m_a$ & $B_a$ & $m_v$ & $B_v$ \\ \hline $P_{a_{O_2}} < 55$ & 0.00211 & 0.0662 & & \\ $55 \le P_{a_{O_2}} < 70$ & 0.00067 & 0.1434 & & \\ $70 \le P_{a_{O_2}}$ & 0.00025 & 0.1728 & & \\ $P_{V_{O_2}} < 55$ & & & 0.00211 & 0.0662 \\ $55 \le P_{V_{O_2}} < 70$ & & & 0.00067 & 0.1434 \\ $70 \le P_{V_{O_2}}$ & & & 0.00025 & 0.1728 \\ \hline \end{tabular} \end{center} \end{table} \begin{thebibliography}{99} \bibitem{Bellman} Bellman, R., Cooke, K. L., {\em Differential-Difference Equations}, Academic Press, New York, London, 1963. \bibitem{Carpentier} Carpentier, M. P., Dos Santos, A. F., Solution of Equations Involving Analytic Functions, {\em J. Comput. Phys.}, {\bf 45}, pp. 210-220, 1982. \bibitem{CookeTuri} Cooke, K. L., Turi, J., Stability, Instability in Delay Equations Modeling Human Respiration, {\em Journal of Mathematical Biology}, {\bf 32}, pp. 535-543, 1994. \bibitem{Engelborghs} Engelborghs, K., Luzyanina, T., Samaey, G., \textsc{dde-biftool} v. 2.00: a Matlab Package for Bifurcation Analysis of Delay Differential Equations, {\em Technical Report, TW-330}, Department of Computer Science, K.U.Leuven, Leuven, Belgium, 2001. \bibitem{Hartung} Gy\H ori, I., Hartung, F., Turi, J., On Numerical Solutions for a Class of Nonlinear Delay Equations with Time- and State-Dependent Delays, {\em Proc. of the First World Congress of Nonlinear Analysts}, Tampa, Florida, 1992. \bibitem{Hartungconv} Gy\H ori, I., Hartung, F., Turi, J., Numerical Approximations for a Class of Differential Equations with Time- and State-Dependent Delays, {\em Appl.~Math.~Lett.}, Vol.~8, No.~6, pp.~19-24, 1995. \bibitem{thesis} Koll\'ar, L. E., {\em Numerical Stability Analysis of a Respiratory Control System Model}, MSc thesis, Department of Mathematical Sciences, The University of Texas at Dallas, Richardson, Texas, USA, 2002. \bibitem{Khoo} Khoo, M. C. K., Kronauer, R. E., Strohl, K. P., Slutsky, A. S., Factors Inducing Periodic Breathing in Humans: a General Model, {\em J.~Appl.~Physiol.: Respirat.~Environ.~Exercise Physiol.}, {\bf 53(3)}, pp. 644-659, 1982. \bibitem{Kuhn} Kuhn, H. W., A New Proof of the Fundamental Theorem of Algebra, {\em Math. Programming Stud.}, {\bf 1}, pp. 148-158, 1974. \bibitem{Manitius} Manitius, A., Tran, H., Payre, G., Roy, R., Computation of Eigenvalues Associated with Functional Differential Equations, {\em SIAM J. Sci. Stat. Comput.}, Vol. 8, No. 3, pp. 222-247, 1987. \bibitem{Stech85} Stech, H. W., Hopf Bifurcation Calculations for Functional Differential Equations, {\em Journal of Mathematical Analysis and Applications}, Vol. 109, No. 2, pp. 472-491, 1985. \bibitem{Stech87} Stech, H. W., A Numerical Analysis of the Structure of Periodic Orbits in Autonomous Functional Differential Equations, {\em NATO ASI Series}, Vol. F37, Dynamics of Infinite Dimensional Systems, Springer-Verlag, Berlin, Heidelberg, 1987. \bibitem{West} West, J. B., {\em Respiratory Physiology - the Essentials}, Baltimore: Williams and Wilkins, 1974. %\bibitem{Stepan} St\'ep\'an, G., {\em Retarded Dynamical Systems}, %Longman, Harlow, UK, 1989. \end{thebibliography} \end{document}