- Email: [email protected]
A mathematical model of the human respiratory system
A MATHEMATICAL MODEL OF THE HUMAN RESPIRATORY SYSTEM W.F. Fincham and F.T. Tehrani
ABSTRACT A model of the human respiratory system is proposed which has a satisfactory performance under different physiological conditions. The model comprises a continuous plant and a discrete controller which generates and updates the drive signal to the plant at the end of every breath to represent the Hering-Breuer reflex. Arterial and central medullary sensors are included. The lung volume,
dead space volume, cardiac output and cerebral blood flow are time varying. The respiratory work is minimized. The model is examined and simulation results of its performance in hypercapnia, hypoxia, periodic breathing and moderate exercise are presented. The responses presented include the relatively fast transients of Cheyne-Stokes breathing and the slower transients associated with carbon dioxide inhalation,
Keywords: Respiratory system, mathematical model, circulation
INTRODUCTION Mathematical modelling of the respiratory system started with the work of Gray which was published in 1945. He formulated pulmonary ventilation in relation to hydrogen ion concentration, and to carbon dioxide and oxygen tensions of arterial blood in the steady-state and in doing so became the first person to provide a mathematical description of the chemical control of ventilation’. Grodins et al. later followed by presenting the first dynamic model of the respiratory system in 19542. After this relatively simple model, the 1960s saw enormous activity in the development of more sophisticated dynamic models. It was during this period that the digital computer began to take over from the analogue machine as the medium for equation solving. References 3 to 11 are a small selection of the papers published by some stalwarts of modelling which include the names of Defares, Grodins, Horgan, Lange, Milhorn, Guyton, Longobardo, Yamamoto and Dickinson. Although these models include relationships of considerable mathematical detail, none consider the respiratory work output and the Hering-Breuer reflex and both the plant and controller act continuously. In almost all the chemical control models, the breathing cycle is ignored and the lungs are regarded as a rigid, continuously ventilated box. A paper by Priban and Fincham12 outlines a complex stochastic model which includes several optimizing loops and has the novel feature of working in breath time, the objective being to produce a scheme which was less constrained by parameter differences between individual experimental subjects. The model still has not been fully tested. The purpose of this paper is to present a sufficiently descriptive model which embodies at least the omisElectrical Engineering Department, Queen Mary College, Mile End Road, London, El 4NS, UK.
o 1983 Butterworth & Co (Publishers) Ltd. 0141~5245/83/020125-09 803.00
sions stated above. It can reproduce both the slower and faster transient effects of different stimuli on blood chemistry. The controller acts discretely while the plant is modelled as a continuous system. The lung volume varies sinusoidally over the breathing cycle, the ventilation being defined by the chemical states of the blood in the neighbourhoods of the peripheral and central receptors. The breathing frequency and the tidal volume are optimized using the criterion of minimal respiratory work and assumes that the dead-space volume can also be controlled. A Hering-Breuer type reflex is included which allows the controller to modify the rate and depth of breathing only at the end of each breath. The blood flow rates in the body tissues and brain are time varying and are defined in terms of the blood chemistry and metabolic rate. Solution of the equations was achieved using the block oriented ‘Interactive Simulation Language’ (ISL) on an ICL 2980 digital computer.
Description of the model Figure 1 shows the system divided into the two major compartments of the plant and the respiratory controller. Information about the arterial and central medullary receptor gas tensions is continuously available to the controller which updates the drive signal to the plant at the end of every breath. The plant consists of the blocks lungs, brain, cerebrospinal fluid and lumped body tissue. Also included are the central receptors of the medulla and peripheral receptors of the carotid body, a blood transport time delay being interposed on gas concentration changes between the lungs and carotid body receptors. The plant includes blood flow controllers which determine the cerebral blood flow and cardiac output in terms of the arterial tensions of oxygen and carbon dioxide. In addition, cardiac output is dependent on the metabolic rate.
J. Biomed. Eng. 1983, Vol. 5, April
125
model: W.F. Fincham and F.T. Tehrani
Respiratory
mspired arterial blood 3",dtj _ driving signal generator f13), !I
If
frequency optimizer
clock 4831
k-4
11
1
1 Q 1 cardiac
P-
ii k
receptor(9)
P
CSFC02
R
8
73 input dynamics
'dynamics continuous plant
discrete controller
B
4 rain metabolic
02, CO2 main blood circulation; - - -. central medullary CO2 effects;
Figure 1 Model block diagram. -, equation - see text.
Plant equations Equations (1) to (6) give the mass balance conditions of the various compartments. The symbols used in the equations are listed in Appendix 2. Alveolar space. Flow rate equation for carbon dioxide assumes homogeneous mixing of venous blood from the brain and body tissues occurs prior to the lungs.
(CVTC02
-
Gco2)
QT
+ (CVBCO~
-
Gco2)
Similarly (GO2
-
FACTl=FACT2=Ofor
Cvro,)
QT
+ (COO,
-
cwo2)
+
FACT
dt
l
where, for inspiration
and FACT
126
lpAC02
-
1
(1)
body tissue
CVTCO~
l
QT
=
Pb -47
2 =
('102
Cacoz
l
(2)
2
CVTO~
l
QT
= Cao, _
s,
l
.
QT
+MRTco,
dCTcoz
l
(3)
dt
QB QT
-MRTO,
(4)
dCTo, dt
: PIco,)
QB
-ST
flA02
1 =
2
The constraint conditions on FACT 1 and FACT 2 assume that the composition of gas entering the alveolar space at the beginning of inspiration is unaffected by the composition of gas in the dead space volume V, except that it is saturated with water vapour.
for oxygen:
= P,--‘47
FACT
+ FACT
dt
) relevant
and for expiration:
Lumped @ACO, =P+.
i nputs
dv dt
l
- pAO2)
Pb -47
J. Biomed. Eng. 1983, Vol. 5, April
.
Brain tissue 1 dv z J
for
CVBCO~
l
(28 =
~~cozQB l
-SB
+MRBCOz
dCBco2 l
-d-F---
(5)
Respiratory model: W.F. Fincham and F.T. Tehrani
cVB02
&
l
=
co’o, - @3 -SB
(the prime denotes
Metabolism
MRBO,
dCBOz
’
dt
a time delayed
quantity).
The usual approximation of the gas tension in tissue being equal to that in the venous blood is made. In addition, it is further assumed that the alveolar and arterial tensions of CO2 are equal but for oxygen, the alveolar tension is 4 mm Hg higher than the arterial tension.
dynamics
The metabolic rates of the body and brain tissues are included in equations (3) to (6) where the rates associated with oxygen and carbon dioxide are defined independently thus allowing selection of a suitable value for the respiratory quotient. An increase in the level of exercise was simulated by injecting a step change at input RTT (the magnitude of thestep defining the final value of MRT) in Figwe I and time shaping this by a first order lag of time constant r3 to produce a relatively smooth increase in tissue metabolism MRT, i.e.
d(M&-)
RTT-MRT
Blood gas relationships
dt=
= CO2
Brain tissue metabolism constant.
(1
K1
_e-K2p02)2
for oxygen l3 and for carbon dioxide a linear best fit is used for a typical working rangei4.
73
MRB was assumed to remain
The total body metabolism
MR.
MR = MRT +MRB GO,
=
K3
l
Pco,
Central receptor
which was used to determine ratio MRR
equations
The carbon dioxide tension transduced by the central receptors Pcco2 is taken to be the CO2 pressure existing at a depth, d, below the surface of the medulla and is a function of surface CO2 tension, i.e. that of the cerebrospinal fluid P~SFCQ, and that of the deep brain tissue PVBCO~15, viz. pcco2
= pVBC02
+ (PCSFC02
exp W(aB
l
-
pVBC02)
1
K4P
(9)
and the diffusion equation for CO2 across the bloodCSF barrier is taken to be: ‘%SFC02
dt
= k
(pVBC02
(11)
-
pCSFC02)
(10)
K, being the CO, diffusion
time constant. This equation assumes a first-order effect and the numerical value of KS was estimated from the change in the l? of cisternal fluid of dogs following a ventilatory induced change in arterial bloodN. It is assumed that the oxygen tension of the blood has no centrally originated influence on ventilation. Blood flow control Equations (1) to (6) and (9) each contain a blood flow dependent parameter. Algebraic relationships giving the dependence of total cardiac output and and brain blood flow have previously been derived in terms of gas tensions and metabolic rate ratio16. The equations are not repeated here because they are too long. However the relationship numbers R 1 to R6 used ini are included in Figure 1. These relationships give steady state values of the blood flow and in this model the controller equations were completed by adding a first-order lag to each. The time constants used were rr and r2, respectively, for cardiac butput and brain blood flow, their values being determined as described below.
MRR = b
the metabolic
rate
MR (actual) asal level of MR
A metabolically related neural drive component the ventilation MR V was defined in Ref. 18. T4
dww
=
dt
(MRR -1)
-MRV
to
(12)
for MRR > 1 and MR V = 0 for MRR < 1, this being the rest condition. It will be helpful for the reader to refer to Figure 1 to follow the argument used to determine the values of rl, TV, r3 and TV. The approximate value for (r3 + r4) = 80 s was taken from [Ref. 17, p. 231 and supported by Figure 5 in Ref. 31. Representing the change in ventilation to exercise by an exponential, the value of 73 was estimated to be 30 s [Ref. 31, Figure 51. However, from Refs 17 and 31 the value of (rr + r3) N 30 s and so r1 was made reasonably small at 3 s. The work reported in Ref. 32 suggests that the value of r2 is < 10 s and although this relates to goats, it seemed not unreasonable to make r2 N r3. Controller equations The blocks comprising the respiratory controller shown in Figure 1 have their functions synchronized by the clock pulse generator. The breathing frequency optimiser defines the time interval between the pulses which : (1)
defines the intervals over which the mean values of the sensor outputs (P,‘02, P& , and Pcco,) are computed according to the general expression
‘LJ
P(t) dt 0
where D is the duration
of the current breath
J. Biomed. Eng. 1983, Vol. 5, April
127
Respiratory
model:
W.F. Fincham and F.T. Tehrani
(2)
causes the estimated demand level of alveolar ventilation v~ to be held for one breath [see equation (13)] , and
(3)
starts the driving signal generator to produce the neural drive at the commencement of each breath such that : l
$
= 7r
where
l
VA
l
sin 27rft
(13)
(P+4r
+cYv,f-“~~
The subsidiary relationship: VA
Rearranging equation (14) in Ref. 20 yields a cubic in the optimum frequency:
4y(md2f3 +7r2v,
f = l/D
It is this synchronizing feature which represents the inflation-deflation reflexes of the mechanoreceptors of the lungs, i.e. the Hering-Breuer reflex.
VE =
et ~1.~~.The relationship considers the lungs elastance oL,airway viscous resistance & airway turbulent resistance 7 and assumes a sinusoidal flow waveform as given by equation (13). Inertia of the lunes is neglected as is the energy dissipated:during”expiration.
+f* VD
enables comparison to be made between experimental and model results.
fA)f’ =
0
(15)
where the dead space volume is given” by: V,
= 0.1698
l&
+ 0.1587
(16)
The numerical prgcess used in the simulation determines V, , given VA using equation (16) and then solves equation (15) for f. The amplitude of the sinusoidal neural drive signal in equation (13) is scaled to give the correct alveolar tidal volume demanded by VA. This procedure is repeated for every breath.
Ventilation control equation This is developed in Appendix 1, and the final form is
A selection of the simulation results obtained under differing test conditions is presented in Figures 2 to 5 and these are discussed in the following section.
tiA&r&)
DISCUSSION
= 0.2025 Fico2+ 0.2332 &co2 + FACT3+MRV-I(
where FACT 3 = 4.72 x lO+ (104 - p&)4*g for p; < 104 mm Hg and FACT 3 = 0 for ii&
2 104 mm Hg
(14A)
This equation gives the demanded level of ventilation as a ratio relative to the resting level [ VA (,,Q] . There are, however, important conditions under which equation (14A) is invalid. In conditions of acute hypoxia, equation (14A) may, because of the associated low Pcoz, produce a negative value which is clearly inadmissible. Immediately this occurs apnoea results i.e. fi&4 = 0
(14B)
Control is returned to equation 14A only when
(1)
Pa~02 reaches a threshold level of 33 mm Hg and is followed by
(2)
the calculated k~ exceeds 0.0225 l/s. It should be noted that because of condition 2 the level of COz. at which breathing recommences is usually between 33 and 35 mm Hg. These conditions are based on the work of Nielsen and Smithlg .
Frequency optimizer The frequency of breathing is calculated on a minimum respiratory work basis as found by Otis
128
J. Biomed. Eng. 1988, VoL 5, April
The plant equation set constitutes a twelfth order highly non-linear model. There are also the dynamics associated with calculating mean values, the discretely operating controller and transport time delays in the blood. A comprehensive set of simulations were performed in order to test the model under different physiological conditions both in transient and steady-state operation and to compare these results, where possible with available experimental data. The block orientated structure of the simulation language ISL made available all the variables used in the above equations but for reasons of brevity only a selection of these is presented in the discussion. The transient responses obtained were for conditions of hypoxia, hypercapnia, hyperventilation and exercise and are represented by the partial set of variables comprising the arterial vslues of and alveolar ventilation VA. Only the p PaCO ‘&%h’ perbreath values are plotted, the detailed ‘within breath’ waveform being nearly sinusoidal. Stability. Instability was never observed in any test condition indicating general stability of the model. Regional stability was demonstrated by letting the simulation run for 30 min with the metabolic rates of CO2 and 0, set at 226 and 267 ml/min, respectively. The repetitive waveforms generated with mean values shown in the top row of Table 1 gave confidence in both the model form and numerical method for solving the equations. Hypoxia. Figure 2 shows the transient responses to 9% oxygen breathing and the return to normal air. Similarly shaped responses occurred at other stimulus levels.
Respiratory model: W.F. Fincham and F.T. Tehrani
Table 1 Steady-state values during hypoxia and hypercapnia Test condition 10 2
ko,
Model response paCO,
pVTCO,
pYBCO,
pCSFCO,
paO,
pkTO,
Experimental ~VBO,
Q
&
f
33.0
5.0
0.75
12.7
6.2
VE
f
valuesaa19
VE
paCO,
Resting conditions 0.21
0
39.3
42.0
43.7
43.6
99
Hypoxic
42.1
6.3
* 0.3
38.6 f 0.6
conditions
0.15
0
38.9
41.2
42.4
42.7
63.1
39.1
32.1
5.5
0.93
12.9
6.4
0.10
0
34.6
36.7
38.0
39.9
43.7
34.5
25.1
6.8
0.98
15.0
8.1
0.08
0
31.2
33.3
34.5
37.7
36.9
31.4
21.6
7.5
1.02
16.6
9.7
Hypercapnic
13.7
conditions
0.21
0.02
40.5
42.9
44.7
44.7
114
42.6
33.6
5.0
0.75
15.9
8.9
15.7
9.1 f 0.3
0.21
0.03
41.5
43.9
45.8
45.7
122
42.9
33.8
5.0
0.75
18.1
11.3
17.3
12.0 f 0.6
42.6 f 0.6
0.21
0.05
46.0
48.0
49.5
49.0
132
46.2
39.3
6.2
0.93
23.6
20.2
20.5
22
46.9 * 0.5
f 1.5
41
* 0.6
MR = 267 ml 0, /min. Q, QB and VE are in l/min. All entries are rounded to three figures.
At the beginning of hypoxia there was a rapid fall in Paoz causing an increase in the cerebral blood flow. This effect in turn caused washout and a decrease in COz tension in the neighbourhood of the central receptors. Consequently, alveolar ventilation fell slightly below normal (about 6%) during the first 30 s. This early decrease in ventilation was followed by a rapid increase which was due to a further and rapid decrease in arterial oxygen tension. The rise in ventilation in turn lead to a delayed decrease in CO, tension. This antagonistic effect produced the overshoot evident in the ‘on’ transient. As an equilibrium state was approached, arterial tensions of COz and O2 were lower than the initial values whereas the ventilation was higher.
o200
403
600
10 min excitation. The results show large percentage changes in VE and cardiac output as the hypoxic effects became more acute, a result consistent with practical observations. On re-exposure to fresh air (10~ = 0.21), Paoz increased and caused a rapid fall in ventilation leading to apnoea because the Pacm was also lower than normal. During apnoea, Paoz decreased and Pacoz increased to reach the threshold conditions given in equation (14). This activity was repeated causing periodic breathing which persisted for several minutes. The transient responses of the model to 4% COz breathing and recovery are shown in Figure 3. During the first minute there was almost a 100% rise in ventilation and a rise of 2 mm Hg in arterial CO? tension. There was a further and slower increase in these values which reflected the effects of storage volumes of carbon dioxide in the body tissues and brain and the long time constant of 320 s associated with the CSF diffusion equation (10). In these results, changes in ventilation lagged behind changes in Pacoz and it was noted in other simulation results that the rapidity of the responses slowed for increased concentration levels of CO, in the inspired gas. Table 1 shows values of model variables after 18 min of hypercapnia and offers a comparison between model and experimental values2’p29.
Hypercapnia.
l.......‘..,.......<
0
Table 1 shows the values of other variables after
‘Em loo0
12w
1400
1600
ItlM)
I
During recovery from hypercapnia, the rates of change of ventilation and Pacoz were initially very high and resulted in a slight undershoot before the variables subsequently returned to normal values. These features have been observed experimentally22*23. The simulation result of inducing periodic breathing by voluntary hyperventilation is shown in Figure 4. After two minutes of hyperventilation at a level of ei ht times normal ventilation the P,co, fell to the Pow value of 19 mm Hg
Hyperventilation. cl
200
400
600 t1nLe
Figure 2
Response
&xl
1000
1200
1400
1600
Iso0
Csecs)
to 9% oxygen
breathing
(mean
values).
J. Biomed. Eng. 1983, Vol. 5, April
129
Respiratory
model:
b
W.F. Fincham and F.T. Tehrani
503
1000
15&I
I ,
4: co.7
2000
2500
recovery
I0.20
-
0.15
I I
0.10
.05
Figure 4 Response to 2 min voluntary hyperventilation of 8 times normal ventilation (mean values).
-
. 0
SC0
1000
I I 1500
2003
1 2500
-‘--v-----
time (secl
Figure 3 values).
Response to 4% carbon dioxide breathing (mean
I
whereas Pao2 rose to 139 mm Hg. Cessation of hyperventilation was followed by apnoea for 145 s after which periodic breathing resulted. The apnoeic oscillations were followed by a few cyclic (non-apnoeic) variations in ventilation as the levels of CO, and 0, tensions of the blood increased. Finally normal breathing was resumed. These results are very similar to the clinical examples obtained by Douglas and Haldane24. Exercise. This was simulated by increasing the metabolic rate of the body tissues over .a period of 150 s. Figure 5 shows that ventilation also rose to three times normal in unison with the increasing metabolic rate. The changes in PaGo and Paoz are interesting because there was an initial slight fall in Pa~2, presumably because oxygen was metabolically consumed at a faster rate than it could be replaced in the arterial blood by ventilation. Matching this was a rise in Pacoz. After about 500 s a near equilibrium state was reached with gas tensions having returned to near their resting levels. These features show general agreement with experimental data25-27. During recovery, ‘opposite’ effects occurred, there being a small increase in Pao2 as the high rate of ventilation temporarily caused over oxygenation of the blood. The period of recovery appeared to be longer than the time taken to reach a steady-state
130 J. Biomed. Eng. 1983, Vol. 5, April
38
. 0
I
,
500
I
I
1000
t
1500
2mo
2500
I 1500
zoo0
, 2500
I
exercise
,
recovery
I
I I
0.25 ' 0.2
.05
.
I 0
I 500
I I 1000 Time
Figure 5 values).
(set)
Response to exercise of 3 times normal (mean
Respiratory
model:
W.F. Fincham and F.T. Tehrani
Table 2 Steady state values during moderate exercise
Test value
MR (MRR)
Experimental
Model response
paCO,
PVTCO,
pVBCO,
pCSFCO,
PaO,
PVTO,
values’”
FVBO,
4
QB
f
VE
a-v0,
a-v0,
267 (I)
39.3
42.0
43.7
43.6
99
42.1
33.0
5.0
0.75
12.7
6.2
5.0
5.0 f 1
544 (2)
39.3
42.8
43.6
43.6
102
35.9
33.2
7.6
0.75
18.0
11.6
6.5
6.5 f 0.6
822 (3)
39.2
43.5
43.5
43.6
101
30.5
33.0
9.5
0.75
21.8
16.4
8.2
7.7 f 0.7
1100 (4)
39.2
44.2
43.5
43.6
100
25.8
32.9
10.8
0.74
24.0
21.3
9.9
9.0 f 1
1377 (5)
39.2
45.0
43.5
43.7
99
21.7
32.7
11.8
0.74
25.8
26.0
11.5
10.3 f 1.3
*0 * =0.21,1~0,
=O,MRmlO,/min;a-v0,
= arterio-venousoxygen difference%. MRR refers to CO, ratios. Q, QB and VE are in I/min.
during the ‘on’ transient, mainly because the perfusion rate of blood in the tissues is dependent on the metabolic activityI and decreases substantially during recovery. The detail of the transient shapes could be further investigated if the values of the metabolic activity time constants 73 and 7q were changed, thus offering considerable scope for matching model and experimental results. Table 2 summarizes the results after 10 min of moderate exercise. It includes experimental values of the arterio-venous gas differences (u - u O2 %) compiled from the literature14 and shows good agreement between model and practice..
In Tables 1 and 2 both arterial and venous blood levels of oxygen and carbon dioxide in different parts of the system have been given. These clearly indicate the success of the control model as a regulator of these variables for the different test conditions even when these are greatly different from the rest condition. Concluding
tion solving, is most conducive to interactive work and so offers painless assessment of the effects of changes in model structure and parameter values.
REFERENCES 1 2
3
4
5
6
thoughts
The detailed refinement of models of complex physiological systems is only worthwhile if reliable experimental results are available. The question of repeatability in physiological testing is to a large extent still unresolved and so a model must, in part at least, relate to an hypothetical normal person. This paper shows the potential of the model to realistically describe the human respiratory responses for a variety of test conditions. The model is relatively sophisticated in that it includes both the peripheral and central receptors and a relatively detailed discrete controllei which functions on the basis of optimising the mechanical work of breathing. The effects of a Hering-Breuer type reflex are embodied to accomplish respiratory synchronization. Included is the autonomic control of cardiac output and cerebral blood flow and breathing is of a cyclical nature as opposed to that of a constant through-flow type. The long time constant associated with carbon dioxide diffusion in the cerebrospinal fluid adds a distinctive element in the response shapes. The shaping of the ‘on’ and ‘off’ transients of metabolic rate certainly needs further investigation to determine how response shapes change. The use of ISL as the simulation medium, although not the most efficient computer method for equa-
7
8
9
10
11
I2
13 14 15
Gray, J.S. in Pulmonary
ventilation and its physiological regulation, C.C. Thomas, Illinois, 1949. Grodins, F.S., Gray, J.S., Schroeder, K.R., Norins, A.L., and Jones, R.W., Respiratory responses to CO2 inhalation. A theoretical study of a non-linear biological regulator, J. Appl. Physiol. 1954, 7,283~308. Defares, J.G., Derksen, H.E. and Duyff, J.W., Cerebral blood flow in the regulation of respiration, Acta. Physiol. Pharmacol. Neerkandica, 1960.9.327-360. Horgan, J.D., and Lange, R.L., Digital computer simulation of the human respiratory system, IEEE Int. Conv. Rec., 1963, part 9, 149-157. Milhorn, H.T., Benton, R., Ross, R. and Guyton, A.C., A mathematical model of the human respiratory control system Biophys. J. 1965,5,27-46. Longobardo, G.S., Chcrniack, N.S., and Fishman, A.P., CheyneStokes breathing produced by a model of the human respiratory system,]. Appl. Physiol. 1966, 21(6), 1839-1846. Grodins, F.S., Buell, J. and Bart, A J., Mathematical analysis and digital simulation of the respiratory control system, J. Appl. Physiol. 1967, 22, 260-275. Yamamoto, W. and Hori, T., Phasic air movement model of respiratory regulation of carbon dioxide balance, Comput. Biomed. Res. 1971,3,699-717. Milhorn, H.T., Reynolds, W J. and Holloman, G.H., Digital simulation of the ventilatory responses to CO2 inhalation and CSF perfusion, Comput. Biomed. Res. 1972,5,301-314. Damokosh-Giordano, A., Longobardo, G.S. and Cherniack, N.S., The effects of controlled system (plant) dynamics on vent&tory responses to disturbances in CO2 balance, in Regulation and Control in physiologiCa systems, (Eds. Iberall, A.S. and Guyton, AX.) Instrument Sot Amer, Pittsburgh, 1973. Dickinson, C J., A computer model of human respiration, MTP Press Ltd. England, 1977. Priban, I.P., and Fincham, W.F., Self-adaptive control and the respiratory system, Nature, 1965, 208, 339-343. Visser, B.F., Pulmonary diffudon of carbon dioxide Phys. Med. Biol. 1960-61,5 (2), 155-166. Ruth, T.C. and Patton, H.D., (Eds.) in Physiology and Biophysics, W.B. Saunders, 1966, p. 766. Mitchell, R.A., Loeschcke, H.H., Massion, W.H. and Severinghaus, J.W., Respiratory responses mediated through superficial chemosensitive areas on the medulla,j. Appf. Physiol. 1963,18,523-533.
J. Biomed. Eng. 1983, Vol. 5, April
131
Respiratory
16
17
18
19
20 21
22
2.3
24
25
26
27 28
29
30
31
32
model:
W.F. Fincham and F.T. Tehrani
Fir&am, W.F. and Tehrani, F.T., On the regulation of cardiac output and cerebral blood flow. J Biomed. Eng. 1983,5,73-75. Wigertz, O., Dynamics of respiratory and circulatory adaptation to muscular exercise in man, Acta Physiol. Stand. suppl. 1971,363. Tehrani, F.T., Dynamic modelling of the human respiratory system, Ph. D. thesis, 198 1, University of London. Nielsen, M. and Smith, H., Studies on the regulation of respiration in acute hypoxia, Acta Physiol. Stand. 1952, 24, 293513. Otis, A.B ., Fenn, W.O. and Rahn, H., Mechanics of breathing in man, J. Appl. Physiol. 1950,2, 592-607. Lambertsen, C J., Factors in the simulation of respiration by carbon dioxide, in The regulation of human respiration, (Eds. D J.C. Cunningham and B.B. Lloyd) Symp. Proc., 1963, p 257-276, Blackwell. Nielsen, M., Unter suchungen uber die Atem regulation beim Menschen, Skandinav. Arch. Physiol. 1936,74 (suppl. 10) p 87-208. Reynolds, W J., Milhorn, H.T. and Holloman, G.H., Transient vent&tory responses to graded hypercapnia in man, J. Appl. Physiol. 1972,33(l), 47-54. Douglas, C.G., and Haldane, J.S., The causes of periodic or Cheyne-Stokes breathing, J. Physiol. (Lond.) 1909,38,401-419. Astrand, P.O., Human physical fitness with special reference to sex and age, Physiol. Rev. 1956,36, 307-329. Wasserman, K., Van Kessel, A.L. and Burton, G.G., Interaction of physiological mechanisms during exercise, J. Appl. Physiol. 1967,22, 71-85. Comroe, J.H., The hyperpnea of muscular exercise, Physiol. Rev. 1944,24,319-333. Hey, E., Lloyd, B.B., Cunningham, D J.C., Jukes, M.G.M., and Bolton, D.P.G., Effects of various respiratory stimuli on the depth and frequency of breathing in man, Respir. Physiol. 1966, 1, 193-205. Lambertsen, C J., Carbon dioxide and respiration in acid-base homeostatis, Anesthesiology, 1960, 21, 642-651. Brooks, C.M., Kao, F.F. and Lloyd, B.B. (Eds.), Cerebro spinalfluid and the regulation of ventilation, Oxford, Blackwell, 1965. (see Dynamic response characteristics of several C&-reactive components of the respiratory control system, by Lambertson, C J., Gelfand, R. and Kemp, R.A., p. 2 11-245). Cerretelli, P., Sikand, R. and Farhi, LE., Readjustments in cardiac output and gas change during onset of exercise and recovery, J. Appl. Physiol. 1966, 21(4), 1345-1350. Doblar, D.D., Min, B.G., Chapman, R.W., Harback, E.R., Welkowitz, W. and Edelman, N.H., Dynamic characteristics of cerebral blood flow response to sinusoidal hypoxia, J. Appl. Physiol. 1979,46(4), 721-729.
APPENDIX Derivation
(Al) is extended by assuming that half of the ventilatory response to CO3 and ti results from arterial blood via the peripheral receptors and half from the CSF21. Hence ,
A VR = 0.11 e
+ 0.131
,
Pacoz + 0.11 H&F
of the ventilation
(A2)
+ 0.131 PCSFCQ + FACT 3 -K
Using Gray’s relation (1) between arterial ti and levels of hypoxia and hypercapnia.
P a~~z in different
flu’ = 0.65Pacol + 13.5
(9
Similarly, writing the CSF hydrogen ion concentration:
6.h WCSF
=
(BHC%)CSF
(A4)
’ PCSFCOz
Substitution of equations (A3) and (A4) into (A2), and assuming that (BHCO~)CSF’ 6 and X are constant yields : A VR = 0.2025
&Oz
+ 0.2332
Pcs~co
2 (A5)
+FACT3-K
where the parameter values given in appendix 2 have been used. The final step in this development assumes that PCSFCO? can be replaced by the CO, tension in the neighbourhood of the medullary receptors Pccoz (equation 9) and that metabolic activity also has a relatively direct influence on AVR [equation (12)].
APPENDIX
2
Glossary of symbols and values The main symbols denote the variable or parameter and subscripts indicate the model location with reference to a specific chemical, as appropriate. X
denotes the mean value of x
X’
denotes the time delayed version of x
Model variables and symbols Symbol
Definition
Value/units
AVR
Alveolar ventilation
BHCO,
Bicarbonate
ratio
concentration
0.58
[WPWI
1
C
equation
On the basis of the multiple factor theory of Gray’, the alveolar ventilation ratio A VR is: A VR = 0.22 H+ + 0.262 Pacq
+ FACT 3 -K (Al)
where FACT 3 is the response of the arterial chemoreceptors to hypoxia and K is a constant.
132
this paper, equation
J. Biomed. Eng. 1983, Vol. 5, April
In
Volume concentration gas in the blood
of
[l(STPWl
D
Breath duration
(s)
f
Frequency
of breathing
(breath/s)
F
Fractional of gas
composition
H+
Hydrogen ion concentration
I
Inspired gas concentration (fractional volume)
(nmoles/l)
Respiratory
MRB
MRT
model:
W.F. Fincham and F. 1: Tehrani
Metabolic rate in brain
P/s)
basal rates
0,
(0.000925)
Equivalent gas storage space
co2
(0.0009)
ST (tissue)
Metabolic rate in body tissue
(l/s)
basal rates
(0.00352)
0, co2
s
(0.00287)
S, (brain)
(50 1) (1.1 I)
T
Transport delay: arterial blood
(IO s)
(IL
Respiratory system elastance
MRR
Metabolic rate ratio
MRV
Metabolic neural drive factor
P
Air viscosity resistance
P
Partial pressure or gas pressure
Y
Airway turbulent resistance
tpb
Q
- 47)
t . vA
V
(mmw (mmHiid
Blood flow rate
PI4
basal values QT (tissue)
(0.07083)
(2s (brain) RTT
[MR(actual)/ MR(basal)]
Barometric pressure less water vapour pressure at body temperature
H, O/( 1 (BTPS)s)’ ]
Solubility factor (CO2 in CSF)
h
Carbonic acid dissociation constant Time constants:
(795 nmoles/l)
71
Cardiac output controller
(3 s)
6.783 x 10-4) i I(STPD)/ l/mm &I
(l/s)
72
Brain blood flow controller (3 s)
Time
(s)
73
Exercise metabolism dynamic
(30 s)
74
Metabolically derived neural drive
(50 s)
Alveolar ventilation (mean value per breath) basal value VA (rest)
(l/s) (0.0673
Volume
(1)
Depth of central receptor below surface of the medulla
Ventilation controller constant blood gas dissociation constants :
Kl
Oxygen
K2
Oxygen
K3
Carbon dioxide (dynamic)
K4
Central receptor constant
KS
[ 1.285 cm
Metabolic input function to tissue
I/s)
CO, diffusion time constant
Subscr$ ts a
Arterial blood
A
Alveolar gas
b
Barometric condition
B
Brain
c
A point in the medulla
(332
Carbon dioxide
(17.4)
CSF
Cerebrospinal fluid
(0.2) [0.046 (mm Hg)-’ I [0.016 (mm Hg)-’ 3
D
Dead space
E
External ventilation
I
Inspired gas
02
Oxygen
T
Body tissue (excluding brain)
V
Venous blood
VB
Brain venous blood
VT
Tissue venous blood
(15 x 10m3cm)
FACT N Factor N as specified in text K
(8.55 cm Hz O/l) (3.1 cm H, O/l (BTPS)/s)
6
(0.0125)
Model fmumeters d
(1)
(346 x lo3 sec. cme2 .I-’ ) (320 s)
J. Biomed. Eng. 1983, Vol. 5, April
133
ABSTRACT A model of the human respiratory system is proposed which has a satisfactory performance under different physiological conditions. The model comprises a continuous plant and a discrete controller which generates and updates the drive signal to the plant at the end of every breath to represent the Hering-Breuer reflex. Arterial and central medullary sensors are included. The lung volume,
dead space volume, cardiac output and cerebral blood flow are time varying. The respiratory work is minimized. The model is examined and simulation results of its performance in hypercapnia, hypoxia, periodic breathing and moderate exercise are presented. The responses presented include the relatively fast transients of Cheyne-Stokes breathing and the slower transients associated with carbon dioxide inhalation,
Keywords: Respiratory system, mathematical model, circulation
INTRODUCTION Mathematical modelling of the respiratory system started with the work of Gray which was published in 1945. He formulated pulmonary ventilation in relation to hydrogen ion concentration, and to carbon dioxide and oxygen tensions of arterial blood in the steady-state and in doing so became the first person to provide a mathematical description of the chemical control of ventilation’. Grodins et al. later followed by presenting the first dynamic model of the respiratory system in 19542. After this relatively simple model, the 1960s saw enormous activity in the development of more sophisticated dynamic models. It was during this period that the digital computer began to take over from the analogue machine as the medium for equation solving. References 3 to 11 are a small selection of the papers published by some stalwarts of modelling which include the names of Defares, Grodins, Horgan, Lange, Milhorn, Guyton, Longobardo, Yamamoto and Dickinson. Although these models include relationships of considerable mathematical detail, none consider the respiratory work output and the Hering-Breuer reflex and both the plant and controller act continuously. In almost all the chemical control models, the breathing cycle is ignored and the lungs are regarded as a rigid, continuously ventilated box. A paper by Priban and Fincham12 outlines a complex stochastic model which includes several optimizing loops and has the novel feature of working in breath time, the objective being to produce a scheme which was less constrained by parameter differences between individual experimental subjects. The model still has not been fully tested. The purpose of this paper is to present a sufficiently descriptive model which embodies at least the omisElectrical Engineering Department, Queen Mary College, Mile End Road, London, El 4NS, UK.
o 1983 Butterworth & Co (Publishers) Ltd. 0141~5245/83/020125-09 803.00
sions stated above. It can reproduce both the slower and faster transient effects of different stimuli on blood chemistry. The controller acts discretely while the plant is modelled as a continuous system. The lung volume varies sinusoidally over the breathing cycle, the ventilation being defined by the chemical states of the blood in the neighbourhoods of the peripheral and central receptors. The breathing frequency and the tidal volume are optimized using the criterion of minimal respiratory work and assumes that the dead-space volume can also be controlled. A Hering-Breuer type reflex is included which allows the controller to modify the rate and depth of breathing only at the end of each breath. The blood flow rates in the body tissues and brain are time varying and are defined in terms of the blood chemistry and metabolic rate. Solution of the equations was achieved using the block oriented ‘Interactive Simulation Language’ (ISL) on an ICL 2980 digital computer.
Description of the model Figure 1 shows the system divided into the two major compartments of the plant and the respiratory controller. Information about the arterial and central medullary receptor gas tensions is continuously available to the controller which updates the drive signal to the plant at the end of every breath. The plant consists of the blocks lungs, brain, cerebrospinal fluid and lumped body tissue. Also included are the central receptors of the medulla and peripheral receptors of the carotid body, a blood transport time delay being interposed on gas concentration changes between the lungs and carotid body receptors. The plant includes blood flow controllers which determine the cerebral blood flow and cardiac output in terms of the arterial tensions of oxygen and carbon dioxide. In addition, cardiac output is dependent on the metabolic rate.
J. Biomed. Eng. 1983, Vol. 5, April
125
model: W.F. Fincham and F.T. Tehrani
Respiratory
mspired arterial blood 3",dtj _ driving signal generator f13), !I
If
frequency optimizer
clock 4831
k-4
11
1
1 Q 1 cardiac
P-
ii k
receptor(9)
P
CSFC02
R
8
73 input dynamics
'dynamics continuous plant
discrete controller
B
4 rain metabolic
02, CO2 main blood circulation; - - -. central medullary CO2 effects;
Figure 1 Model block diagram. -, equation - see text.
Plant equations Equations (1) to (6) give the mass balance conditions of the various compartments. The symbols used in the equations are listed in Appendix 2. Alveolar space. Flow rate equation for carbon dioxide assumes homogeneous mixing of venous blood from the brain and body tissues occurs prior to the lungs.
(CVTC02
-
Gco2)
QT
+ (CVBCO~
-
Gco2)
Similarly (GO2
-
FACTl=FACT2=Ofor
Cvro,)
QT
+ (COO,
-
cwo2)
+
FACT
dt
l
where, for inspiration
and FACT
126
lpAC02
-
1
(1)
body tissue
CVTCO~
l
QT
=
Pb -47
2 =
('102
Cacoz
l
(2)
2
CVTO~
l
QT
= Cao, _
s,
l
.
QT
+MRTco,
dCTcoz
l
(3)
dt
QB QT
-MRTO,
(4)
dCTo, dt
: PIco,)
QB
-ST
flA02
1 =
2
The constraint conditions on FACT 1 and FACT 2 assume that the composition of gas entering the alveolar space at the beginning of inspiration is unaffected by the composition of gas in the dead space volume V, except that it is saturated with water vapour.
for oxygen:
= P,--‘47
FACT
+ FACT
dt
) relevant
and for expiration:
Lumped @ACO, =P+.
i nputs
dv dt
l
- pAO2)
Pb -47
J. Biomed. Eng. 1983, Vol. 5, April
.
Brain tissue 1 dv z J
for
CVBCO~
l
(28 =
~~cozQB l
-SB
+MRBCOz
dCBco2 l
-d-F---
(5)
Respiratory model: W.F. Fincham and F.T. Tehrani
cVB02
&
l
=
co’o, - @3 -SB
(the prime denotes
Metabolism
MRBO,
dCBOz
’
dt
a time delayed
quantity).
The usual approximation of the gas tension in tissue being equal to that in the venous blood is made. In addition, it is further assumed that the alveolar and arterial tensions of CO2 are equal but for oxygen, the alveolar tension is 4 mm Hg higher than the arterial tension.
dynamics
The metabolic rates of the body and brain tissues are included in equations (3) to (6) where the rates associated with oxygen and carbon dioxide are defined independently thus allowing selection of a suitable value for the respiratory quotient. An increase in the level of exercise was simulated by injecting a step change at input RTT (the magnitude of thestep defining the final value of MRT) in Figwe I and time shaping this by a first order lag of time constant r3 to produce a relatively smooth increase in tissue metabolism MRT, i.e.
d(M&-)
RTT-MRT
Blood gas relationships
dt=
= CO2
Brain tissue metabolism constant.
(1
K1
_e-K2p02)2
for oxygen l3 and for carbon dioxide a linear best fit is used for a typical working rangei4.
73
MRB was assumed to remain
The total body metabolism
MR.
MR = MRT +MRB GO,
=
K3
l
Pco,
Central receptor
which was used to determine ratio MRR
equations
The carbon dioxide tension transduced by the central receptors Pcco2 is taken to be the CO2 pressure existing at a depth, d, below the surface of the medulla and is a function of surface CO2 tension, i.e. that of the cerebrospinal fluid P~SFCQ, and that of the deep brain tissue PVBCO~15, viz. pcco2
= pVBC02
+ (PCSFC02
exp W(aB
l
-
pVBC02)
1
K4P
(9)
and the diffusion equation for CO2 across the bloodCSF barrier is taken to be: ‘%SFC02
dt
= k
(pVBC02
(11)
-
pCSFC02)
(10)
K, being the CO, diffusion
time constant. This equation assumes a first-order effect and the numerical value of KS was estimated from the change in the l? of cisternal fluid of dogs following a ventilatory induced change in arterial bloodN. It is assumed that the oxygen tension of the blood has no centrally originated influence on ventilation. Blood flow control Equations (1) to (6) and (9) each contain a blood flow dependent parameter. Algebraic relationships giving the dependence of total cardiac output and and brain blood flow have previously been derived in terms of gas tensions and metabolic rate ratio16. The equations are not repeated here because they are too long. However the relationship numbers R 1 to R6 used ini are included in Figure 1. These relationships give steady state values of the blood flow and in this model the controller equations were completed by adding a first-order lag to each. The time constants used were rr and r2, respectively, for cardiac butput and brain blood flow, their values being determined as described below.
MRR = b
the metabolic
rate
MR (actual) asal level of MR
A metabolically related neural drive component the ventilation MR V was defined in Ref. 18. T4
dww
=
dt
(MRR -1)
-MRV
to
(12)
for MRR > 1 and MR V = 0 for MRR < 1, this being the rest condition. It will be helpful for the reader to refer to Figure 1 to follow the argument used to determine the values of rl, TV, r3 and TV. The approximate value for (r3 + r4) = 80 s was taken from [Ref. 17, p. 231 and supported by Figure 5 in Ref. 31. Representing the change in ventilation to exercise by an exponential, the value of 73 was estimated to be 30 s [Ref. 31, Figure 51. However, from Refs 17 and 31 the value of (rr + r3) N 30 s and so r1 was made reasonably small at 3 s. The work reported in Ref. 32 suggests that the value of r2 is < 10 s and although this relates to goats, it seemed not unreasonable to make r2 N r3. Controller equations The blocks comprising the respiratory controller shown in Figure 1 have their functions synchronized by the clock pulse generator. The breathing frequency optimiser defines the time interval between the pulses which : (1)
defines the intervals over which the mean values of the sensor outputs (P,‘02, P& , and Pcco,) are computed according to the general expression
‘LJ
P(t) dt 0
where D is the duration
of the current breath
J. Biomed. Eng. 1983, Vol. 5, April
127
Respiratory
model:
W.F. Fincham and F.T. Tehrani
(2)
causes the estimated demand level of alveolar ventilation v~ to be held for one breath [see equation (13)] , and
(3)
starts the driving signal generator to produce the neural drive at the commencement of each breath such that : l
$
= 7r
where
l
VA
l
sin 27rft
(13)
(P+4r
+cYv,f-“~~
The subsidiary relationship: VA
Rearranging equation (14) in Ref. 20 yields a cubic in the optimum frequency:
4y(md2f3 +7r2v,
f = l/D
It is this synchronizing feature which represents the inflation-deflation reflexes of the mechanoreceptors of the lungs, i.e. the Hering-Breuer reflex.
VE =
et ~1.~~.The relationship considers the lungs elastance oL,airway viscous resistance & airway turbulent resistance 7 and assumes a sinusoidal flow waveform as given by equation (13). Inertia of the lunes is neglected as is the energy dissipated:during”expiration.
+f* VD
enables comparison to be made between experimental and model results.
fA)f’ =
0
(15)
where the dead space volume is given” by: V,
= 0.1698
l&
+ 0.1587
(16)
The numerical prgcess used in the simulation determines V, , given VA using equation (16) and then solves equation (15) for f. The amplitude of the sinusoidal neural drive signal in equation (13) is scaled to give the correct alveolar tidal volume demanded by VA. This procedure is repeated for every breath.
Ventilation control equation This is developed in Appendix 1, and the final form is
A selection of the simulation results obtained under differing test conditions is presented in Figures 2 to 5 and these are discussed in the following section.
tiA&r&)
DISCUSSION
= 0.2025 Fico2+ 0.2332 &co2 + FACT3+MRV-I(
where FACT 3 = 4.72 x lO+ (104 - p&)4*g for p; < 104 mm Hg and FACT 3 = 0 for ii&
2 104 mm Hg
(14A)
This equation gives the demanded level of ventilation as a ratio relative to the resting level [ VA (,,Q] . There are, however, important conditions under which equation (14A) is invalid. In conditions of acute hypoxia, equation (14A) may, because of the associated low Pcoz, produce a negative value which is clearly inadmissible. Immediately this occurs apnoea results i.e. fi&4 = 0
(14B)
Control is returned to equation 14A only when
(1)
Pa~02 reaches a threshold level of 33 mm Hg and is followed by
(2)
the calculated k~ exceeds 0.0225 l/s. It should be noted that because of condition 2 the level of COz. at which breathing recommences is usually between 33 and 35 mm Hg. These conditions are based on the work of Nielsen and Smithlg .
Frequency optimizer The frequency of breathing is calculated on a minimum respiratory work basis as found by Otis
128
J. Biomed. Eng. 1988, VoL 5, April
The plant equation set constitutes a twelfth order highly non-linear model. There are also the dynamics associated with calculating mean values, the discretely operating controller and transport time delays in the blood. A comprehensive set of simulations were performed in order to test the model under different physiological conditions both in transient and steady-state operation and to compare these results, where possible with available experimental data. The block orientated structure of the simulation language ISL made available all the variables used in the above equations but for reasons of brevity only a selection of these is presented in the discussion. The transient responses obtained were for conditions of hypoxia, hypercapnia, hyperventilation and exercise and are represented by the partial set of variables comprising the arterial vslues of and alveolar ventilation VA. Only the p PaCO ‘&%h’ perbreath values are plotted, the detailed ‘within breath’ waveform being nearly sinusoidal. Stability. Instability was never observed in any test condition indicating general stability of the model. Regional stability was demonstrated by letting the simulation run for 30 min with the metabolic rates of CO2 and 0, set at 226 and 267 ml/min, respectively. The repetitive waveforms generated with mean values shown in the top row of Table 1 gave confidence in both the model form and numerical method for solving the equations. Hypoxia. Figure 2 shows the transient responses to 9% oxygen breathing and the return to normal air. Similarly shaped responses occurred at other stimulus levels.
Respiratory model: W.F. Fincham and F.T. Tehrani
Table 1 Steady-state values during hypoxia and hypercapnia Test condition 10 2
ko,
Model response paCO,
pVTCO,
pYBCO,
pCSFCO,
paO,
pkTO,
Experimental ~VBO,
Q
&
f
33.0
5.0
0.75
12.7
6.2
VE
f
valuesaa19
VE
paCO,
Resting conditions 0.21
0
39.3
42.0
43.7
43.6
99
Hypoxic
42.1
6.3
* 0.3
38.6 f 0.6
conditions
0.15
0
38.9
41.2
42.4
42.7
63.1
39.1
32.1
5.5
0.93
12.9
6.4
0.10
0
34.6
36.7
38.0
39.9
43.7
34.5
25.1
6.8
0.98
15.0
8.1
0.08
0
31.2
33.3
34.5
37.7
36.9
31.4
21.6
7.5
1.02
16.6
9.7
Hypercapnic
13.7
conditions
0.21
0.02
40.5
42.9
44.7
44.7
114
42.6
33.6
5.0
0.75
15.9
8.9
15.7
9.1 f 0.3
0.21
0.03
41.5
43.9
45.8
45.7
122
42.9
33.8
5.0
0.75
18.1
11.3
17.3
12.0 f 0.6
42.6 f 0.6
0.21
0.05
46.0
48.0
49.5
49.0
132
46.2
39.3
6.2
0.93
23.6
20.2
20.5
22
46.9 * 0.5
f 1.5
41
* 0.6
MR = 267 ml 0, /min. Q, QB and VE are in l/min. All entries are rounded to three figures.
At the beginning of hypoxia there was a rapid fall in Paoz causing an increase in the cerebral blood flow. This effect in turn caused washout and a decrease in COz tension in the neighbourhood of the central receptors. Consequently, alveolar ventilation fell slightly below normal (about 6%) during the first 30 s. This early decrease in ventilation was followed by a rapid increase which was due to a further and rapid decrease in arterial oxygen tension. The rise in ventilation in turn lead to a delayed decrease in CO, tension. This antagonistic effect produced the overshoot evident in the ‘on’ transient. As an equilibrium state was approached, arterial tensions of COz and O2 were lower than the initial values whereas the ventilation was higher.
o200
403
600
10 min excitation. The results show large percentage changes in VE and cardiac output as the hypoxic effects became more acute, a result consistent with practical observations. On re-exposure to fresh air (10~ = 0.21), Paoz increased and caused a rapid fall in ventilation leading to apnoea because the Pacm was also lower than normal. During apnoea, Paoz decreased and Pacoz increased to reach the threshold conditions given in equation (14). This activity was repeated causing periodic breathing which persisted for several minutes. The transient responses of the model to 4% COz breathing and recovery are shown in Figure 3. During the first minute there was almost a 100% rise in ventilation and a rise of 2 mm Hg in arterial CO? tension. There was a further and slower increase in these values which reflected the effects of storage volumes of carbon dioxide in the body tissues and brain and the long time constant of 320 s associated with the CSF diffusion equation (10). In these results, changes in ventilation lagged behind changes in Pacoz and it was noted in other simulation results that the rapidity of the responses slowed for increased concentration levels of CO, in the inspired gas. Table 1 shows values of model variables after 18 min of hypercapnia and offers a comparison between model and experimental values2’p29.
Hypercapnia.
l.......‘..,.......<
0
Table 1 shows the values of other variables after
‘Em loo0
12w
1400
1600
ItlM)
I
During recovery from hypercapnia, the rates of change of ventilation and Pacoz were initially very high and resulted in a slight undershoot before the variables subsequently returned to normal values. These features have been observed experimentally22*23. The simulation result of inducing periodic breathing by voluntary hyperventilation is shown in Figure 4. After two minutes of hyperventilation at a level of ei ht times normal ventilation the P,co, fell to the Pow value of 19 mm Hg
Hyperventilation. cl
200
400
600 t1nLe
Figure 2
Response
&xl
1000
1200
1400
1600
Iso0
Csecs)
to 9% oxygen
breathing
(mean
values).
J. Biomed. Eng. 1983, Vol. 5, April
129
Respiratory
model:
b
W.F. Fincham and F.T. Tehrani
503
1000
15&I
I ,
4: co.7
2000
2500
recovery
I0.20
-
0.15
I I
0.10
.05
Figure 4 Response to 2 min voluntary hyperventilation of 8 times normal ventilation (mean values).
-
. 0
SC0
1000
I I 1500
2003
1 2500
-‘--v-----
time (secl
Figure 3 values).
Response to 4% carbon dioxide breathing (mean
I
whereas Pao2 rose to 139 mm Hg. Cessation of hyperventilation was followed by apnoea for 145 s after which periodic breathing resulted. The apnoeic oscillations were followed by a few cyclic (non-apnoeic) variations in ventilation as the levels of CO, and 0, tensions of the blood increased. Finally normal breathing was resumed. These results are very similar to the clinical examples obtained by Douglas and Haldane24. Exercise. This was simulated by increasing the metabolic rate of the body tissues over .a period of 150 s. Figure 5 shows that ventilation also rose to three times normal in unison with the increasing metabolic rate. The changes in PaGo and Paoz are interesting because there was an initial slight fall in Pa~2, presumably because oxygen was metabolically consumed at a faster rate than it could be replaced in the arterial blood by ventilation. Matching this was a rise in Pacoz. After about 500 s a near equilibrium state was reached with gas tensions having returned to near their resting levels. These features show general agreement with experimental data25-27. During recovery, ‘opposite’ effects occurred, there being a small increase in Pao2 as the high rate of ventilation temporarily caused over oxygenation of the blood. The period of recovery appeared to be longer than the time taken to reach a steady-state
130 J. Biomed. Eng. 1983, Vol. 5, April
38
. 0
I
,
500
I
I
1000
t
1500
2mo
2500
I 1500
zoo0
, 2500
I
exercise
,
recovery
I
I I
0.25 ' 0.2
.05
.
I 0
I 500
I I 1000 Time
Figure 5 values).
(set)
Response to exercise of 3 times normal (mean
Respiratory
model:
W.F. Fincham and F.T. Tehrani
Table 2 Steady state values during moderate exercise
Test value
MR (MRR)
Experimental
Model response
paCO,
PVTCO,
pVBCO,
pCSFCO,
PaO,
PVTO,
values’”
FVBO,
4
QB
f
VE
a-v0,
a-v0,
267 (I)
39.3
42.0
43.7
43.6
99
42.1
33.0
5.0
0.75
12.7
6.2
5.0
5.0 f 1
544 (2)
39.3
42.8
43.6
43.6
102
35.9
33.2
7.6
0.75
18.0
11.6
6.5
6.5 f 0.6
822 (3)
39.2
43.5
43.5
43.6
101
30.5
33.0
9.5
0.75
21.8
16.4
8.2
7.7 f 0.7
1100 (4)
39.2
44.2
43.5
43.6
100
25.8
32.9
10.8
0.74
24.0
21.3
9.9
9.0 f 1
1377 (5)
39.2
45.0
43.5
43.7
99
21.7
32.7
11.8
0.74
25.8
26.0
11.5
10.3 f 1.3
*0 * =0.21,1~0,
=O,MRmlO,/min;a-v0,
= arterio-venousoxygen difference%. MRR refers to CO, ratios. Q, QB and VE are in I/min.
during the ‘on’ transient, mainly because the perfusion rate of blood in the tissues is dependent on the metabolic activityI and decreases substantially during recovery. The detail of the transient shapes could be further investigated if the values of the metabolic activity time constants 73 and 7q were changed, thus offering considerable scope for matching model and experimental results. Table 2 summarizes the results after 10 min of moderate exercise. It includes experimental values of the arterio-venous gas differences (u - u O2 %) compiled from the literature14 and shows good agreement between model and practice..
In Tables 1 and 2 both arterial and venous blood levels of oxygen and carbon dioxide in different parts of the system have been given. These clearly indicate the success of the control model as a regulator of these variables for the different test conditions even when these are greatly different from the rest condition. Concluding
tion solving, is most conducive to interactive work and so offers painless assessment of the effects of changes in model structure and parameter values.
REFERENCES 1 2
3
4
5
6
thoughts
The detailed refinement of models of complex physiological systems is only worthwhile if reliable experimental results are available. The question of repeatability in physiological testing is to a large extent still unresolved and so a model must, in part at least, relate to an hypothetical normal person. This paper shows the potential of the model to realistically describe the human respiratory responses for a variety of test conditions. The model is relatively sophisticated in that it includes both the peripheral and central receptors and a relatively detailed discrete controllei which functions on the basis of optimising the mechanical work of breathing. The effects of a Hering-Breuer type reflex are embodied to accomplish respiratory synchronization. Included is the autonomic control of cardiac output and cerebral blood flow and breathing is of a cyclical nature as opposed to that of a constant through-flow type. The long time constant associated with carbon dioxide diffusion in the cerebrospinal fluid adds a distinctive element in the response shapes. The shaping of the ‘on’ and ‘off’ transients of metabolic rate certainly needs further investigation to determine how response shapes change. The use of ISL as the simulation medium, although not the most efficient computer method for equa-
7
8
9
10
11
I2
13 14 15
Gray, J.S. in Pulmonary
ventilation and its physiological regulation, C.C. Thomas, Illinois, 1949. Grodins, F.S., Gray, J.S., Schroeder, K.R., Norins, A.L., and Jones, R.W., Respiratory responses to CO2 inhalation. A theoretical study of a non-linear biological regulator, J. Appl. Physiol. 1954, 7,283~308. Defares, J.G., Derksen, H.E. and Duyff, J.W., Cerebral blood flow in the regulation of respiration, Acta. Physiol. Pharmacol. Neerkandica, 1960.9.327-360. Horgan, J.D., and Lange, R.L., Digital computer simulation of the human respiratory system, IEEE Int. Conv. Rec., 1963, part 9, 149-157. Milhorn, H.T., Benton, R., Ross, R. and Guyton, A.C., A mathematical model of the human respiratory control system Biophys. J. 1965,5,27-46. Longobardo, G.S., Chcrniack, N.S., and Fishman, A.P., CheyneStokes breathing produced by a model of the human respiratory system,]. Appl. Physiol. 1966, 21(6), 1839-1846. Grodins, F.S., Buell, J. and Bart, A J., Mathematical analysis and digital simulation of the respiratory control system, J. Appl. Physiol. 1967, 22, 260-275. Yamamoto, W. and Hori, T., Phasic air movement model of respiratory regulation of carbon dioxide balance, Comput. Biomed. Res. 1971,3,699-717. Milhorn, H.T., Reynolds, W J. and Holloman, G.H., Digital simulation of the ventilatory responses to CO2 inhalation and CSF perfusion, Comput. Biomed. Res. 1972,5,301-314. Damokosh-Giordano, A., Longobardo, G.S. and Cherniack, N.S., The effects of controlled system (plant) dynamics on vent&tory responses to disturbances in CO2 balance, in Regulation and Control in physiologiCa systems, (Eds. Iberall, A.S. and Guyton, AX.) Instrument Sot Amer, Pittsburgh, 1973. Dickinson, C J., A computer model of human respiration, MTP Press Ltd. England, 1977. Priban, I.P., and Fincham, W.F., Self-adaptive control and the respiratory system, Nature, 1965, 208, 339-343. Visser, B.F., Pulmonary diffudon of carbon dioxide Phys. Med. Biol. 1960-61,5 (2), 155-166. Ruth, T.C. and Patton, H.D., (Eds.) in Physiology and Biophysics, W.B. Saunders, 1966, p. 766. Mitchell, R.A., Loeschcke, H.H., Massion, W.H. and Severinghaus, J.W., Respiratory responses mediated through superficial chemosensitive areas on the medulla,j. Appf. Physiol. 1963,18,523-533.
J. Biomed. Eng. 1983, Vol. 5, April
131
Respiratory
16
17
18
19
20 21
22
2.3
24
25
26
27 28
29
30
31
32
model:
W.F. Fincham and F.T. Tehrani
Fir&am, W.F. and Tehrani, F.T., On the regulation of cardiac output and cerebral blood flow. J Biomed. Eng. 1983,5,73-75. Wigertz, O., Dynamics of respiratory and circulatory adaptation to muscular exercise in man, Acta Physiol. Stand. suppl. 1971,363. Tehrani, F.T., Dynamic modelling of the human respiratory system, Ph. D. thesis, 198 1, University of London. Nielsen, M. and Smith, H., Studies on the regulation of respiration in acute hypoxia, Acta Physiol. Stand. 1952, 24, 293513. Otis, A.B ., Fenn, W.O. and Rahn, H., Mechanics of breathing in man, J. Appl. Physiol. 1950,2, 592-607. Lambertsen, C J., Factors in the simulation of respiration by carbon dioxide, in The regulation of human respiration, (Eds. D J.C. Cunningham and B.B. Lloyd) Symp. Proc., 1963, p 257-276, Blackwell. Nielsen, M., Unter suchungen uber die Atem regulation beim Menschen, Skandinav. Arch. Physiol. 1936,74 (suppl. 10) p 87-208. Reynolds, W J., Milhorn, H.T. and Holloman, G.H., Transient vent&tory responses to graded hypercapnia in man, J. Appl. Physiol. 1972,33(l), 47-54. Douglas, C.G., and Haldane, J.S., The causes of periodic or Cheyne-Stokes breathing, J. Physiol. (Lond.) 1909,38,401-419. Astrand, P.O., Human physical fitness with special reference to sex and age, Physiol. Rev. 1956,36, 307-329. Wasserman, K., Van Kessel, A.L. and Burton, G.G., Interaction of physiological mechanisms during exercise, J. Appl. Physiol. 1967,22, 71-85. Comroe, J.H., The hyperpnea of muscular exercise, Physiol. Rev. 1944,24,319-333. Hey, E., Lloyd, B.B., Cunningham, D J.C., Jukes, M.G.M., and Bolton, D.P.G., Effects of various respiratory stimuli on the depth and frequency of breathing in man, Respir. Physiol. 1966, 1, 193-205. Lambertsen, C J., Carbon dioxide and respiration in acid-base homeostatis, Anesthesiology, 1960, 21, 642-651. Brooks, C.M., Kao, F.F. and Lloyd, B.B. (Eds.), Cerebro spinalfluid and the regulation of ventilation, Oxford, Blackwell, 1965. (see Dynamic response characteristics of several C&-reactive components of the respiratory control system, by Lambertson, C J., Gelfand, R. and Kemp, R.A., p. 2 11-245). Cerretelli, P., Sikand, R. and Farhi, LE., Readjustments in cardiac output and gas change during onset of exercise and recovery, J. Appl. Physiol. 1966, 21(4), 1345-1350. Doblar, D.D., Min, B.G., Chapman, R.W., Harback, E.R., Welkowitz, W. and Edelman, N.H., Dynamic characteristics of cerebral blood flow response to sinusoidal hypoxia, J. Appl. Physiol. 1979,46(4), 721-729.
APPENDIX Derivation
(Al) is extended by assuming that half of the ventilatory response to CO3 and ti results from arterial blood via the peripheral receptors and half from the CSF21. Hence ,
A VR = 0.11 e
+ 0.131
,
Pacoz + 0.11 H&F
of the ventilation
(A2)
+ 0.131 PCSFCQ + FACT 3 -K
Using Gray’s relation (1) between arterial ti and levels of hypoxia and hypercapnia.
P a~~z in different
flu’ = 0.65Pacol + 13.5
(9
Similarly, writing the CSF hydrogen ion concentration:
6.h WCSF
=
(BHC%)CSF
(A4)
’ PCSFCOz
Substitution of equations (A3) and (A4) into (A2), and assuming that (BHCO~)CSF’ 6 and X are constant yields : A VR = 0.2025
&Oz
+ 0.2332
Pcs~co
2 (A5)
+FACT3-K
where the parameter values given in appendix 2 have been used. The final step in this development assumes that PCSFCO? can be replaced by the CO, tension in the neighbourhood of the medullary receptors Pccoz (equation 9) and that metabolic activity also has a relatively direct influence on AVR [equation (12)].
APPENDIX
2
Glossary of symbols and values The main symbols denote the variable or parameter and subscripts indicate the model location with reference to a specific chemical, as appropriate. X
denotes the mean value of x
X’
denotes the time delayed version of x
Model variables and symbols Symbol
Definition
Value/units
AVR
Alveolar ventilation
BHCO,
Bicarbonate
ratio
concentration
0.58
[WPWI
1
C
equation
On the basis of the multiple factor theory of Gray’, the alveolar ventilation ratio A VR is: A VR = 0.22 H+ + 0.262 Pacq
+ FACT 3 -K (Al)
where FACT 3 is the response of the arterial chemoreceptors to hypoxia and K is a constant.
132
this paper, equation
J. Biomed. Eng. 1983, Vol. 5, April
In
Volume concentration gas in the blood
of
[l(STPWl
D
Breath duration
(s)
f
Frequency
of breathing
(breath/s)
F
Fractional of gas
composition
H+
Hydrogen ion concentration
I
Inspired gas concentration (fractional volume)
(nmoles/l)
Respiratory
MRB
MRT
model:
W.F. Fincham and F. 1: Tehrani
Metabolic rate in brain
P/s)
basal rates
0,
(0.000925)
Equivalent gas storage space
co2
(0.0009)
ST (tissue)
Metabolic rate in body tissue
(l/s)
basal rates
(0.00352)
0, co2
s
(0.00287)
S, (brain)
(50 1) (1.1 I)
T
Transport delay: arterial blood
(IO s)
(IL
Respiratory system elastance
MRR
Metabolic rate ratio
MRV
Metabolic neural drive factor
P
Air viscosity resistance
P
Partial pressure or gas pressure
Y
Airway turbulent resistance
tpb
Q
- 47)
t . vA
V
(mmw (mmHiid
Blood flow rate
PI4
basal values QT (tissue)
(0.07083)
(2s (brain) RTT
[MR(actual)/ MR(basal)]
Barometric pressure less water vapour pressure at body temperature
H, O/( 1 (BTPS)s)’ ]
Solubility factor (CO2 in CSF)
h
Carbonic acid dissociation constant Time constants:
(795 nmoles/l)
71
Cardiac output controller
(3 s)
6.783 x 10-4) i I(STPD)/ l/mm &I
(l/s)
72
Brain blood flow controller (3 s)
Time
(s)
73
Exercise metabolism dynamic
(30 s)
74
Metabolically derived neural drive
(50 s)
Alveolar ventilation (mean value per breath) basal value VA (rest)
(l/s) (0.0673
Volume
(1)
Depth of central receptor below surface of the medulla
Ventilation controller constant blood gas dissociation constants :
Kl
Oxygen
K2
Oxygen
K3
Carbon dioxide (dynamic)
K4
Central receptor constant
KS
[ 1.285 cm
Metabolic input function to tissue
I/s)
CO, diffusion time constant
Subscr$ ts a
Arterial blood
A
Alveolar gas
b
Barometric condition
B
Brain
c
A point in the medulla
(332
Carbon dioxide
(17.4)
CSF
Cerebrospinal fluid
(0.2) [0.046 (mm Hg)-’ I [0.016 (mm Hg)-’ 3
D
Dead space
E
External ventilation
I
Inspired gas
02
Oxygen
T
Body tissue (excluding brain)
V
Venous blood
VB
Brain venous blood
VT
Tissue venous blood
(15 x 10m3cm)
FACT N Factor N as specified in text K
(8.55 cm Hz O/l) (3.1 cm H, O/l (BTPS)/s)
6
(0.0125)
Model fmumeters d
(1)
(346 x lo3 sec. cme2 .I-’ ) (320 s)
J. Biomed. Eng. 1983, Vol. 5, April
133