The definition and separation of pressure

Pressure is a such a familiar concept that we frequently forget how it is defined scientifically. In thermodynamics pressure P is defined as

where U is the free energy of the system, V is its volume, S is the entropy and N_k are the mole numbers of the different chemical components of the system. This fundamental definition is essentially useless in the clinic because of the impossibility of measuring or controlling the intrinsic and extrinsic parameters.

In mechanics the formal definition of P is

where T = σ_ij is the stress tensor where i,j = [1,2,3] indicate the three cartesian coordinates and we use the summation convention for indices. In the jargon of mechanical analysis, P is the average of the trace of the stress tensor or, more accessibly, the normal component of force per unit area. Definition of the stress tensor is not straightforward in a system involving blood, an extremely complex fluid, and distensible arterial walls. However, this expression is the basis of all of the clinically useful definitions of pressure. Amazingly, despite much effort, no one has been able to show that the thermodynamic and mechanical pressures are equivalent.

Pressure is frequently divided into component parts. Probably the most common division of pressure is the gauge pressure

It is possible to define an absolute pressure P_absolute. However, this is frequently inconvenient because we usually function in a sea of atmospheric pressure. For this reason we generally use pressure to mean the pressure relative to some reference pressure, i.e. a gauge pressure. This is common practice in the catheter lab where the pressure transducer is calibrated to some pressure relative to the heart which includes the atmospheric pressure. I do not know of any clinic that routinely records the atmospheric pressure which means that it is practically impossible to explore the effect of absolute pressure in hemodynamics.

The most famous separation of pressures into different components is undoubtedly the Bernoulli equation

where P is the hydrodynamic pressure, ρ is the density and U is the velocity of the fluid, g is the gravitational constant, H is the height and P₀ is a constant generally called the total pressure. This equation holds along a streamline in steady flow of an inviscid fluid. In unsteady flows it is necessary to add another term involving the velocity potential that accounts for the effects of acceleration and is difficult to evaluate except in the simplest of flows. This is an energy equation (pressure has the units of energy per unit volume) that divides the pressure into potential and kinetic energy. Despite the formal restrictions in its derivation, this equation is very useful clinically and is the basis of estimates of pressure in the cardiac chambers in echocardiographical investigations. Even though blood is not inviscid and the arterial system is highly dynamic, this equation provides a very useful way of interpreting various observations of arterial hemodynamics.

In the context of this meeting, undoubtedly the most common separation of pressure into different components is the separation of the arterial pressure waveform into its forward and backward components shown in Fig. 1. This follows from the work by Westerhof and his colleagues who showed that simultaneous measurements of the pressure and flow waveform could be used for the separation through calculation of the reflection coefficient.¹ A few years later Laximinarayan, working in Westerhof’s group, showed that the separation could be made more conveniently using the characteristic impedance.² I suspect that everyone attending this meeting has made use of this result in their work. In the separation, it is unclear how to apportion the zeroth component (the steady pressure and flow) between the forward and backward waveforms. Westerhof et al. cleverly got around this problem by letting the forward and backward waveforms drift relative to the scales of their measured counterparts. Laximinarayan resolved the problem by not showing any scales at all. This observation, seemingly trivial, is actually important and has a bearing on the definition of the reservoir pressure.

Figure 1

The separation of measured pressure P_M and flow I_M waveforms into their forward P_F and I_F and backward P_B and I_B components. (a) The original results which used the reflection coefficient.¹ (b) The results using the now standard method using the characteristic impedance.²

More than a decade ago, we formulated the reservoir-wave hypothesis that it might be useful to divide the measured arterial pressure into a reservoir pressure and an excess pressure defined as the difference between the measured pressure and the reservoir pressure.³ Our argument was based on the success of the Windkessel model in describing the diastolic pressure waveform and in the original paper we separated the pressure into a ‘Windkessel’ pressure and a ‘wave’ pressure. After publishing that paper we realised that the Windkessel pressure was, by definition, uniform throughout the arterial system and could not describe the observed propagation of this component of the pressure down the aorta. For this reason in subsequent papers we defined the ‘reservoir’ P_r and ‘excess’ P_x pressures.⁴

where τ(x) is the wave travel time from the aortic root to location x. Note that P_r varies in time but is uniform throughout the arteries. This hypothesis was tested extensively in experiments in dogs and there is now a rather extensive literature on the subject.⁵

One of the unexpected observations of our experimental work was that the excess pressure waveform derived by subtracting the reservoir pressure from the measured pressure was very similar in shape to the flow waveform at the aortic root. Assuming that this similarity is always true, it is possible to derive the reservoir pressure from the measured pressure alone.^a

This ‘pressure-only’ algorithm has been used to retrospectively calculate the reservoir and excess pressure waveforms from measured arterial pressure waveforms in a number of clinical studies which have shown that various parameters related to these pressures are significant risk factors for various cardiovascular events.^7–9 These studies indicate that reservoir pressure is, indeed, a useful concept, but that is not the topic of this talk.

Reservoir and excess pressure in the human aorta

One of the primary assumptions in the definition of reservoir pressure is that it is uniform throughout the arterial system and is determined by the global compliance and resistance. The excess pressure, on the other hand, varies at different locations in the arteries and is dominated by local waves. The uniformity of the reservoir pressure is observed in the measurements of pressure and flow at different locations in the dog. It has not, however, been tested in man. Recently we had the opportunity to analyse high fidelity pressure measurements made at 5 different arterial locations in 40 patients. This work has recently been written up and submitted and I will only discuss a few of the relevant results here.¹⁰

The pressure was measured sequentially at 5 locations in the aorta; 1-the ascending aorta, 2-the transverse aortic arch, 3-the descending aorta at the level of the diaphragm, 4-at the level of the renal arteries and 5-at the aortic bifurcation; by withdrawing the pressure transducer tipped catheter to the appropriate location determined by angiography. The pressure was measured for 25 s at each site, Fig. 2. These data were ensemble averaged by an algorithm that excluded outlier beats automatically (based on the magnitude of the correlation of the beat to the ensemble average beat). The ensemble average waveforms P were used to calculate the reservoir P_r and excess P_x pressure using the pressure-only algorithm. Figure 3 shows P − P_d, P_r − P_d and P_x at the 5 measurement sites in one of the patients, where P_d is the diastolic (minimum) pressure measured at each site, The thin black line connects the diastolic points at the different locations and its slope in the x − t plane is an indication of the wave speed in the aorta.

Figure 2

The pressure measured over 25 s in the ascending aorta of one patient (left) was ensemble averaged to give the measured pressure waveform used in the analysis. The standard deviation at each measurement time is not shown but was generally of the order of 4 mmHg, mainly due to the respiratory variation in the pressure seen in the measured data on the left. Some patients exhibited missing or ectopic beats and these were ignored in the calculation of the ensemble average pressure waveform.

Figure 3

The ensemble average pressure P − P_d (solid line), the calculated reservoir pressure P_r − P_d (dashed line) and excess pressure P_x (dash dot line) at the 5 measurement sites in one of the patients. P_d is the diastolic pressure. Since the measurement sites were determined by angiography they are indicated by their index rather than a distance along the aorta.

There are several notable features of this plot, which are common to all of our measurements. First, in all cases, P_r is very close to P during diastole. This is not surprising because the algorithm that calculates P_r derives two of the fitting parameters by fitting an exponential model to the diastolic portion of P. Second, P_x is very close to P during very early systole. This follows from the model of P_r as a solution of the equation for overall mass conservation in the arteries CdPdt=Qin−Qout where Q_in is the inflow from the ventricle at the aortic root and the outflow is determined by simple resistance relationship Qout=(Pr−P∞)R , where P_∞ is a fitted parameter that describes the pressure at which flow through the microcirculation goes to zero and R is the resistance of the microcirculation downstream of the terminal arteries. This outflow during diastole, when Q_in = 0, generates the diastolic decrease in pressure. The outflow will carry on during systole and during the short period at the start of systole when Q_in < Q_out, P_r will continue to fall reaching a minimum when Q_in = Q_out. This decrease and minimum in P_r can be seen during very early systole at each measurement site. This simple physical argument based on conservation of mass in the arterial system explains why P_x is very close to P during early systole. Third, the well-known phenomenon of pressure peaking in the distal aorta can be seen in P − P_d at the different sites. The average pulse pressure P_s − P_d increased from 66.6 mmHg at the aortic root to 77.8 mmHg at the aortic bifurcation. From the figure, it is apparent that this increase in the pulse pressure is accounted for almost entirely by the increase in the maximum of P_x, the maximum of P_r being effectively constant.

These observations were tested by applying intra-class correlation analysis to the data at all of the sites in all of the patients. The results shown in Fig. 4 confirm the conclusions drawn from the results for a single patient shown in Fig. 3. The systolic pressure P_s is shown by the black squares and the black dots indicate the 95% CI in the mean at each site. Mean P_s increases from 135.7 mmHg at the aortic root to 144.3 mmHg at the aortic bifurcation. The peak of P^r decreases very slightly from 51.4 mmHg at the aortic root to 48.7 mmHg at the aortic bifurcation, a difference smaller than the 95% CI for this variable. The peak excess pressure P^x , however, increases significantly from 24.2 mmHg at the aortic root to 40.6 mmHg at the aortic bifurcation. Note that although P = P_r + P_x at every time, the peaks occur at different times in the cardiac cycle and so the sum of the peaks of P_r and P_x do not have to equal the peak of P.

Figure 4

The results of the ICC analysis of patient data at the different sites of measurement: 1-asc aorta, 2-transverse aortic arch 3-desc aorta at the level of the diaphragm, 4-desc aorta at the level of the renal arteries and 5-aortic bifurcation. The solid symbols indicate the mean values and the smaller dots indicate the 95% confidence intervals. (squares, solid line) systolic pressure, (circles, dashed line) peak reservoir pressure and (diamonds, dash dot line) peak excess pressure.

We studied the variation of many other parameters describing the reservoir and excess pressures calculated from the clinical measurements of aortic pressure at different sites in the aorta and anyone who is interested should refer to the paper when it is published. For the purposes of this talk, these measured results support the assumption that P_r is uniform along the aorta. Unfortunately, the flow at the aortic root was not measured and so we could not test the second basic assumption of the pressure-only calculation of P_r that the P_x waveform is proportional to Q_in in man.

To conclude this section of the talk, we note that the definition of reservoir pressure does not rely on periodicity of the aortic pressure. This is in contradistinction to impedance analysis of aortic pressure where periodicity is implicit in the calculation of the Fourier components of the waveform. For this reason, the results of impedance analysis cannot be extended to unsteady phenomena.

Figure 5 shows the instantaneous P, P_r and P_x calculated beat by beat for measurements made in the ascending aorta of one patient exhibiting regular cardiac behaviour. There are variations in the beat by beat measurements of P (solid line) that are most likely due to respiration. P_r (dashed line) and P_x (dash dot line) calculated beat by beat also show some variation from beat to beat. We have also calculated the ensemble average of the beat by beat P_r and P_x and find that they are almost identical to the P_r and P_x calculated from the ensemble average P.

Figure 5

The instantaneous aortic pressure P (solid line), and P_r (dashed line) and P_x (dash dot line) calculated from P beat by beat during the 25 s measurement in the ascending aorta of a patient with a regular cardiac cycle. There are variations from beat to beat, probably due to respiration. These difference may be clearer in P_x where the peak values vary in a very regular way at a frequency that is commensurate with the (unmeasured) respiratory cycle.

Figure 6 shows the pressure measurements in the ascending aorta of a patient who exhibited a number of long beats during the period of measurement. The main message of this plot is that it is possible to calculate P_r and P_x during highly irregular cardiac behaviour. We have not explored the utility of this separation but it seems likely that the excellent fit of P_r to P during extended diastole is telling us something about arterial mechanics during irregular beats and that this mode of analysis my be useful in developing a better understanding of arterial mechanics during these periods.

Figure 6

P, P_r and P_x (format identical to Figure 5) for a patient with a highly irregular cardiac cycle. Note that it is possible to calculate P_r and P_x during both the ‘normal’ and ‘extended’ cardiac cycles. The ability to analyse irregular beats suggests that this mode of analysis may be useful for studying patients who are not exhibiting regular cardiac behaviour.

To summarise this section, we have shown:

1. (1)

   Pressure has been defined and divided in many ways, some useful, some not.

2. (2)

   Reservoir and excess pressure can be calculated from pressure measurements in the human aorta.

3. (3)

   Reservoir pressure propagates down the aorta at the same speed at the measured pressure waveform.

4. (4)

   Reservoir pressure is not non-uniform down the aorta.

5. (5)

   Reservoir and excess pressure can be calculated in irregular beats and may be helpful in understanding these unsteady phenomena.

Asymptotic analysis of the input impedance

For the analysis I assume that the arterial network is a bifurcating tree (not necessarily a symmetrical tree) where all of the properties of the individual arteries are uniform and known: length L, area A, compliance C, wave speed c, wave travel time T and characteristic impedance z. The nodes in the tree are assumed to be bifurcations (internal nodes) where three vessels connect or termini (external nodes) where the terminal arteries are connected to simple known resistances and the root is connected to the ventricle.

We assume that the wave speed c in each vessel is real and given by the Bramwell–Hill relationship c=1ρD where the distensibility D=1AdAdP , the fractional change in area with pressure. For this simple case there are some interesting interrelationships between the variables that are not as well-known might be expected. For a uniform cylindrical tube with constant length the compliance C=dVdP=LdAdP . Substituting from the equation for wave speed C=LAρc2 . We also assume that the vessels are elastic so that the characteristic impedance z=ρcA . Using these results we find that zC=(ρcA)(LAρc)=Lc=T , the time it takes for a wave to travel the length of the artery.

The asymptotic analysis is based on the equations derived by Westerhof relating the impedance Z_0in at the input of vessel n to the reflection coefficient Γ_n at the terminal end of the vessel.¹²

where ω is the frequency of the Fourier component, Z_n is the characteristic impedance and T_n is the wave travel time of vessel n. The reflection coefficient depends on the terminal impedance at the outlet of vessel Z_Tn and the characteristic impedance Z_n

Since trees are parallel in their structure, it is most convenient to recast these equations into terms of admittances Y = 1/Z.

and

To complete the analytical formulation for a bifurcating network, Westerhof showed that the outlet admittance of the parent vessel is equal to the sum of the input admittances of the daughter vessels (Kirchoff’s law).

If we define the maximum wave travel time over all of the vessels in the tree T = max(T_n) and tn=TnT , then we can take ωT ≪ 1 as an asymptotic variable. This low-frequency limit is formally equivalent to the long-wavelength limit. In this limit we can expand the expression for the input admittance neglecting terms of O(ωT)2 .

If we assume that the admittances have the form Y = σ + iωTη where σ and η are real, we can derive the asymptotic recursion relations

where D_n is the set of vessels downstream of vessel n. We assume that terminal artery k is connected to a known terminal resistance R_k. Thus we know for the termini that σk=Sk=1Rk and η_k = 0. By recursion starting at the terminal arteries back to the root of the tree (n = 1), we find after a lot of algebra where and where K is the set of all termini and N is the set of all of the vessels in the network. The first result says that in the low-frequency limit the real part of the input admittance of the arterial tree is equal to the sum of all of the terminal admittances. This is the net admittance of the arterial tree which determines the steady (ω = 0) behaviour of the network and so this result is expected. The second result is less expected. It is usually assumed that the first order behaviour of the system will depend on the sum of the compliances of the vessels in the network whereas this analysis indicates that it is the sum of the individual vessel compliances weighted by a factor that depends on the square of the product of the characteristic impedance of the vessel and the net admittance of all of the termini downstream of the vessel.

I have not been able to think of a simple physical explanation for the form of this effective, weighted compliance, but there is a simple argument that indicates the need for some form of weighting. Consider a terminal vessel whose characteristic impedance matches its terminal resistance. The reflection coefficient at the terminal end of this vessel will be zero (S_nz_n = 1) and so the wave in this vessel will not be reflected back into the rest of the tree. This means that the rest of the network is not influenced by the compliance of that vessel and we would not expect it to affect the first order behaviour of the tree. The weighting (1−S¯n2zn2) suggests that the influence of the compliance of vessel n is influenced by the degree to which the wave in that vessel is reflected back into the network.

Reverting to impedance, these results indicate that the low-frequency behaviour of the tree is

This is familiar as the input impedance of a 2-element Windkessel with net arterial resistance R and compliance equal to the effective compliance C¯1* defined above. The grey dashed lines in Fig. 7 are fits of this relationship to the low-frequency input impedance calculate by Taylor. Unfortunately he did not give enough information about the lengths of the arteries in each realisation of his random models to enable us to test the accuracy of the estimate of effective compliance. In identical calculations that I have performed where the distribution of lengths is known, the use of the effective compliance provided an excellent fit to the low-frequency amplitude of the calculated input impedance.

Low-frequency flow and pressure

Having derived the low-frequency input impedance of the arterial tree, we can now calculate the pressure and flow throughout the network starting at the root and working downward through the tree. We first observe that flow in the parent vessel (index a) at a bifurcation distributes into the daughter vessels (indices b & c) depending on their input impedances.

where (•)˜ indicates the Fourier transform of (•. The last equation follows from Kirchoff’s law at the bifurcation. In a tree there is a unique path from the root of the tree to vessel n which can be written as a list of edges A_n = [e₁, e₂,e₃,…,e_j]. Using the above relationship, the division of flow along this path is

Substituting the low-frequency admittances

where we have defined the effective time constant Ti*=C¯i*S¯i . The last equation comes from the asymptotic approximation to the product neglecting terms of order O(ωT)2 .

By definition of admittance, Q˜0n=Y0nP˜0n , which allows us to express the last equation in terms of the Fourier pressure component

The summation in this expression is simply the summation of the effective wave travel times over all of the vessels on the path from the root of the arterial system to vessel n. Using the bar notation for net quantities previously used for S¯n and C¯n* , we define this sum as T¯n* . We now apply the inverse Fourier transform and the shift theorem to conclude that

in the low-frequency limit. In words, this result indicates that the low-frequency pressure waveform at any point in the arterial tree is equal to the low-frequency pressure waveform at the aortic root delayed by the effective wave travel time from the root to that point. The dependence on the effective wave travel time (based on the effective compliances) instead of the wave travel time is another unexpected result of this analysis.

We conclude this section by noting that the low-frequency component of the pressure shares the fundamental properties of the reservoir pressure, it is uniform throughout the arterial system but delayed by a wave travel time. In the present company I am reticent to claim that it is the reservoir pressure; this will require further analysis and study. However it is clear, whatever your biases about impedance and wave intensity analysis, that the low-frequency component of the pressure has some interesting properties and certainly warrants further study.

To summarise this section, we have shown:

1. (1)

   The low-frequency input impedance is Windkessel-like: Z01=R1+iωRC¯1* .

2. (2)

   The relevant compliance is an effective compliance: C¯1*=∑n∈NCn(1−S¯n2zn2)) .

3. (3)

   The low-frequency pressure is uniform throughout the arterial system, but delayed by a wave travel time: P0n(t−T¯n*)=P01(t) .

4. (4)

   The wave travel time is the effective wave travel time: Tn*=Tn(1−S¯n2zn2)S¯nzn .

The generation and evolution of wave fronts in the arterial tree

The 1-D conservation equations in each vessel together with a suitable tube law describing the change in area of the vessel with pressure give a hyperbolic system of ODEs. They are amenable to solution by the method of characteristics introduced by Riemann.¹⁸ The salient feature of this solution for our model which neglects viscous effects is that any perturbation of conditions in the vessel will propagate as a wave in both the forward and backward direction. Following the example of gas dynamics where the 1-D model is well-developed, we will assume that the basic form of wave is a wave front, a step change, either positive or negative, in the pressure and the flow across the wave. It is easy to build any waveform by a sequence of discrete wave fronts which makes them convenient as a basis for any waveform. We make the further simplification that the wave fronts are linear in the sense that they are additive and do not interact with each other when they cross. This assumption is essentially the acoustic limit in fluid dynamics and can be ensured in the absence of shock waves simply by considering wave fronts whose amplitudes are sufficiently small.

With all of the simplifying assumptions we can think of each wave front as a ‘particle’ that propagate unchanged through the edge (vessel) until it encounter a node. When it encounters an external node the wave front is reflected with a magnitude that is given by a reflection coefficient. When it encounters an internal node the wave front produces three wave fronts, a reflected wave front with a magnitude given by a reflection coefficient and two transmitted waves with magnitudes given by a transmission coefficient. The problem thus reduces to following the evolution of the single wave front introduced at the root at t = 0.

Because the interaction of a wave front with an internal node produces three waves, the total number of waves in the system will grow exponentially. The amplitude of the new waves decreases on average and we will see that the amplitude gets exponentially small as the number of waves gets exponentially large. The net effect of the myriad numbers of vanishing small waves is not immediately apparent and we have to be careful in the analysis not to truncate results inappropriately.

To facilitate the tracking of the waves in the system we introduce the idea of a wave history W_n = {e₀,e₁,e₂, …,e_i,e_i+1,…,e_N−1,e_N} which is the sequence of edges visited by wave n. For waves introduced at the root, e₀ = 1. N is the number of edges visited by the wave which leads to the concept of a step, which is defined as going from one node of the edge to the other. Because each edge has a wave travel time associated with it, the time it takes a wave to execute one step depends upon the vessel it is in. This leads to difficulties in calculating the location of a wave at a particular time whereas locating a wave after a given number of steps is simple. The synchronisation of steps with time is a final step in the analysis.

After N steps the total time (duration) of wave T_n is the sum of all of the wave travel times in all of the edges it has traversed

The amplitude of the wave A_n is the product of all of the reflection and transmission coefficients at all of the nodes it has encountered

where Γ_eiei+1 is the reflection/transmission coefficient as the wave passes from edge e_i to edge e_i+1 in its wave history.

The exponential growth of the number of wave fronts and hence wave histories is an important feature of this theory. The exact number depends on the topology of the network. For large N the number of waves is given by the recursion formula Nn+1waves=3NiNnwaves , where N_i is the number of internal nodes. We also know that mean distance moved during a random walk on a tree, like random walks on a line or a square grid, is proportional to the square root of the number of steps. Thus, for a tree with altitude A the average number of steps needed to traverse the whole tree is A². For the 55 artery model, A = 14 and so approximately 200 steps are needed for the waves to spread throughout the tree. After 200 steps there are approximately 1.4 × 10⁸¹ waves.^b The vast number of waves indicates the necessity and the validity of a statistical approach to the problem.

Fortunately the importance of large networks on modern life has stimulated many advances in the field of network theory and we will adopt some of these methods into our analysis. Networks can be described in many ways but for our purposes the most convenient representation is the connectivity matrix. For a network with N nodes and E edges the connectivity matrix is an N × E matrix with each row representing one node and each column representing one vessel (edge). By definition each edge is connected to 2 nodes, one identified as the input node and the other the output node. In rooted binary trees, it is straightforward to describe the node closest to the root as the input node. In more general networks with loops, this is not possible and the designation of input and output nodes can be arbitrary. Each column in the connectivity matrix has −1 in the row of the input node and +1 in the row of the output node. All other entries in the row are zero. This matrix provides a complete description of the topology of the network. In particular summing the number of non-zero entries along a row gives the degree of the node defined as the number of edges that are connected at that node. We assume that the arterial network can be modelled by nodes of degree 3 (internal bifurcations) or of degree 1 (external nodes). Real arterial networks sometimes include trifurcations but these are rare and it is almost always possible to model them as two bifurcations separated by a very short edge.

Given the topology of the network, it is necessary to find a compact way to represent the reflection and transmission coefficients at the nodes. For an external node k connecting a resistance R_k to an edge with a characteristic impedance z_k, the reflection coefficient Γ_k is

This is the same as the reflection coefficient defined in the impedance analysis in Section Reservoir and excess pressure along the human aorta. The reflection and transmission coefficients for an internal node depend upon the edge in which the wave front approaches the node. For an internal node i connecting edges a, b and c the reflection and transmission coefficients are given by the 3 × 3 matrix

where the reflection coefficient for a wave in a can be written most conveniently in terms of the characteristic admittances y=1z γ_b and γ_c can be found by permutation of the indices. The transmission coefficient for a wave with a reflection coefficient γ is 1 + γ. Note that − 1 ≤ γ ≤ 1 and that 0 ≤ 1 + γ ≤ 2. It can also be shown with a little algebra that the sum of the reflection coefficients over all of the edges (γ_a + γ_b + γ_c) = − 1. With even more algebra it can be shown that Γi2=I where I is the identity matrix. This means that the matrix of reflection/transmission coefficients is its own inverse or, in other words, it is idempotent of degree 2. Physically this means that with all of the simplifying assumptions that we have made, particularly the neglect of viscous effects and the linear behaviour of the wave fronts, no energy is lost when a wave front interacts with an internal node. This is not true for interactions with external nodes where the amplitude of the reflected wave is less than one which means that some energy is lost.

Knowing the reflection and transmission coefficients at each node, we can now construct a global matrix of coefficients for the whole network. For a matrix with E edges, this will be a 2E × 2E matrix where each edge has a separate row for the forward and backward waves in the edge. For a large network this can be a rather large matrix, for example the matrix for the 55-artery model is 110 × 110, but the matrices are sparse (mainly zero elements) and they can be easily accommodated by current computers. We term this matrix the ‘scattering matrix’ by analogy to scattering theory in acoustics and microwave transmission theory. Calculating the amplitude of the waves in the system now reduces to a problem of matrix multiplication. Given an initial wave at the root with amplitude a₀ =[1,0,0,0,…,0]^T, i.e. a unit step wave at the root) the amplitude of the waves in the system after N steps is given by

Thus, the history of the waves in the system after N steps is reduced to evaluating the powers of Γ.

Eigen-modes and time constants

Although computers can calculate the power of a large matrix very quickly, the representation and interpretation of the results is difficult. Fortunately linear algebra provides a way to deal with the problem that provides results that are easy to interpret; eigenvalues and eigenvectors. An eigenvalue λ and eigenvector v of matrix Γ are defined by the equation

That is, multiplying an eigenvector by the matrix results in a scale factor times the eigenvector. This is very convenient when we are taking powers of the matrix because successive multiplications by the matrix yields

Thus, the power of the matrix times an eigenvector is just the power of the scalar eigenvalue times the eigenvector, and the problem of matrix multiplication further reduces to multiplication by a scalar.

Because they are central to matrix calculations, the theory of eigenvalues and eigenvectors is a well-studied branch of linear algebra. In general an N × N non-degenerate matrix can have N eigenvalues each of them with an associated eigenvector. The eigenvalues determined by finding the roots of an N th order polynomial, however, are not necessarily unique and it is difficult to determine the multiplicity of eigenvalues without calculating them. Also, for a real matrix the eigenvalues can be complex in which case they occur as conjugate pairs, and these have complex conjugate pairs of eigenvectors associated with them.

Traditionally the eigenvectors are taken to be unit vectors because it is the direction of the eigenvector, not its amplitude, that is important. Since the sum of the square of the absolute value of the coefficients of a unit vector is equal to 1, the eigenvectors can be thought of as eigen-modes associated with each eigenvalue. The evolution of these eigen-modes as the number of steps in increased is determined by the absolute value of the eigenvalue. If we define A_m(N) to be the amplitude of mode m, after N steps then

where we have written the power in its exponential form. If |λ| > 1, log | λ | > 0 and the amplitude of the mode will grow exponentially with the number of steps. If | λ | < 1, log |λ | < 0 and the amplitude of the mode will decrease exponentially. The rate of decrease being determined by the value of log | λ |.

Physically we are interested in the behaviour of the system with time, not with the number of steps taken by the wave fronts. To relate the results in terms of steps to time forces us to turn to statistical methods. Given the wide range of wave travel times in the arterial network (in the 55 artery model they range from 1.9 ms for the celiac artery to 71.5 ms for the two brachial arteries). Fortunately, we can apply statistical arguments with a high level of confidence, given the large number of waves that are generated.

For our analysis, we assume that the time it takes to execute N steps is given in a probabilistic sense by t = N〈T〉, where 〈T〉 is the appropriate average wave travel time. There are a number of candidates for the appropriate average, but we believe that best choice is the weighted average wave travel time for each particular mode. As discussed previously, this weighting is the square of the magnitude of the coefficients of the eigenvector which defines the mode. That is, for mode m

where T_i is the wave travel time of vessel i and v_mi is the i th component of the eigenvector v_m associated with mode m. Substituting for N in the equation for the amplitude of a mode where we have defined the time constant for the mode

The minus sign in the definition is usual since we are only interested in eigenvalues with |λ| ≤ 1 with negative logarithms.

The conclusion of this analysis is that the for large times, when we can expect t = N〈T〉 to be a good estimate, each of the modes of the system will contribute an exponentially decaying term with different time constants for each of the modes. The long-time behaviour of the system will depend on the mode with the longest time constants, since the other modes will die away more quickly. Ideally this time constant would correspond to the time constant that is usually cited for the diastolic pressure decay. For the limited number of cases which I have studied, this is not the case and the largest time constants are commonly two or three times larger than the expected diastolic time constant. This could be due to deficiencies in the models I have used, an extraordinary number of parameters is needed to define a realistic arterial network. Alternatively, it could be related to problems associated with fitting a single exponential to a multi-exponential function.

Results for the 55 artery model

The above analysis is abstruse and undoubtedly difficult to follow without a background in large dynamical systems. I will try to make the important points clearer by presenting the results of this analysis applied to the 55 artery model of the human circulation. This model has been used extensively in model calculations and I have based my work on a slightly simplified version of the model used by Alastruey and his colleagues in their numerical computations (refer to Ref.¹³ for an example of their work). The main simplifications are the assumption of uniform arteries (no taper) and real terminal impedances (resistances with no compliance). With these assumptions it is possible to calculate the characteristic impedance of each vessel and, knowing the connectivity of the network, the reflection and transmission coefficients at all of the nodes. These were used to construct the 110 × 110 scattering matrix Γ.

The eigenvalues and eigenvectors of Γ were calculated using MatLab. The level of multiplicity of the eigenvalues was fairly large and there were only 30 unique absolute values of the eigenvalues. The eigenvectors associated with each of the unique eigenvalues defined 30 different modes for the network. Table 1 shows the time constant τ in s, the magnitude of the eigenvalue |λ | and the mean wave travel time 〈T 〉 for the modes with the 6 largest time constants and for the mode with the smallest time constant, ranked according to the time constant.

Representations of the modes is difficult because the indices of the arteries in the model are somewhat arbitrary (see Fig. 8). We show them graphically in Fig. 9 where the grey level of each artery is an indication of its weight in the mode. The mode with the largest time constant (upper left) is dominated by the arteries in the arm. The magnitude of the eigenvalue for this mode is only the 4th largest but the largest wave travel time is for the brachial artery and this is reflected in the time constant for the mode. The modes with the next 3 largest time constants are dominated by the legs and lower arms. The mode with the 5th largest time constant corresponds to the mode with the most uniform distribution of arteries. The smallest time constant (not shown) is interesting because it is almost completely dominated by the celiac artery which has the smallest wave travel time in the model.

Figure 9

The eigen-modes corresponding to the 6 largest time constants. The grey scale of each artery is an indication of its contribution to the mode.

Implications for clinical measurements of arterial pressure

These results for the evolution of a wave front introduced at the root of the arterial network show that the input wave front will excite a large number of modes in the network which decay with different time constants. Calculating these modes and their time constants requires full knowledge about the connectivity and physical properties of all of the arteries in the network – god-like knowledge that is not available to the clinician. We have to conclude, therefore, that the theory does not offer direct information about the measured waveform. It does, however, offer insights into the physical processes involved which can increase our understanding of the interaction of the heart and the arteries.

The theory indicates that there are many modes that will be stimulated by waves generated by the heart and, after sufficient time, they will decay exponentially with different time constants. Fitting measured data to a multi-exponential function is “notoriously difficult”.^20,21 This means that it is very difficult to determine the multiple time constants from measurements of the diastolic pressure decay because there are many combinations of coefficients and time constants that can fit the same data over the limited span of diastole.

This difficulty is illustrated in Fig. 10 where we have fit the decaying pressure during a single diastolic period to (dashed line) a single exponential model with a zero asymptote, (dash dot line) a single exponential model with a free asymptote and (dotted line) a double exponential model with a zero asymptote. The fit of the single exponential without an asymptote (R² = 0.9663) is significantly worse than the fit of the other models which have identical R² values (R² = 0.9995) but the time constant (τ = 0.68s) is closer to the textbook value of the diastolic time constant. The single exponential model with a free asymptote gives a significantly shorter time constant (τ = 0.28s) but a surprisingly high value of the asymptote. The asymptote of the double exponential model is zero but the largest time constant (τ = 5.58s) means that this asymptote is reached very slowly. The smaller time constant for the double exponential model (τ = 0.24s) is close to the time constant for the single exponential with an asymptote model. Fitting higher order exponential models to the data do not lead to any significant increase in the fitting statistics and the results are unpredictable and generally uninterpretable.

Figure 10

Comparison of different exponential models to the diastolic pressure measured in the ascending aorta. (solid line) P measured for 3 beats in the ascending aorta using an invasive pressure catheter calibrated to atmospheric pressure. (dashed line) data fitted by a single exponential with zero asymptote P = P₀e^−t/τ₁. (dot dash line) data fitted by a single exponential with a free asymptote P =(P₀−P_∞)e^−t/τ₂+P_∞. (dotted line) data fitted by a double exponential with zero asymptote P=P₀(αe^−t/τ₃+(1−α)e^−t/τ₄). The single exponent with zero asymptote does not fit the data during diastole as well as the other models. The fits using the single exponential model with asymptote and the double exponential model with zero asymptote are statistically indistinguishable (R²=0.9995 for both cases) but diverge significantly for times longer than the diastolic period. The fitted time constants are τ₁=0.68s, τ₂=0.28s, τ₃=0.24s and τ₄=5.58s.

The only message we can draw from this example is that there are many exponential models that will fit the decaying diastolic pressure. The fitted time constants depend on the choice of model and so it is probably wrong to cite the time constant as a unique parameter describing diastole. The second message is that using any of the models fitted during the diastolic period can be used with confidence during diastole but should be used with caution when extrapolating to longer times.

The implications of this to the determination of the reservoir pressure which involves fitting the diastolic pressure with a single exponential with a free asymptote and extrapolating this into the systolic period are uncertain and warrant further research. The work described here suggests that the reservoir pressure should have a multi-exponential form but the uncertainty in fitting such a model to diastole limits the practical usefulness of this result. It is clear, however, that the long-time effect of the wave fronts generated by the input from the heart to the arteries is not negligible. This can be illustrated very clearly by looking at some computational results for the simplified 55 artery model (kindly generated by J. Alastruey) where the pressure in the ascending aorta was calculated for a given periodic flow from the heart. The calculations were started at arbitrary initial conditions and continued until the calculated pressure became periodic, at which point the input from the heart was stopped, corresponding to an indefinitely extended diastole. The result of that calculation is shown in Fig. 11.

Figure 11

The pressure in the ascending aorta calculated by J. Alastruey for the simplified 55-artery model used in this study. The black line shows the periodic pressure that is generated by a given periodic inflow from the heart. Note that the relatively large fluctuations in pressure seen in late systole and mid-diastole and persist for at least 2 s into the extended diastole do not appear in his calculations which include a significant level of compliance the terminal impedances in the model. The grey lines are the results obtained when the input from the heart was stopped to model the effect of an extended diastole. Because the pressure was periodic, with a period of 1 s, the light grey lines show the pressure that would result if inflow from the heart was stopped after 1, 2 and 3 beats. Looking at the fourth beat in the sequence, we see that nearly half of the diastolic pressure (the minimum pressure at the start of the next systole) is due to the pressure that would be present in the artery if the inflow was stopped at the end of the second beat. Furthermore, about half of that pressure is due to the pressure remaining in the arterial system due to the first beat. This exponentially decreasing contribution from previous beats further and further away from the current beat is typical of exponential processes and is clearly not negligible. It is this contribution to the arterial pressure that is modelled by the reservoir pressure.

The light grey lines indicate the pressure in the ascending aorta calculated assuming that the flow from the heart into the arteries was stopped before the start of the next systole. Because the calculated pressures are periodic, we can plot the result of extended diastole occurring after first, second or third beat. In this way we can see the effect of the previous systoles on subsequent beats. In this example, approximately half of the diastolic pressure at the start of a beat is due to pressure that was generated by the previous systole. The pressure due to the systole two beats before the current beat contributes about a quarter of the diastolic pressure and this contribution falls of exponentially with all of the previous beats. For this calculation we know that the pressure will asymptote to zero and fitting a single exponential model with a zero asymptote give a time constant to the extended diastole τ = 1.02s. Looking at Table 1 we see that this fitted time constant is in the middle of the 6 largest time constants calculated for the this model. This suggests that the calculated data could be fitted equally well with a single exponential model or a high order exponential model with time constants given by the theory.