Euler, Father of Hemodynamics

Sylvio R. Bistafa

University of Sao Paulo, Sao Paulo, Brazil.

DOI: 10.4236/ahs.2018.72007                

Abstract

In 1742, Euler presented an essay on the “Principles for determining the motion of the blood through arteries”. This is the first known work on the mechanics of flows in elastic tubes, in which Euler applied his equations to analyze the flow of blood through arteries, driven by a piston pump simulating the heart. However, Euler did not recognize the wave nature of his equations, which led him to a dead end on trying to find a closed form solution. Nonetheless, it will be shown that the hemodynamic equations developed by Euler about 275 years ago, still undergird the most advanced numerical methods in use today for blood flow analysis in arterial networks. Therefore, Euler’s pioneering and seminal work in the area of blood flow justifies he be called the father of hemodynamics.

Keywords

Hemodynamics, Arterial Wave Mechanics, Cardiovascular Mechanics, Modeling the Features of the Arterial Network, Application of the Method of Characteristics

1. Introduction

In 1742, the Dijon Academy launched its first contest on a subject with the title “To determine the difference in velocities between a liquid that flows through elastic and rigid tubes” (Déterminer la différence des vitesses d’un liquide qui passe par des tuyaux inflexibles et de celui qui passe par des tuyaux élastiques), in which Leonhard Euler (1707-1783) submitted a manuscript (presumably) with the title Principia pro motu sanguinis per arterias determinando. Truesdell (1955) traces the mishaps of this publication, indicating that: “[∙∙∙] While it was the academy’s principle to retain all manuscripts, the box on whose label 1742 appears contains no memoir earlier than 1765”. Principia pro motu sanguinis per arterias determinando appeared in print much later, in 1862, and in fact dated 1775.

As originally published in Euler’s Opera Postuma (Euler, 1862) , E 855 Principia pro motu sanguinis per arterias determinando of 1775 is an incomplete manuscript, lacking paragraphs 1 to 14, whereas the complete version appeared in Euler’s Opera Omnia (Euler, 1979) (its Preface contains details about the missing parts). In fact, although in the introduction to the problem (paragraphs 1-8), Euler focuses on the discussion on blood flow through arteries, to such a degree to even proposing adhoc models for the behavior of their elastic cross-sections, surprisingly, and also somehow disappointing, is that the main body of the manuscript (paragraphs 9-34) is devoted to the modeling of flows through rigid tubes, driven by a piston pump. It is only in the remaining sections of the manuscript (paragraphs 35-43) that Euler investigates the formulas for the motion of fluids in elastic tubes, ending the work with a very pessimistic statement that “[∙∙∙] since there appears to be no way in which this can be completely resolved and the investigation can be considered to transcend human powers, the work will end here”.

As Truesdell (1955) pointed out “[∙∙∙] It is interesting that Euler’s first paper on hydraulics should concern this extremely difficult subject [∙∙∙] since (E 855) employs concepts Euler was not to develop until 1755”. In fact, what is considered to be Euler’s first hydraulic paper appeared in printing only in 1754 (Euler, 1754) , where Euler dealt with the problem of raising water to an elevated reservoir with a piston pump. As Cerny & Walawender (1974) noted “[...] This work (on blood flow) may be the first mathematical treatment of circulatory physiology and hemodynamics”. Euler perhaps could very well be called the “Father of Haemodynamics”. In this publication, these authors essentially translated the incomplete form of E 855, not on equal footing with the original material, to omit most of the verbose parts, and by adding comments and interpretations from a more modern point of view. Nonetheless, all the mathematical developments, assumptions and justifications were included. In this manner, the authors succeeded in providing a very objective translation, covering the essential parts of Euler’s manuscript on blood flow. A complete English translation of E 855 by the author¹ differs from the one above, not only because it is based on the complete form of manuscript, but also because, it had been accomplished by a pari passu type of translation, which reflects with more fidelity the original Latin form of the manuscript in English.

It is well known that the modern understanding of the cardiovascular system undoubtedly starts with the work of William Harvey (1578-1657) who published his discovery of the circulation of blood in 1628. A historical review on the subject (Parker, 2009) indicates that before Euler, other investigators had already been concerned with blood circulation. Giovanni Borelli (1608-1679), which is seen by many as the father of bioengineering, studied the contraction of the heart and its interaction with the arteries, where he understood the capacitive effect of the elastic arteries on smoothing the flow of blood (now known as the Windkessel effect to be later discussed). His works were posthumously published in 1680 in De Motu Animalium. Stephan Hales (1677-1746) published in 1733 a series of papers that he had presented before the Royal Society as Statical Essays: Containing Haemastaticks. It is full of original observations about the mechanics of the cardiovascular system including the first measurements of in vivo blood pressure.

As we shall see, the one-dimensional modeling of the human arterial system introduced by Euler in 1775, results in a system of hyperbolic partial differential equations expressing the conservation of mass and momentum for inviscid flow. In order to close the problem, he also suggested two possible, but unrealistic, constitutive equations describing the behavior of the elastic wall with changes in the internal pressure. Apparently, Euler did not recognize the wave-like nature of the flow and was not able to find a solution for his equations, citing “insurmountable difficulties”. In fact, as we shall see, solutions to Euler’s hemodynamic equations in the arterial system were only possible rather recently, with the advances of digital computer simulations in the late 20th century. Although Euler’s equations had remained essentially the same, model refinements recently introduced, such as the inclusion of wall friction, more realistic non one-dimensional velocity profiles and pressure-area relationships, together with newer numerical schemes, have improved the ability to capture the main features of pressure, flow and area waveforms in arteries.

Being a rather multidisciplinary and specialized field of study, and because of the vast amount of work that has been done in the area of blood flow analysis with all its involvements and complexity, the present work does not intend to trace the evolution of the human arterial system analysis, but to show the pervasiveness of Euler’s model and his hemodynamic equations as a basis for the most advanced numerical models in use today. As noted by Alastruey, Parker, & Sherwin (2012) , despite the complexity of the vascular system, the one-dimensional (1-D) “reduced” Euler’s model is commonly applied to simulate the changes in blood flow and cross-sectionally averaged blood pressure and velocity in time and only along the axial direction of larger arteries, with reasonable accuracy and computational cost. Interesting to note is that Euler’s hemodynamic equations describing flow in elastic arteries have the same mathematical structure as the compressible gas dynamics and shallow water equations; the elasticity of the artery vessel wall taking the place of the compressibility of the gas in the former, and the channel geometry in the latter equations.

2. Euler’s Hemodynamic Equations

In E 855, Euler applied the principles of mass conservation and momentum conservation to the one dimensional flow of an incompressible fluid through an elastic tube driven by a piston pump to obtain, in Euler’s notation, the equation

$\left(\frac{\partial s}{\partial t}\right)+\frac{\partial \left(vs\right)}{\partial z}=0$ (1)

from the mass conservation principle, and

$2g\left(\frac{\partial p}{\partial z}\right)+v\left(\frac{\partial v}{\partial z}\right)+\left(\frac{\partial v}{\partial t}\right)=0$ (2)

from the momentum conservation principle, which it is recognized as the one-dimensional form of what is now known as the Euler equations of fluid dynamics, which are used for inviscid flows calculations.

In these equations $s=s\left(z,t\right)$ is the cross-sectional area of the artery, $v=v\left(z,t\right)$ is the velocity, and $p=p\left(z,t\right)$ is the pressure (in fact, the pressure head), t is time, and z is the coordinate along the tube. Equation (2) has units of acceleration (force per unit of mass), hence p has units of length, and 2 g is actually the gravity, due to Euler’s peculiar form of writing equations involving forces.

To solve for s, v and p, Euler had to add an algebraic pressure-area relationship, for which he proposes two expressions

$s=\frac{\Sigma p}{c+p}$ , or $s=\Sigma \left(1-{\text{e}}^{-\frac{p}{c}}\right)$ (3a, b)

where $\Sigma$ is the cross-sectional area of the tube at infinite pressure, and c is a constant quantity that depends on the degree of elasticity of the tube. By adopting Equation (3b), it would result in a simpler expression for the pressure gradient term in Equation (2). Nonetheless, he makes no further use of either in the analysis.

Since Euler had succeeded in finding a solution for p and v for the rigid tube case, he then tried to do the same for the elastic tube case, but, with no avail. Nonetheless, as far as the solution for rigid tubes is concerned, his approach is impeccable. By invoking the continuity equation and the momentum equation, and with great analytical skills, he was able to reduce the problem to the integration of a single differential equation, for the flow of liquids through rigid tubes with variable cross sections. From this, he gave closed form expressions for the velocity and pressure for tubes of uniform cross sections. Euler will then return to this matter again in the memoir of 1754 (Euler, 1754) , where he considers the problem of raising water to an elevated reservoir with a piston pump, where he derives an expression for the pressure in any location along the pipe, which, with the exception of the gravity terms and different notation, is exactly in the same form of his general equation developed in E 855.

The Wave Nature Solution of Euler’s Hemodynamic Equations―The Linear Solution

An approximate solution to Euler’s hemodynamic equations can be obtained from the linear theory, which is valid for small amplitude disturbances. Let $s=s_0+a$ , such that $\left|a\right|\ll s_0$ , $\left|v\right|\ll 1$ , and $p\equiv p_e-p_0=P\left(s\right)$ , the transmural

pressure. Then, $\frac{\partial p}{\partial t}=P^{\prime }\left(s_0\right)\frac{\partial a}{\partial t}+O\left(a^2\right)$ .

Under the assumption that the artery vessel has uniform cross-section (s is a constant along z), and the blood has density ρ, then Equations (1) and (2) give

$\frac{1}{s_0{P}^{\prime }\left(s_0\right)}\frac{\partial p}{\partial t}+\frac{\partial v}{\partial z}=0$ (4)

$\frac{\partial v}{\partial t}+\frac{1}{\rho }\frac{\partial p}{\partial t}=0$ (5)

on neglecting small quantities. Here p is actually the pressure in pascal.

Eliminating v gives

$\frac{{\partial }^{2}p}{\partial t^2}=c_0^2\frac{{\partial }^{2}p}{\partial z^2}$ (6)

where $c_0^2={\left(\rho D_0\right)}^{-1}$ , in which $D_0$ is the distensibility of the artery vessel, defined as $\frac{\text{d}}{\text{d}p}\left(\frac{a}{s_0}\right)=D_0^{-1}$ , or equivalently $D_0=\frac{1}{s_0{P}^{\prime }\left(s_0\right)}$ .

Equation (6) is recognized as the wave equation, which admits the so-called D’Alembert solution² in the forms

$\begin{array}{l}p=f\left(z\pm c_{0}t\right),\\ v={\left(\rho c_0\right)}^{-1}f\left(z\pm c_{0}t\right),\\ a={\left(\rho c_0^2\right)}^{-1}f\left(z\pm c_{0}t\right),\end{array}$ (7)

where $f\left(z\pm c_{0}t\right)$ is the wave form, which propagates with speed $c_0$ ; the − sign indicates waves travelling in the positive z-direction, and the + sign indicates waves travelling in the negative z-direction. This linearization is valid provided $\left|v\right|\ll c_0$ , with $c_0$ given by the Moens-Korteweg formula, to be presented next.

After Euler’s first attempt to apply the principles of hydraulics to the flow of blood through arteries, Thomas Young (1808) derived the wave speed using an analogy between waves in a compressible fluid in a rigid tube and waves in a incompressible fluid in an elastic tube, where the speed of propagation follows from Newton’s theory of speed of sound in air. For 70 years, Young’s rather confusing derivation of the wave speed remained forgotten, until Moens in 1877/1878 and Korteweg in 1878 found independently the correct formula for the speed of propagation of waves given by

$c_0=\sqrt{\frac{E}{\rho }\frac{h_0}{D_0}}$

where E is the modulus of elasticity of the material, $h_0$ and $D_0$ are, respectively, the wall thickness and the inside diameter of the artery vessel, and ρ is the density of the blood. This formula for the wave speed became known as the Moens-Korteweg formula. Tijsseling & Anderson (2012) gives a detailed historical account on the development of this formula by both authors.

3. Refinements to Euler’s Hemodynamic Model

As we saw earlier, Euler’s pressure-area relationship is considered a vulnerable point of Euler’s hemodynamic model, and, therefore, more realistic constitutive wall-laws have been introduced. Also viscous flow effects have been added to his equations, as well as non one-dimensional velocity profiles have been considered.

3.1. Pressure-Area Relationship

The adhoc pressure-area relationships given by Euler (Equations (3a) and (3b)), were not derived from any theory of elasticity. In fact, arterial walls exhibit elastic and viscous behavior (viscoelastic behavior), and as the pressure wave propagates inside the tube, it distends (elastic effect of the walls), with the wave being progressively attenuated (dissipation effect due to the viscosity of the walls), with its amplitude decreasing exponentially during propagation. The attenuation is primarily caused by the viscosity of the walls. Furthermore, the stress-strain relation for flexible tubes such as arteries is nonlinear (Saito et al., 2011) . However, as a first approximation, we can assume that the artery is a linearly elastic material that obeys Hooke’s law, where the stress is related to the strain in both the longitudinal and circumferential direction through the Young’s modulus E.

The artery viscoelastic behavior can be modeled with different mathematical forms. One of such formulations leading to a pressure-area relationship considers a single linear term and an attenuation term according to the Kelvin-Voigt model. This model, also called the Voigt model, can be represented by a purely viscous damper and a purely elastic spring connected in parallel. Since the two components of the model are arranged in parallel, the strains in each component are identical, and then ${ϵ}_s={ϵ}_d$ , where the subscript s indicates stress-strain in the spring, and the subscript d indicates the stress-strain in the damper. The total stress in the system will be the sum of the stress in each component ${\sigma }_{\text{total}}={\sigma }_s+{\sigma }_d$ .

From these we get that in a Kelvin-Voigt material, stress σ, strain ε, and their rates of change with respect to time t are governed by an equation of the form

$\sigma (t)=E\epsilon \left(t\right)+\eta \frac{\text{d}\epsilon \left(t\right)}{\text{d}t}$ (8)

where E is the modulus of elasticity and η is the viscosity of the artery vessel.

If we consider the artery vessel as a cylindrical tube, the resultant circumferential strain ${ϵ}_{\theta }$ takes the form ${ϵ}_{\theta }=\frac{1}{E}\left({\sigma }_{\theta }+\nu {\sigma }_z\right)$ , where ν is Poisson’s ratio.

Similarly, the resultant longitudinal strain ${ϵ}_z$ is given by ${ϵ}_z=\frac{1}{E}\left({\sigma }_z+\nu {\sigma }_{\theta }\right)$ . If

we stipulate that the artery vessel is tethered in the longitudinal direction, then ${ϵ}_z=0$ , and the circumferential strain becomes ${ϵ}_{\theta }=\frac{1}{E}\left[{\sigma }_{\theta }\left(1-{\nu }^2\right)\right]$ . Since the circumferential strain is defined as

${ϵ}_{\theta }=\frac{R-R_0}{R_0}=\frac{\sqrt{A}-\sqrt{A_0}}{\sqrt{A_0}}$

where R is the tube radius and A its corresponding area (the subscript 0 indicating conditions where no stresses are imposed), then

$\frac{1}{E}\left[{\sigma }_{\theta }\left(1-{\nu }^2\right)\right]=\frac{\sqrt{A}-\sqrt{A_0}}{\sqrt{A_0}}$ (9)

But, the circumferential stress acting in the tube wall of length l and thickness h, is due to the difference between the internal pressure $P_i$ and external pressure $P_e$ , and under the assumption that the arterial vessel wall is thin, that is, $R_e\approx R_i$ , then

${\sigma }_{\theta }\approx \frac{\left(P_i-P_e\right)R_i}{h}$ (10)

Substituting Equation (10) into Equation (9), and recognizing that $A_0=\text{π}{R}_i^2$ , yields

$P_i=P_e+\frac{\beta }{A_0}\left(\sqrt{A}-\sqrt{A_0}\right)$ (11)

where $\beta =\frac{\sqrt{\text{π}}{h}_{0}E}{\left(1-{\nu }^2\right)}$ .

Equation (11) represents the elastic behavior of the artery vessel, and corresponds to the spring in the Kelvin-Voigt model. Then, $E\epsilon \left(t\right)$ in Equation (8) is given by

$E\epsilon \left(t\right)=P_e+\frac{\beta }{A_0}\left(\sqrt{A}-\sqrt{A_0}\right)$ (12)

and

$\eta \frac{\text{d}\epsilon \left(t\right)}{\text{d}t}=\frac{\Gamma }{A_0\sqrt{A}}\frac{\partial A}{\partial t}$ (13)

where $\Gamma =\frac{h_0}{2\left(1-{\nu }^2\right)}\sqrt{\text{π}}\eta$ .

Finally, substituting Equations (12) and (13) into Equation (8) gives the resultant viscoelastic tube law as

$P_i=P_e+\frac{\beta }{A_0}\left(\sqrt{A}-\sqrt{A_0}\right)+\frac{\Gamma }{A_0\sqrt{A}}\frac{\partial A}{\partial t}$ (14)

3.2. Velocity Profile and Wall Friction

In 1-D modelling the velocity profile is commonly assumed to be constant in shape and axisymmetric. A typical profile satisfying the no-slip condition at the artery vessel wall $\left[v\left(r=R\right)=0\right]$ is (Hughes & Lubliner, 1973)

$v\left(z,r,t\right)=U\frac{\zeta +2}{\zeta }\left[1-{\left(\frac{r}{R}\right)}^{\zeta }\right]$ (15)

where r is the radial coordinate, $R\left(z,t\right)$ is the radius of the artery vessel, ζ is a velocity profile parameter, and U is the average velocity.

Integration of the incompressible Navier-Stokes equation for axisymmetric vessels yields

$f\left(z,t\right)=2\mu \text{π}R{\left[\frac{\partial v}{\partial r}\right]}_{r=R}$

where f is the frictional force per unit of length, and μ is the viscosity of the blood. For the velocity profile given by Equation (15), we have $f=-2\left(\zeta +2\right)\mu \text{π}U$ . Following Smith, Pullan, & Hunter (2002) , $\zeta =9$ provides a good compromise to experimental data obtained at different points in the cardiac cycle, which then gives

$f=-22\mu \text{π}U$ (16)

Then, Euler’s momentum equation (Equation (2)), corrected for velocity profile, with average velocity U, and wall friction f would read as

$\frac{\partial U}{\partial t}+U\frac{\partial U}{\partial z}+\frac{1}{\rho }\frac{\partial p}{\partial z}=\frac{f}{\rho A}$ (17)

which should be solved together with the continuity equation (Equation (1)), rewritten here as

$\left(\frac{\partial A}{\partial t}\right)+\frac{\partial \left(UA\right)}{\partial z}=0$ (18)

4. Numerical Solutions of Euler’s Hemodynamic Equations

For a long time, Euler’s hemodynamic model was mainly considered an historical curiosity, because practical applications of his equations require rather advanced numerical computing. It was only in the last decades of the last century that we see a rapid growth in the blood flow analysis, mainly due to the developments in digital computer simulations. In fact, the major advances have been occurring since the beginning of the current century, and it is possible to say that this is an area of considerable research interest, with several publications each year.

It is possible to solve Euler’s hemodynamic equations using different numerical schemes. These include the method of characteristics, finite volumes, finite differences, and finite elements. A comparison of these methods has shown a good agreement in their ability to capture the main features of pressure, velocity and area waveforms in large arteries (Boileau et al., 2015) . In the next section, we analyze their solution using the method of characteristics.

4.1. Method of Characteristics Analysis

Under physiological conditions, the elastic term in the tube law (Equation (14)) is dominant over the viscous term. Then, neglecting the viscous term, the pressure gradient term in the momentum equation (Equation (17)) takes the form

$\frac{1}{\rho }\frac{\partial p}{\partial z}=\frac{1}{\rho }\left(\frac{\partial p}{\partial A}\frac{\partial A}{\partial z}+\frac{\partial p}{\partial \beta }\frac{\partial \beta }{\partial z}+\frac{\partial p}{\partial A_0}\frac{\partial A_0}{\partial z}\right)$ (19)

Following Alastruey, Parker, & Sherwin (2012) , Equations (17) and (18) form a system of hyperbolic partial differential equations that can be written in non-conservative form with constant A as

$\frac{\partial U}{\partial t}+H\frac{\partial U}{\partial z}=S$ (20)

where

$U=\left[\begin{array}{c}A\\ U\end{array}\right]$ , $H=\left[\begin{array}{cc}U& A\\ \frac{1}{\rho }\frac{\partial p}{\partial A}& U\end{array}\right]$ , $S=\left[\begin{array}{c}0\\ \frac{1}{\rho }\left(\frac{f}{A}-\frac{\partial p}{\partial \beta }\frac{\text{d}\beta }{\text{d}z}-\frac{\partial p}{\partial A_0}\frac{\text{d}{A}_0}{\text{d}z}\right)\end{array}\right]$

This system can be analyzed using Riemann’s method of characteristics³. Since $A>0$ , and in normal physiological conditions $\frac{1}{\rho }\frac{\partial p}{\partial A}>0$ , then, H has two real

and distinct eigenvalues, ${\lambda }_{f,b}=U\pm c$ , where the subscript f applies to the forward travelling wave, and the subscript b applies to the backward travelling

wave, and $c=\sqrt{\frac{A}{\rho }\frac{\partial p}{\partial A}}$ .

Additionally, the matrix H is diagonalizable, since there exists an invertible matrix L such that

$H=L^{-1}\Lambda L$ (21)

where $L=\delta \left[\begin{array}{cc}\frac{c}{A}& 1\\ -\frac{c}{A}& 1\end{array}\right]$ , $\Lambda =\left[\begin{array}{cc}{\lambda }_f& 0\\ 0& {\lambda }_b\end{array}\right]$ , and δ is a scaling factor.

Substitution of Equation (21) into Equation (20), and premultiplication of Equation (20) by L yields

$L\frac{\partial U}{\partial t}+\Lambda L\frac{\partial U}{\partial z}=LS$ (22)

Taking $\frac{\partial W}{\partial U}=L$ , where $W={\left[W_f,W_b\right]}^{\text{T}}$ is the vector of characteristic variables, Equation (22) reduces to

$\frac{\partial W}{\partial t}+\Lambda \frac{\partial W}{\partial z}=LS$ (23)

For any path $x=\hat{x}\left(t\right)$ in the $\left(x,t\right)$ space, the variation of W along $\hat{x}\left(t\right)$ can be written as

$\frac{\text{d}W\left(\hat{x}\left(t\right),t\right)}{\text{d}t}=\frac{\partial W}{\partial t}+\frac{\text{d}\hat{x}}{\text{d}t}I\frac{\partial W}{\partial \hat{x}}$ (24)

Comparison of Equation (23) and Equation (24) shows that if the path $\hat{x}\left(t\right)$ is chosen such that $\frac{\text{d}\hat{x}}{\text{d}t}I=\Lambda$ , then

$\frac{\text{d}W}{\text{d}t}=LS=\left[\begin{array}{c}\frac{1}{\rho }\left(\frac{f}{A}-\frac{\partial p}{\partial \beta }\frac{\text{d}\beta }{\text{d}z}-\frac{\partial p}{\partial A_0}\frac{\text{d}{A}_0}{\text{d}z}\right)\\ \frac{1}{\rho }\left(\frac{f}{A}-\frac{\partial p}{\partial \beta }\frac{\text{d}\beta }{\text{d}z}-\frac{\partial p}{\partial A_0}\frac{\text{d}{A}_0}{\text{d}z}\right)\end{array}\right]$ (25)

Thus, for any point $\left(X,T\right)$ in the $\left(x,t\right)$ space there are two characteristic paths, $C_f$ and $C_b$ , along which $W_f$ and $W_b$ propagate at speeds ${\lambda }_f$ and ${\lambda }_b$ , respectively, changing their values due to fluid viscous dissipation and spatial variations of wall distensibility and reference luminal area.

The characteristic variables $W_f$ and $W_b$ are then determined by integration of $\frac{\partial W}{\partial U}=L$ . However, before doing that, it should be noted that to satisfy the Cauchy-Riemann condition $\frac{{\partial }^2{W}_{f,b}}{\partial A\partial U}=\frac{{\partial }^2{W}_{f,b}}{\partial U\partial A}$ , the value of δ in L must be

constant, which, without loss of generality, can be assumed to be equal to one. Then

$W_{f,b}={\displaystyle {\int }_{U_0}^U\text{d}U}\pm {\displaystyle {\int }_{A_0}^A\frac{c}{A}\text{d}A}$ (26)

where $U_0$ and $A_0$ are reference values.

For the pressure-area relationship given by Equation (14), we can write an explicit form for $W_{f,b}$ . Recalling that $c=\sqrt{\frac{A}{\rho }\frac{\partial p}{\partial A}}$ , and by neglecting the viscous dissipation term in Equation (14), we have that $c=\sqrt{\frac{\beta }{2\rho A_0}}{A}^{1/4}$ ; then, the integration of Equation (26) results in

$W_{f,b}=U-U_0\pm 4\left(c-c_0\right)$ (27)

where $c_0=\sqrt{\frac{\beta }{2\rho }}{A}_0^{-1/4}=\sqrt{\frac{E{h}_0}{\rho D_0\left(1-{\nu }^2\right)}}=\sqrt{\frac{4}{3}\frac{E{h}_0}{\rho D_0}}$ , for $\nu =0.5$ . As we saw earlier, Moens and Korteweg independently derived the equation $c_0=\sqrt{\frac{E{h}_0}{\rho D_0}}$ for

the pulse wave speed, which is similar to the equation for $c_0$ obtained from Equation (14).

This shows that $W_f$ propagates changes in area (and, hence, pressure and velocity) in the positive x-direction in an arterial segment; that is, forwards. On the other hand, $W_b$ propagates changes in the negative x-direction; that is, backwards. Therefore, blood pressure and flow rate at any point in an artery vessel may be described as the combination of forward- and backward-travelling waves.

4.2. Boundary Conditions

It is now necessary to prescribe boundary conditions at both the inlet and outlet of each arterial segment. These are often classified as inflow, junction and terminal boundary conditions.

Inflow boundary condition

This boundary condition, is usually accomplished by prescribing the volume flow rate measured in vivo, $Q_{in}\left(t\right)$ . In a healthy adult at rest, the heart rate is about 70 beats/min (bpm), giving a cardiac period of just less than 1 s. Each ventricle ejects about 70 - 100 mL of blood per stroke; the so-called stroke volume. The net volume of blood ejected from the left ventricle to the ascending aorta per unit of time, the cardiac output, is around 6 l/min (Alastruey, Parker, & Sherwin, 2012) .

Junction matching conditions

Junction matching conditions allow the connection of individual arteries to form an arterial network. There are two types of junctions: splitting flow and merging flow. Splitting flow junctions are the most common arrangement in large human systemic arteries, whereas merging flow junctions are the most common arrangement in the venous system. Energy losses in junctions are usually neglected (Alastruey, Parker, & Sherwin, 2012) .

Terminal boundary conditions

Any arterial 1-D model has to be truncated after a relatively small number of generations of bifurcations. According to Alastruey, Parker, & Sherwin (2012) , in the most peripheral vessels (small arteries, arterioles and capillaries), fluid resistance dominates over wall compliance and fluid inertia, which are both dominant in large arteries. The effect of peripheral resistance, compliance and fluid inertia on pulse wave propagation in large 1-D model arteries is commonly simulated using linear lumped parameter models (or zero-dimensional (0-D) models) coupled to the 1-D model terminal branches.

Windkessel theory was developed to explain how the pulsatile motion of blood from the heart is transformed to a continuous steady flow at the peripheral blood vessels. The 3-Element Windkessel Model simulates the characteristic impedance of the proximal aorta. The aorta is the largest artery in humans, originating from the heart’s left ventricle and extending down to the abdomen, where it branches into smaller arteries.

In the tree-element lumped parameter Windkessel model, the chamber with compliance C is considered to be filled at the inlet by a pulsatile flow of blood $Q_{in}\left(t\right)$ , and with an outflow $Q_{out}\left(t\right)$ , and resistance to flow at the inlet and outlet, $R_1$ and $R_2$ , respectively. The compliance C of the elastic Windkessel chamber with blood volume $V\left(t\right)$ , and blood pressure $p\left(t\right)$ , is defined as

$C=\frac{\text{d}V}{\text{d}p}$ , and represents the change in vessel volume for a given pressure change.

This model relates pressure and the flow at the end point of a 1-D domain through

$Q\left(1+\frac{R_1}{R_2}\right)+C{R}_1\frac{\text{d}Q}{\text{d}t}=\frac{p-p_v}{R_2}+C\frac{\text{d}p}{\text{d}t}$ (28)

where $R_1$ is given by the characteristic impedance $Z_0$ as seen by the wave travelling down the aorta, given by $Z_0=\frac{\rho c_0}{A_0}$ , $p_v$ and $R_2$ are the pressure and resistance of the peripheral vasculature.

4.3. Wave Reflections

According to Alastruey, Parker, & Sherwin (2012) , in normal conditions, the junctions in the aorta and first generation of bifurcations are close to well-matched for the propagation of waves travelling from the heart to the periphery. This means that the same junctions are necessarily poorly-matched for waves travelling from the periphery to the heart. Thus, as waves travel progressively further down the generations of bifurcations of the arterial tree, their reflections effectively become “trapped”, with an ever diminishing proportion of their amplitude making the way back to proximal arteries. This wave trapping mechanism prevents distal changes in pressure and velocity from being seen in the proximal aorta.

The reflection coefficient at the outlet of a terminal branch $R_f$ coupled to a single resistance $R_1$ is

$R_f=\frac{R_1-Z_0}{R_1+Z_0}$ . (29)

By assuming $R_1=Z_0$ yields $R_f=0$ , which is a condition often used in the terminal branches of the artery system. This also means that since the transmission coefficient T is given by $T=1+R_f$ , the amplitude of the incident wave is fully transmitted to the peripheral vasculature. Finally, it should be noted that $R_1=Z_0$ means that any incoming wave is completely absorbed by the 0-D model described in the last section.

4.4. Simulation Tests and Results

Equations (14), (17), and (18) (the refined forms of Euler’s hemodynamic equations) were tested by Alastruey, Parker, & Sherwin (2012) , using a discontinuous Galerkin scheme, by comparison against in vitro data in a 1:1 replica of the 37 largest conduit arteries made of distensible silicone tubes (Figure 1 (center)). The simulated aorta was connected to a pump, which simulates the left ventricle, and the terminal branches are connected through resistance elements to a returning circuit, which simulates the venous return. The comparisons shown in Figure 1 demonstrate the ability of the 1-D formulation to simulate pulse wave propagation in large arterial networks with reasonable accuracy.

The 1-D equations were also solved by Alastruey, Parker, & Sherwin (2012) in the 55 larger arteries in the human (Figure 2). The periodic flow rate shown in Figure 2 (bottom left) was prescribed as the boundary condition at the inlet of the ascending aorta (Segment 1) and couple RCR models to each terminal branch. They exhibit the characteristics features observed in vivo under normal conditions (for a discussion about the behavior of these curves see Alastruey, Parker, & Sherwin, 2012 ).

Figure 3 compares visco-elastic ad purely-elastic results of simulated area-pressure

Figure 1. Experimental (exp) and simulated (num) pressure and flow waveforms in the right carotid artery (left), thoracic aorta (right) and right femoral artery (bottom) of a 1:1 replica of the 37 largest systemic arteries (top middle). 1: Pump (left heart); 2: catheter access; 3: aortic valve; 4: peripheral resistance tube; 5: stiff plastic tubing (veins); 6: venous over flow; 7: venous return conduit; 8: buffering reservoir; 9: pulmonary veins (Alastruey, Parker, & Sherwin, 2012) .

Figure 2. Pressure (top) and flow rate (bottom) with time at the (left) aortic root (Root, segment 1) and midpoint of the aortic arch B (Arch, segment 14), thoracic aorta B (Tho, segment 27), and abdominal aorta D (Abd, segment 39), and (right) midpoint of the left common carotid (CCA, segment 15), left brachial (Bra, segment 21), right renal (Ren, segment 38) and left femoral (Fem, segment 46) arteries. They were simulated using a nonlinear visco-elastic 1-D model of pulse wave propagation in the larger 55 systemic arteries in the human (centre). The flow rate at the root (d) was measured in vivo and prescribed as the inflow boundary condition. In red, it is shown the uniform Windkessel pressure p_w, and Windkessel flow rate q_w, out of the arterial system into the microcirculation (Alastruey, Parker, & Sherwin, 2012) .

Figure 3. Simulated area-pressure curve (a); pressure with time (b) and flow rate with time (c) in the midpoint of the right radial artery of the 55-artery “normal young” model (Segment 12) (Alastruey, Parker, & Sherwin, 2012) .

curve, pressure with time, and flow rate with time in the midpoint of the right radial artery of the 55-artery normal young model (Segment 12), where it is seen that the area-pressure curve exhibits hysteresis (Figure 3(a)) due to wall visco-elasticity.

5. Conclusion

It has been demonstrated the pervasiveness of Euler’s hemodynamic model, and that the hemodynamic equations developed by him in 1742, about 275 years ago, still undergird the most advanced numerical methods in use today for blood flow analysis in arterial networks. At the time Euler wrote his essay, the knowledge on flow in elastic tubes had not been yet subjected to mathematical analysis, being Euler the first to propose a plausible model to the problem. Nonetheless, by not recognizing the wave nature of the hemodynamic equations, led Euler to a dead end on trying to find a closed form solution to his equations. It was only in 1759 that Euler himself obtained the wave equation and its associated general solution in an essay on the propagation of sound. The propagation of waves in an elastic tube requires an adequate constitutive relation for the viscoelastic behavior of its walls, which Euler was unable to establish at that time. Then, for more 200 years, Euler’s hemodynamic equations remain practically dormant, and it was only in the last decades of the 20th century that blood flow analysis became possible due to the advances in numerical computing. Today, because of the importance on better understanding cardiovascular diseases, blood flow in human arteries is a thriving area of research, and it is possible to say that all the particular features of arterial network can be now rather adequately modeled, which allow the simulation of pulse wave propagation in all the cardiac cycle phases in large arterial networks with very good accuracy. As noted earlier, Euler’s pioneering and seminal work in the area of blood flow justifies he be called the father of hemodynamics.

NOTES

¹Available in the Euler Archive http://eulerarchive.maa.org/.

²In 1747, D’Alembert (1747) derived the equation for the shape of an stretched string subject to vibrations in the form $y=\psi \left(t+s\right)+\Gamma \left(t-s\right)$ , without explicitly deriving the wave equation. In 1759, Euler (1759) instead, not only derived the wave equation, but also gave its solution in the form $y=\varphi \left(x+ct\right)+\psi \left(x-ct\right)$ , in an essay on the propagation of sound. Nonetheless, the general solution of the wave equation became associated with D’Alembert.

³The method of Riemann characteristics has been used for more than a century to describe linear and nonlinear waves propagating in solids, liquids, and gases. The method can be also applied to one-dimensional wave propagation in uniform and nonuniform rigid and elastic ducts.

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Alastruey, J., Parker, K. H., & Sherwin, S. J. (2012). Arterial Pulse Wave Haemodynamics. In S. Anderson (Ed.), 11th International Conference on Pressure Surges (pp. 401-443). Lisbon: Virtual PiE Led t/a BHR Group.
[2]	Boileau, E., et al. (2015). A Benchmark Study of Numerical Schemes for One-Dimensional Arterial Blood Flow Modelling. International Journal for Numerical Methods in Engineering, 31, Article ID: e02732. https://doi.org/10.1002/cnm.2732
[3]	Cerny, L. C., & Walawender, W. P. (1974). Leonhardi Euleri’s Principia pro motu sanguinis per arterias determinando. Journal of Biological Physics, 2, 41-56. https://doi.org/10.1007/BF02309372
[4]	D’Alembert, J. L. R. (1747). Recherches sur la courbe que forme une corde tendue mise en vibration. Histoire de l’Académie Royale des Sciences et Belles-Lettres de Berlin pour l’année, 1750, 214-219.
[5]	Euler, L. (1754). E 206—Sur le mouvement de l’eau par des tuyaux de conduite. Mémoires de l’académie des sciences de Berlin, 8, 111-148. (Also in Opera Omnia: Series 2, Volume 15, pp. 219-250. According to C. G. J. Jacobi, a Treatise with This Title Was Presented to the Berlin Academy on October 23, 1749).
[6]	Euler, L. (1759). E 305—De la propagation du son. Mém Acad Sci Berlin, 15, 185-209.
[7]	Euler, L. (1862). E 855—Principia pro motu sanguinis per arterias determinando. Opera Postuma, 2, 814-823.
[8]	Euler, L. (1979). E 855—Principia pro motu sanguinis per arterias determinando. Opera Omnia: Series 2, 16, 178-196.
[9]	Hughes, T. J. R., & Lubliner, J. (1973). On the One-Dimensional Theory of Blood Flow in the Larger Vessels. Math Biosciences, 18, 161-170. https://doi.org/10.1016/0025-5564(73)90027-8
[10]	Parker, K. H. (2009). A Brief History of Arterial Wave Mechanics. Medical & Biological Engineering & Computing, 47, 111-118. https://doi.org/10.1007/s11517-009-0440-5
[11]	Saito, M., et al. (2011). One-Dimensional Model for Propagation of a Pressure Wave in a Model of the Human Arterial Network: Comparison of Theoretical and Experimental Results. Journal of Biomechanical Engineering, 133, Article ID: 121005. https://doi.org/10.1115/1.4005472
[12]	Smith, N. P., Pullan, A. J., & Hunter, P. J. (2002). An Anatomically Based Model of Transient Coronary Blood Flow in the Heart. SIAM Journal on Applied Mathematics, 62, 990-1018. https://doi.org/10.1137/S0036139999355199
[13]	Tijsseling, A. S., & Anderson, A. (2012). A Isebree Moens and DJ Korteweg: On the Speed of Propagation of Waves in Elastic Tubes [Keynote Paper]. http://www.win.tue.nl/~atijssel/pdf_files/CASA-12-42.pdf
[14]	Truesdell, C. A. (1955). Leonhardi Euleri Opera. Commentationes Mechanica ad Theoriam Corporum Fluidorum Pertinentes, Volumen Posterius, Lausannae, Orell Fussli Turici (p. LXXVII).
[15]	Young, T. (1808). Hydraulic Investigations, Subservient to an Intended Croonian Lecture on the Motion of the Blood. Philosophical Transactions of the Royal Society of London, 98, 164-186. https://doi.org/10.1098/rstl.1808.0014