2 Experimental data

We base our work on the in vivo measurements of pressure $p$ and velocity $v$ in brain vessels that were taken during endovascular surgery on AVM embolisation in the Meshalkin Research Institute of Circulation Pathology (Chupakhin et al. Reference Chupakhin, Cherevko, Khe, Telegina, Krivoshapkin, Orlov, Panarin and Baranov2012). All the measurements were done with the use of Volcano ComboMap© Pressure and Flow System in different points of the afferent artery and efferent vein (i.e. feeding and draining vessels). The 0.356 mm thick Volcano ComboWire guide wire of the ComboMap system was used to measure pressure and flow velocity. Typical diameters of the vessels under consideration were 2–3 mm.

Measurements were performed for more than 40 cases, mostly for the arterial part of the vascular system. The data obtained during operation are pressure and velocity as functions of time. The typical graphs for velocity and pressure in the feeding artery and draining vein are shown in figure 1(a). Note that arterial and venous measurements were never done simultaneously, although graphs for arterial and venous vessels in figure 1(a), correspond to one patient. It is convenient to represent the data on a $vp$ -diagram (see figure 1 b). Dashed purple and dash-dotted green diagrams correspond to the measurements taken during one treatment in the afferent of the AVM before and after the embolisation. The $vp$ -diagram corresponding to the draining vein is shown in solid blue.

Figure 1.

Typical graphs of in vivo measured velocity and pressure of the blood flow in a feeding artery (afferent) and a draining vein of an AVM. (a) Time series. (b)  $vp$ -diagram in the afferent of the AVM before and after the embolisation, and in the draining vein. Arrows show the direction of movement along the diagrams with time increase.

Several typical features of the $vp$ -diagram have been described (Chupakhin et al. Reference Chupakhin, Cherevko, Khe, Telegina, Krivoshapkin, Orlov, Panarin and Baranov2012):

1. (i) characteristic semicircular shape of the diagram;

2. (ii) conservation of the diagram’s geometry for different measurement points along the vessel;

3. (iii) shift of the diagrams in the upper-left direction on $vp$ -plane in the course of the embolisation;

4. (iv) counterclockwise circulation of a point drawing the diagram in arteries during the heartbeat;

5. (v) $vp$ -diagrams in veins have a similar shape but are reflected with respect to the $v=p$ line with clockwise circulation of a drawing point during the heartbeat.

The semicircular shape of the diagram indicates the presence of the time shift between the velocity and pressure waves in the vessels. The velocity wave takes the lead over the pressure wave in the arterial part and falls behind in the venous part. The lower smooth part of the arterial diagram in figure 1(b) corresponds to the systole (contraction of the heart muscle), whereas the wavy upper part reflects the diastole (relaxation of the heart muscle). The aim of the paper is to reproduce this behaviour of the $vp$ -diagram in the numerical model.

3 How to explain the phase shift?

In the literature concerning the cardiovascular system, the phase shift between the velocity and pressure waves is described in terms of the non-dimensional Womersley number (Caro et al. Reference Caro, Pedley, Schroter and Seed1978). This number governs the relationship between the unsteady and viscous forces for a non-steady flow of viscous fluid in a cylindrical tube.

The velocity of fully developed non-stationary laminar flow of a viscous fluid driven by a time-dependent pressure gradient in a rigid cylindrical tube of fixed radius $R$ is given by the following exact solution of Navier–Stokes equations (Gromeka Reference Gromeka1882, Reference Gromeka1883; Womersley Reference Womersley1955):

Here $w$ is the only non-zero velocity component directed along the tube axis, $l$ is the tube length and $r$ is the distance from the centre of the tube. Function $W(y)$ is the real part of the complex function

where $\text{J}_{0}$ is the Bessel function of the first kind and ${\it\rho}$ is the fluid density. The non-dimensional number

is known as the Womersley number (Caro et al. Reference Caro, Pedley, Schroter and Seed1978). Here ${\it\nu}$ is the kinematic viscosity of the fluid.

According to the solution (3.1) the phase shift between the velocity and the pressure wave appears due to the non-zero argument of the complex function $W$ as it follows from the representation

For ${\it\alpha}\ll 1$ (narrow tube, high viscosity) the velocity profile is close to Poiseuille flow:

The phase shift is negligible since function $W$ is real up to the order $O({\it\alpha}^{4})$ .

For large ${\it\alpha}$ (wider tube, small viscosity) the function $W$ has a non-zero imaginary part and the phase shift becomes finite.

In our problem $R\approx 0.1$  cm, ${\it\omega}\approx 1~\text{s}^{-1}$ , ${\it\nu}\approx 3\times 10^{-2}~\text{cm}^{2}~\text{s}^{-1}$ . Hence,

This calculation implies that in the observed, relatively narrow, blood vessels the velocity profile across the tube can be averaged, as in Poiseuille flow, and the inertial effects described by Gromeka & Womersley’s solution do not seem to play a significant role.

Below, we present a model that gives another explanation of the velocity–pressure time shift and qualitatively describes all of the mentioned features of the $vp$ -diagrams.

4 Mathematical model

Let us assume that the brain tissue occupies a region denoted by ${\it\Omega}$ . The outer boundary of ${\it\Omega}$ is impermeable. The vascular system within the tissue consists of two basins: arterial and venous. There could be several input and output points at the vessel networks, however for simplicity we observe the case of one input $O$ at the arterial side and one output $E$ at the venous side.

The vascular system plays the role of inner boundary ${\it\Gamma}\subset \partial {\it\Omega}$ for the area ${\it\Omega}$ . Each vessel’s segment (from one bifurcation to another) is denoted ${\it\Gamma}_{j}$ , $j=1,\ldots ,N$ , so there are a total of $N$ segments in the vessels network and ${\it\Gamma}=\bigcup _{j=1}^{N}{\it\Gamma}_{j}$ .

In our model we explicitly resolve the vessel network only up to a certain level of detail. All the vessels of smaller size, as well as the capillary network, are modelled uniformly as a filtration flow within the porous brain tissue. It is assumed that there is a difference between the blood pressure $U$ within the vessel and the brain fluid pressure $p$ in the capillary network. Pressure $U$ is determined only over the vessel network.

5 Numerical experiment

In order to reproduce numerically the described qualitative characteristics of the cerebral blood flow we use a 2-D model of the vessel net. As shown below, the proposed model demonstrates the phase shift between the pressure and velocity waves, as well as a suitable interval for the pressures and velocities. More accurate simulations would require reconstruction of the full 3-D model of the vessels up to a certain level of detail, selection of parameters of the vessels and of the brain tissue, etc.

6 Theoretical explanation of the time shift

In order to give a theoretical explanation of the observed time shift between the velocity and pressure waves in the vessel we observe the simplest configuration of one infinite vessel within a cylindrical porous medium. The vessel of radius $R_{0}$ is located along the positive $Ox$ -axis. The porous tissue forms a cylinder of radius $R_{1}$ . On the outer wall of the tissue, the no-flow condition is specified (see figure 4).

Figure 4.

Modelling configuration for the analytical prediction of the time shift between velocity and pressure waves. A single vessel of radius $R_{0}$ is immersed in a cylindrical porous medium of radius $R_{1}$ . The fluid pressure in the vessel and in the medium are $U$ and $p$ respectively.

The mathematical statement of the problem is the following. The equation for the fluid pressure in the tissue is

Here ${\it\tau}=k/({\it\mu}m_{0}{\it\varepsilon})$ . The boundary conditions on the vessel’s wall are

The no-flow condition on the outer cylinder reads

We construct the solution in the form of a single wave packet:

where $k$ is some complex wavenumber. Generally, $k$ has a non-zero real part $\text{Re}(k)\neq 0$ . In the physically meaningful solution $\text{Re}(k)\leqslant 0$ for $x>0$ (upstream source) and $\text{Re}(k)\geqslant 0$ for $x<0$ (downstream sink). In the cerebral vessels the overall flow is directed from the arteries to the veins, hence the upstream source corresponds to the arterial part and the downstream sink describes the venous part.

Representation of the solution for $p$ in (6.4) can be written as follows:

This implies that $\text{arg}(P)/{\it\omega}$ specifies the phase shift for the wave in the porous tissue with respect to the wave in the vessel.

The fluid velocity in the vessel is given by formula (4.2) which is

Therefore, the phase shift of the velocity wave with respect to the pressure wave is determined by the $\arg (k)/{\it\omega}$ . A positive value of $\arg (k)$ implies a shift of the graph of function $v_{v}$ to the left, i.e. the velocity wave takes the lead over the pressure wave (arterial scenario). On the contrary, a negative value of $\arg (k)$ shifts the graph of $v_{v}$ to the right and the pressure wave runs ahead of the velocity wave (venous scenario).

Combining the two observations concerning signs of $\text{Re}(k)$ and $\arg (k)$ we conclude that there should be $0\leqslant \arg (k)\leqslant {\rm\pi}/2$ in the arterial part and ${\rm\pi}\leqslant \arg (k)\leqslant 3{\rm\pi}/2$ in the venous part. It remains to show the existence of solutions of the problem with wavenumber $k$ in the required sectors of the complex plane.

Substitution of representation (6.4) into (6.1) gives

Integration of the latter ordinary differential equation gives the following solution for $P$ :

Henceforth, $\text{J}_{i}$ and $\text{Y}_{i}$ , $i=0,1$ , are Bessel functions of the first and second kind. The no-flow condition (6.3) after cancellation of insufficient multiplier yields

Hence, $P$ can be represented as

with an arbitrary constant $c$ . Substitution into the boundary conditions (6.2) leads to the equations

The combination of these equations allows determination of the constant $c$ in the form

Finally, we obtain the equation for the complex wavenumber $k$ :

Note that the function $g$ depends only on $k^{2}$ , hence solutions of (6.33) are symmetrical with respect to the origin $k=0$ , i.e. there are two sets of solutions for $\text{Re}(k)>0$ and $\text{Re}(k)<0$ that coincide after rotation on ${\rm\pi}$ in the complex plane.

Equation (6.33) for the complex wavenumber $k$ was solved numerically using Wolfram Mathematica©. The following set of parameters was used:

Zero-level lines for real and imaginary parts of the function $g$ are shown in figure 5. Only half of the complex plane is shown, graphs in the lower half-plane are symmetrical with respect to the origin $k=0$ . The qualitative picture of zero-level lines in figure 5 is stable to variation of the problem parameters. Intersections of the lines give the roots of (6.33) for $k$ . From this numerical experiment one can conclude the existence of a (probably, countable) set of roots $k$ of (6.33) with $0<\arg (k)<{\rm\pi}/2$ and a symmetrical set of roots with ${\rm\pi}<\arg (k)<3{\rm\pi}/2$ . Numerical values of the several first roots of (6.33) are given in table 2. Graphs of the real and imaginary parts of function $P(r)$ for the corresponding roots $k_{j}$ are shown in figure 6.

Figure 5.

Zero-level lines for real (dashed red) and imaginary (solid green) parts of the function $g(k)$ . Intersections of the lines give the complex roots $k$ of (6.33).

Table 2.

Numerical values of some roots of (6.33).

Figure 6.

Real (dashed red) and imaginary (solid green) parts of the function $P(r)$ for some roots $k_{j}$ of (6.33): (a) $k=0.798876+0.0528743i$ ; (b) $k=0.763257+0.656872i$ ; (c) $k=0.898953+0.550532i$ ; (d) $k=1.13493+0.421902i$ ; (e)  $k=1.40215+0.351927i$ ; (f) $k=1.67104+0.290608i$ ; (g) $k=2.0127+0.255911i$ ; (h)  $k=3.24419+0.155055i$ .

From the observations above we conclude that the roots $k=k_{i}$ shown in figure 5 correspond to the upstream source with the lead of the velocity wave (arterial behaviour), whereas the symmetrical roots reflect the downstream sink with the lead of the pressure wave (venous behaviour). This remark explains the difference in the direction of circulation of the arterial and venous diagrams.

7 Discussion and conclusion

Mathematical description of blood flow is often based on the 1-D model derived from the Navier–Stokes equations by averaging over the cross-section of the blood vessel. Analytical investigations of the equations governing flows in a tube with rigid and elastic walls can be found in works by Gromeka (Reference Gromeka1882, Reference Gromeka1883) and Womersley (Reference Womersley1955). A number of mathematical models for blood flow simulation have since been developed (Caro et al. Reference Caro, Pedley, Schroter and Seed1978; Pedley Reference Pedley1980).

Models based on the mass and momentum balance laws for viscous fluid flow in a net of collapsible tubes in a 1-D approximation were proposed in Sherwin et al. (Reference Sherwin, Formaggia, Peiro and Franke2003a ,Reference Sherwin, Franke, Peiró and Parker b ), Formaggia, Quarteroni & Veneziani (Reference Formaggia, Quarteroni and Veneziani2009) and Müller and Toro (Reference Müller and Toro2014). In Reymond et al. (Reference Reymond, Merenda, Perren, Rüfenacht and Stergiopulos2007, Reference Reymond, Merenda, Perren, Rüfenacht and Stergiopulos2009, Reference Reymond, Bohraus, Perren, Lazeyras and Stergiopulos2011) a 1-D model of the cardiovascular system was developed. The authors obtained a physiologically relevant model of the entire system, including the cerebral artery network, which is coupled to a heart model. A comparison of pulse wave propagation in a 1-D nonlinear model governing blood flow in Voigt-type viscoelastic arteries with in vitro measurements was carried out by Alastruey et al. (Reference Alastruey, Khir, Matthys, Segers, Sherwin, Verdonck, Parker and Peiró2011). Being suitable for description of the blood flow in relatively large vessels, these models do not account for the flow in the surrounding tissue.

Mathematical modelling of AVM haemodynamics has been tackled in several works. Local haemodynamics in the AVM nidus is considered in works by Smith and Jones (Reference Smith and Jones2000, Reference Smith and Jones2003). Various branching configurations were studied, including one-to-many branching with different flow rates. Propagation of the embolisation agent was studied in White and Smith (Reference White and Smith2013). The authors used a two-phase flow model to simulate viscous droplet motion through a bifurcation or multifurcation point.

Global hydrodynamic properties of the AVM were studied in papers by Guglielmi (Reference Guglielmi2008), Litao, Pilar-Arceo & Legaspi (Reference Litao, Pilar-Arceo and Legaspi2012) and Telegina, Chupakhin & Cherevko (Reference Telegina, Chupakhin and Cherevko2013). Their approach was based on the electric circuit or hydraulic network analogy and allowed investigation of the stationary hydraulic characteristics of the blood flow. In particular, based on an electrical network analogy, Guglielmi (Reference Guglielmi2008) analysed the embolisation treatment procedure for two models of AVMs: a low-flow and a high-flow model. Their study shows that occlusion of the arterial feeders results in a decrease in both the arterial pressure and the flow in the nidus, while occlusion of the draining vein leads to a significant increase in blood pressure in the AVM nidus. Moreover, occlusion of the high-flow fistulous component in high-flow AVMs considerably increases the blood flow through the healthy brain.

Using the electrical network model, Litao et al. (Reference Litao, Pilar-Arceo and Legaspi2012) found the influence of the presence of the multiple compartments in the AVM on a pressure drop across the nidus. They showed that the transnidal pressure decreases with an increase in the number of the AVM compartments. In the paper by Telegina et al. (Reference Telegina, Chupakhin and Cherevko2013) different scenarios for embolisation based on a hydraulic analogy, which is similar to the electrical analogy, were studied. The general recommendations were that a treatment should start with occlusion of a fistula compartment (i.e. a shunt represented by a vessel of relatively large diameter) if present, following the embolisation of racemose AVM compartments (consisting of a chaotic net of small-diameter vessels). The above-mentioned findings agree with the common medical view on the AVM heamodynamics.

In Papapanayotou, Cherruault & De La Rochefoucauld (Reference Papapanayotou, Cherruault and De La Rochefoucauld1990), by using 1-D mass and momentum balance laws, the authors developed a model describing the influence of the AVM on the distribution of fluxes within the vessels of the Willis circle. Numerical experiments show that the volume flux in the AVM could be significantly reduced by reduction of the feeding artery and obstruction of the anterior communicating artery.

Simulation of the AVM haemodynamics was performed by Simakov et al. (Reference Simakov, Gorodnova, Chupakhin and Khe2013). The authors applied a previously developed global blood circulation model to flow in the AVM described as a graph of vessels with some arbitrary connections (Abakumov et al. Reference Abakumov, Ashmetkov, Esikova, Koshelev, Mukhin, Sosnin, Tishkin, Favorskii and Khrulenko2000; Kholodov Reference Kholodov, Belotserkovskii and Kholodov2001). The governing system was a 1-D hyperbolic system of two differential equations: conservation laws for mass and momentum. The compliance of the vessel wall was incorporated via the ‘equation of state’: a relation between the pressure and the cross-sectional area. The embolisation procedure was modelled by an increase of the hydraulic resistance of the afferent (i.e. incoming) vessels.

One-dimensional models were successfully applied to both arterial and venous vessel networks. In Qureshi et al. (Reference Qureshi, Vaughan, Sainsbury, Johnson, Peskin, Olufsen and Hill2014) such a model was used for a description of the pulmonary circulation. Blood flow in ‘large’ vessels up to the fourth generation from the inlet was governed by a 1-D model, while the mirror terminal ends of the arterial and venous nets were connected through a vessel bed. Exploiting the smallness of the vessels and linearising the 1-D governing equation, the authors constructed $2\times 2$ admittance matrices relating the flows and pressures at the outlet of the arterial and the inlet of the venous networks. Thus, the blood flow through the capillary bed was effectively resolved by relations at the terminal ends of the arterial and venous trees. We note that there should be a one-to-one mapping between corresponding terminal ends.

An approach very close to ours was applied by Cattaneo and Zunino (Reference Cattaneo and Zunino2014) for modelling of the flow in capillaries and tissue interstitium. Their model comprised a 1-D network of microvessels immersed in a 3-D porous interstitial medium. In contrast to our work, the authors addressed the problem of microcirculation and interstitial flow interaction, including lymphatic drainage. Keeping in mind mass transport and drug delivery, the authors showed the ability of the model to describe the flow patterns in healthy and tumorous tissues (i.e. in 3-D domains of the problem’s geometry). Mathematical correctness of a similar problem with the coupling of 3-D and 1-D flows was studied by D’Angelo and Quarteroni (Reference D’Angelo and Quarteroni2008), where the well-posedness of the coupled problem was established in suitable functional spaces. This allows one to perform a convergence analysis of finite element schemes.

In our paper we propose a model of blood propagation within two trees of relatively narrow vessels (arterial and venous) connected by a capillary network. The model is applicable to the numerical simulation of local cerebral haemodynamics and, in particular, can be used for modelling of the AVM anomaly. Numerical simulation shows qualitative agreement with in vivo endovascular measurements of the flow velocity and the fluid pressure. The experimentally observed phase shift between the velocity and pressure waves is reproduced in the numerical experiments. Moreover, the numerical model predicts the proper order of the waves: the lead of the velocity wave in the arterial system and of the pressure wave in the venous system.

The AVM is modelled as a shunt between the arterial and venous vessel trees. The numerical experiments demonstrate the local equalization of pressure in the arterial and venous parts, which can be treated as a possible cause of rupture. In the course of the embolisation process, pressure in the afferent of AVM increases, whereas the velocity decreases. This feature is also observable in the numerical experiments, where the embolisation is modelled as a decrease of the shunt radius.

The phase shift is investigated theoretically by construction of an analytical solution in the form of a single wave packet to a simplest model consisting of one vessel surrounded by a capillary network. Analysis of the dispersion relation reveals the existence of two sets of roots corresponding to the waves in the arterial and venous parts with the proper time shifts. Thus, the presented model suggests another (in comparison to the inertia-based one described by the Gromeka & Womersley solution) explanation of the time shift between the velocity and the pressure waves in a relatively narrow blood vessel within a porous tissue. The model can further be used for simulation of various scenarios of AVM treatment and the solution of similar problems regarding blood flow within a network of narrow blood vessels connected by a capillary network.