\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}