Introduction

Arterial beating is the fluid mechanical signature of human life’s beginning and end. Pressure pulse in arterial vascular system is also one of the first clinical sign used from ancient time to present. Not only permitting the qualitative assessment of hypertension from pulse’s hardening, but also subject’s aging from pulse’s damping and speed in line with Sydenham’s comment that “Man is as old as his arteries”. Nevertheless, beside many indicators and metrics can be extracted from arterial pulse beating signal for clinically used — e.g amplitude, wave speed, variability, etc.1,2,3, clear, physiological-based, interpretations of this complex signal are difficult to decipher4. This is due to the complexity of the underlying physiology and mechanics involved. Pressure pulse wave propagation in arteries is coupled with an extensional elastic stress wave within their walls, a phenomenon called Fluid-Structure-Interaction (FSI) in water-hammer literature5,6. One possible approach to analyze FSI’s effects are to solve them using direct 3D finite element numerical computations, for exemple using Arbitrary Lagrangian-Eulerian (ALE) formulation7. Another one consists in considering simplified one-dimensional approximation of them6,8,9. One dimensional pressure pulse wave model is indeed popular—e.g10,11,12,13,14,15— either approximating or even disregarding12,14,16,17 FSI couplings which however become significant when solid walls are soft and thin as in vascular networks5,6,18.

The rationale behind one-dimensional model is a long-wavelength asymptotic approximation recently clarified19. FSI effects are associated with two waves having distinct speeds: a fast elastic wave within the arterial wall and a slower pressure pulse wave propagation in the blood. These coupled waves are partially transmitted and reflected at each bifurcation of the arterial network, leading to a myriad of echoes interfering together. Waves damping also arises from both wave dissipation within boundary layers19, but also from viscoelastic dissipation within the solid from FSI coupling4,18,20 adding additional complexity. Finally, the heart being at the pulse’s origin, its cardiac valve closure dynamics also brings a specific influence onto the observed signal.

This structural, mechanical and geometrical complexity has been a hindrance for both modeling and interpretation of blood pressure pulse waves. Hence, being able to derive simplified, and relevant modeling in this context is still a challenge at the present state of the art. Let-us now introduce in more detail the one-dimensional long-wavelength FSI wave formulation developed in refs. 5,6,8,21 (among others). Considering an elastic tube —i.e., a vessel— of internal radius \({R}_{0}^{\star }\), thickness e⋆, permits to define the thickness to radius ratio \(\alpha ={e}^{\star }/{R}_{0}^{\star }\)8. found the two main velocities of interest for wave propagation: the pulse wave velocity \({c}_{p}^{\star }\) and the solid wall elastic velocity \({c}_{s}^{\star }\)

$${c}_{p}^{\star 2}=\frac{\frac{{K}^{\star }}{{\rho }_{f}^{\star }}}{1+\frac{2{K}^{\star }}{\alpha {E}^{\star }}\left(\frac{2(1-{\nu }^{2})}{2+\alpha }+\alpha (1+\nu )\right)},\,\,{c}_{s}^{\star 2}=\frac{{E}^{\star }}{{\rho }_{s}^{\star }}.$$
(1)

where K⋆ is blood’s isentropic compressibility modulus, ρf the blood density, E⋆ the arterial wall elastic Young modulus and ν its Poisson coefficient. In biomechanical context where 2K⋆/αE⋆ ≫ 1 and α ≪ 1, \({c}_{p}^{\star 2}\) in (1) simplifies to the Moens-Korteweg wave velocity4,18 given by,

$${c}_{p,MK}^{\star 2}=\frac{\alpha {E}^{\star }}{2{\rho }_{f}^{\star }(1-{\nu }^{2})}$$
(2)

The coupled wave propagating problem involves the velocity ratio \({{\mathcal{C}}}_{s}\) and dimensionless density ratio

$${{\mathcal{C}}}_{s}=\frac{{c}_{s}^{\star }}{{c}_{p}^{\star }},\,\,\,\,\,\,\,\,\,\,{\mathcal{D}}=\frac{{\rho }_{f}^{\star }}{{\rho }_{s}^{\star }}.$$
(3)

Reference 8 also established the four-equations long-wavelength FSI waves formulation which details the coupled, first order, hyperbolic equations for longitudinal fluid velocity W*, pressure P*, longitudinal solid deformation velocity \({\dot{\zeta }}^{* }\) and longitudinal stress \({\sigma }_{zz}^{* }\equiv {\sigma }^{* }\) as

$${\partial }_{{t}^{\star }}{W}^{* }=-\frac{1}{{\rho }_{f}^{\star }}{\partial }_{{z}^{\star }}{P}^{* },$$
(4)
$${\partial }_{{z}^{\star }}{W}^{* }+\frac{1}{{\rho }_{f}{c}_{p}^{2}}{\partial }_{{t}^{\star }}{P}^{* }=2\nu {\partial }_{{z}^{\star }}{\dot{\zeta }}^{* },$$
(5)
$${\partial }_{{t}^{\star }}{\dot{\zeta }}^{* }=\frac{1}{{\rho }_{s}}{\partial }_{{z}^{\star }}{\sigma }^{* },$$
(6)
$${\partial }_{{t}^{\star }}{\sigma }^{* }-E{\partial }_{{z}^{\star }}{\dot{\zeta }}^{* }=\frac{2\nu }{\alpha \left(2+\alpha \right)}{\partial }_{{t}^{\star }}{P}^{* },$$
(7)

where subscript * refers to dimensional quantities. The dimensionless version of Eqs. (4)–(7) have been considered in ref. 9, where the velocity perturbation is set as \({W}_{0}^{\star }\), so that \({W}^{* }={W}_{0}^{\star }W\), the pressure perturbation P*, is Joukowski’s22 one, \({P}^{* }={\rho }_{f}^{\star }{c}_{p}^{\star }{W}_{0}^{\star }P\), axial perturbed stress, σ* scales as P*, i.e., \({\sigma }^{* }={\rho }_{f}^{\star }{c}_{p}^{\star }{W}_{0}^{\star }\sigma\) from stress continuity at the tube wall. The axial displacement field, ζ*, is set as \({\zeta }^{* }=\alpha {\mathcal{M}}{L}^{\star }\zeta\) to ensure both axial velocity matching at the tube’s wall9,20, where \({\mathcal{M}}={W}_{0}^{\star }/{c}_{p}^{\star }\) is the Mach number, and L⋆ a longitudinal reference length. Scaling time t⋆ on the advection time, i.e \({t}^{\star }={L}^{\star }\tau /{c}_{p}^{\star }\) and axial longitudinal scale as z⋆ = L⋆z, dimensionless Eqs. (4)–(7) reads

$${\partial }_{\tau }W=-{\partial }_{z}P,$$
(8)
$${\partial }_{\tau }P+{\partial }_{z}W=2\alpha \nu {\partial }_{z}\dot{\zeta },$$
(9)
$${\partial }_{\tau }\dot{\zeta }=\frac{{\mathcal{D}}}{\alpha }{\partial }_{z}\sigma ,$$
(10)
$${\partial }_{\tau }\sigma -\frac{\alpha {{\mathcal{C}}}_{s}^{2}}{{\mathcal{D}}}{\partial }_{z}\dot{\zeta }=\frac{2\nu }{\alpha (2+\alpha )}{\partial }_{\tau }P.$$
(11)

This dimensionless long-wavelength FSI wave formulation is valid in a frequency range given by the Korteweg cut-off frequency \({f}_{Kc}^{\star }\) given by9,

$${f}_{Kc}^{\star }=\frac{{{\mathcal{C}}}_{s}}{{\epsilon }_{0}}{f}_{0}^{\star }\,\,\,\text{, with,}\,\,\,{f}_{0}^{\star }=\frac{{c}_{p}^{\star }}{2\pi {L}^{\star }}.$$
(12)

with ϵ0 = R0/L⋆. Eqs. (8)–(11) can be re-casted into a coupled waves propagation problem, associated with the pressure/longitudinal stress vector \({\bf{P}}=\left(\begin{array}{r}P\\ \sigma \end{array}\right)\)9,

$$\left({\partial }_{\tau }^{2}-{{\bf{C}}}_{{\bf{P}}}^{2}{\partial }_{z}^{2}\right){\bf{P}}=0,$$
(13)

with \({{\bf{C}}}_{{\bf{P}}}^{2}\) the coupled-waves matrix

$${{\bf{C}}}_{{\bf{P}}}^{2}=\left(\begin{array}{rc}1&2\nu {\mathcal{D}}\\ \frac{2\nu }{\alpha (2+\alpha )}&\frac{4{\nu }^{2}{\mathcal{D}}}{\alpha (2+\alpha )}+{{\mathcal{C}}}_{s}^{2}\end{array}\right).$$
(14)

Eigenvalues \({c}_{\pm }^{2}\) of the \({{\bf{C}}}_{{\bf{P}}}^{2}\) matrix correspond to coupled vibrating modes wave speeds propagation being solution of the polynomial characteristic problem

$${c}_{\pm }^{4}-\left[1+{{\mathcal{C}}}_{s}^{2}+\frac{4{\nu }^{2}{\mathcal{D}}}{\alpha (2+\alpha )}\right]{c}_{\pm }^{2}+{{\mathcal{C}}}_{s}^{2}=0,$$
(15)

the solution of which reads19

$${c}_{\pm }^{2}=\frac{1+{{\mathcal{C}}}_{s}^{2}+\frac{4{\nu }^{2}{\mathcal{D}}}{\alpha (2+\alpha )}\pm \sqrt{{\left(1+{{\mathcal{C}}}_{s}^{2}+\frac{4{\nu }^{2}{\mathcal{D}}}{\alpha (2+\alpha )}\right)}^{2}-4{{\mathcal{C}}}_{s}^{2}}}{2}.$$
(16)

As in ref. 9, this paper considers the diagonal base of matrix \({{\bf{C}}}_{{\bf{P}}}^{2}\)

$${\boldsymbol{\Pi }}=\left(\begin{array}{rc}\frac{2\nu {\mathcal{D}}}{{c}_{-}^{2}-1}&\frac{2\nu {\mathcal{D}}}{{c}_{+}^{2}-1}\\ 1&1\end{array}\right),\,\,{{\boldsymbol{\Pi }}}^{-1}=\frac{1}{\det \left({\boldsymbol{\Pi }}\right)}\left(\begin{array}{rc}1&-\frac{2\nu {\mathcal{D}}}{{c}_{+}^{2}-1}\\ -1&\frac{2\nu {\mathcal{D}}}{{c}_{-}^{2}-1}\end{array}\right),\,\,$$
(17)

with

$${{\mathcal{C}}}_{{\mathcal{P}}}^{2}={{\boldsymbol{\Pi }}}^{-1}\cdot {{\bf{C}}}_{{\bf{P}}}^{2}\cdot {\boldsymbol{\Pi }}=\left(\begin{array}{rc}{c}_{-}^{2}&0\\ 0&{c}_{+}^{2}\end{array}\right),$$
(18)

and

$${\bf{P}}={\boldsymbol{\Pi }}\cdot {\mathcal{P}},\,\,\,\,{\rm{with}}\,\,\,\,{\mathcal{P}}=\left({{\mathcal{P}}}_{-},{{\mathcal{P}}}_{+}\right),$$
(19)

In this paper, a new theoretical framework for solving FSI coupled waves propagation within complex elastic networks using long-wavelength approximation is developed. Section Method first discusses the long-wavelength approximation validity inside human arterial system. Section Method then establishes the coupling conditions arising at network’s bifurcation for FSI-waves within an elastic deformable network within the long-wavelength approximation. Section Result first provides the theoretical solution for these FSI waves to be decomposed into network standing waves modes generalizing the quantum graph solutions of ref. 23. Finally, Section Result then discusses and illustrates the FSI wave solutions within Willis circle network cerebral arteries network, for which the vibration spectrum is evaluated. Even if wave damping is not considered in this paper, extensions of the approach toward taking into account damping mechanisms10,11 are nevertheless provided in the discussion, also including viscoelastic effects.

Methods

Asymptotic validity for long-wavelength blood-hammer FSI wave modeling

Reference 19 found that the relevant dimensionless parameters associated with FSI wave propagation within elastic tubes are the dimensionless boundary layer thickness δ — also called the water-hammer parameter in the literature16, aspect ratio \({\epsilon }_{0}={R}_{0}^{* }/{L}^{* }\), Reynolds number Re, pulsed Reynolds number Rep and Mach number \({\mathcal{M}}\) defined as

$$\begin{array}{ll}R{e}_{p}=\frac{{c}_{p}^{\star }{R}_{0}^{\star }}{{\nu }_{f}^{\star }}\gg 1,&Re=\frac{{W}_{0}^{\star }{R}_{0}^{\star }}{{\nu }_{f}^{\star }}={\mathcal{M}}R{e}_{p},\end{array}$$
(20)
$$\begin{array}{ll}{\delta }^{2}=\frac{{\nu }_{f}^{\star }{L}^{\star }}{{c}_{p}^{\star }{R}_{0}^{\star 2}}=\frac{1}{{\epsilon }_{0}R{e}_{p}}\ll 1,&{\mathcal{M}}=\frac{{W}_{0}^{\star }}{{c}_{p}^{\star }}\ll 1,\end{array}$$
(21)

where \({\nu }_{f}^{\star }\) is blood’s kinematic viscosity. It is interesting to note that, defining the reference longitudinal length-scale L* as the propagative length over the cardiac cycle period T⋆, i.e., \({L}_{p}^{* }={T}^{* }{c}_{p}^{\star }\), one finds that \({\delta }^{2}={\nu }_{f}^{\star }{L}_{p}^{\star }/{c}_{p}^{\star }{R}_{0}^{\star 2}=1/{{\mathcal{W}}}_{0}\) with \({{\mathcal{W}}}_{0}={{R}_{0}^{\star }}^{2}/{\nu }_{f}^{\star }{T}^{\star }\) the classical Womersley dimensionless parameter used in pulsatile flows. In the context of steady flow within deformable elastic pipe networks, ref. 19 established long-wavelength FSI waves under asymptotic conditions

$${\delta }^{2}\gg {\mathcal{M}} > \frac{{\mathcal{M}}}{{{\mathcal{C}}}_{s}^{2}},\,\,\,\,\,\,\,\,\,\,\,\,\,\,\delta \gg {\epsilon }_{0}^{2},$$
(22)
$$\,\,\,\,\,\,\,\,\,\,\,\,\,\,\delta \gg \alpha {\mathcal{M}},\,\,\,\,\,\,\,\,\,\,\,\,\,\,1\gg {\epsilon }_{0}\gg \alpha {\mathcal{M}}.$$
(23)

As discussed in ref. 18, Eqs. (22)–(23) also provides conditions for which the wave propagation is decoupled from the base flow. For non-stationary base flows —as the one resulting from cardiac beating— the same asymptotic analysis can be applied since the flow acceleration \({\rho }_{f}^{\star }{\partial }_{{t}^{\star }}{W}^{\star }\) is taken into account in Eqs. (4) and (8). For α ≈ 0.0536, ν = 0.495, E⋆ = 72 kPa37, \({\rho }_{f}^{\star }\approx 1060\,kg.{m}^{-3}\), using (2) and (1), one finds \({c}_{p,MK}^{\star }=1.5\,m/s\), \({c}_{s}^{\star }=8.25\,m/s\) and thus \({{\mathcal{C}}}_{s}\approx 5.66\). From (16), one can find explicitly

$${c}_{+}^{2}+{c}_{-}^{2}=1+{{\mathcal{C}}}_{s}^{2}+\frac{4{\nu }^{2}{\mathcal{D}}}{\alpha (2+\alpha )}$$
(24)

so that using \({{\mathcal{C}}}_{s}\approx 5.66\), ν = 0.495, \({\mathcal{D}}=0.95\), one finds \({c}_{+}^{2}+{c}_{-}^{2}=42.12\), \({c}_{-}^{2}=0.135\) and \({c}_{+}^{2}=41.98\), resulting in pressure pulse wave velocity \({c}_{-}{c}_{p,MK}^{\star }\approx 0.55m/s\) and a longitudinal stress wave one \({c}_{+}{c}_{p,MK}^{\star }\approx 9.71m/s\). The latter is within the range 7–23 m/s, of in-vivo reported pulse vessel’s wall wave velocity varying by merely a factor two from one vessel to another (from 7.7 m/s aorta to 20.4 m/s in coronary arteries)38. Using T⋆ = 1s, \({c}_{p,MK}^{\star }=1.5\,m/s\), one finds \({L}_{p}^{\star }={T}^{\star }{c}_{p,MK}^{\star }=1.5m\), which is larger than the maximum length of the considered cerebral arterial system, i.e \({L}_{p}^{\star } > {L}^{\star }\). Hence, for this reason, the actual reference length of Willis circle L⋆ = 0.2m28 is chosen instead.

From (12), using a typical aspect ratio ϵ0 = 1/100 and the previously evaluated parameters, permit to find the Korteweg cut-off frequency \({f}_{Kc}^{\star }=550{f}_{0}^{\star }\). Furthermore using a typical network length L⋆ = 0.2m 28, one finds from (12) that \({f}_{0}^{\star }=1.2\)Hz, so that the possible relevant range of frequencies for cerebro-vascular vessels lies between 1.2 to 770Hz. Nevertheless, the relevance of the presented analysis is more stringently enforced by the asymptotic long-wavelength condition \(\epsilon \equiv k{k}_{0}^{\star }{R}_{0}^{\star }=2\pi f{f}_{0}^{\star }{R}_{0}^{\star }/{c}_{p}^{\star }\ll 1\) for the wave-continuity conditions at bifurcations examined in next Section Method then establishes. From (12) one realizes that ϵ = fϵ0 ≃ f10−2. Hence, taking ϵ = 0.1 as long-wavelength condition shortens the previous available range [1.2–770]Hz to [1.2–12]Hz for relevant physiological frequency range of the hereby developed analysis. It will nevertheless be observed in Section Result then discusses, that this frequency range is the most relevant one for the spectral content of physiological observations.

Let-us now discuss the asymptotic conditions Eqs. (22)–(23) given the physiological parameters of cerebral arteries depicted in Fig. 1A found in the literature and provided in Table 1.

Fig. 1: Willis circle's metric graph description.
figure 1

A Anatomical sketch of Willis circle arterial network using label’s numbering of table B following28. B Section S⋆ and lengths L⋆ of the various Willis circle vessels of A, C from ref. 28. C Vectorization of the Willis circle’s of A with metric graph edge numbering following28. D Node numbering of metric graph C. Black nodes are those where homogeneous Kirchoff conditions are applied whereas colored nodes are those where Kirchoff–Robin boundary conditions are used. The blue node #1 represents downstream networks associated with micro-circulatory modeled by Windkessel-like RCR boundary condition. Node #12, 13, 14, 15 are triggering nodes associated with the cardiac-valve dynamical closure.

Full size image
Table 1 Various physiological parameters from refs. 28 and 49
Full size table

Using \({c}_{p,MK}^{\star }=1.5\,m/s\), \({\nu }_{f}^{\star }=2.51{0}^{-6}\)SI, L⋆ = 2.10−1m, \({R}_{0}^{\star }=2.1{0}^{-3}\,{\rm{m}}\), one finds δ = 0.28, i.e., δ2 = 7.8410−2 which fulfills the first asymptotic condition (22), i.e \({\delta }^{2}\gg {\mathcal{M}}\) (since, the previous estimate of \({{\mathcal{C}}}_{s}\approx 5.66\) leads to the last condition automatically satisfied when the first does) given the Mach number figures of the last column of Table 1. The second asymptotic condition of (22) is also fulfilled by the reference aspect ratio of each vessels ϵ0 = 10−2, so that δ = 0.28 is much larger than \({\epsilon }_{0}^{2}=1{0}^{-4}\). The third asymptotic conditions (23) are also satisfied by the Mach numbers M of the last column of Table 1, also taking α = 0.05 as a reference dimensionless vessel thickness. The fourth asymptotic conditions of (23) is also satisfied by ϵ0 = 10−2 for the Mach numbers of Table 1. We now consider how junction coupling should be written for long-wavelength FSI waves propagation within networks.

Long-wavelength FSI waves continuity conditions within networks

When a wave train encounters an abrupt discontinuities inside a rigid channel (such as abrupt section variations) some localized scattering modes having evanescent longitudinal decay are excited when the wavelength is comparable with the typical transverse length-scale of the pipe system39,40,41,42,43,44. On the contrary, when the wavelength is very large compared to the typical transverse length, those modes are weakly excited, providing only a small correction to the amplitude and phase of the incident wave45,46. Asymptotically, in the long-wavelength approximation, this correction zeroes, and some effective continuity conditions can be found for connecting wave trains propagating into various channels/pipes. As shown in ref. 45, the long-wavelength approximation for acoustic wave propagation within networks leads to the continuity of pressure and zero total acoustic flux. More precisely45 has examined a T-junction between three two-dimensional channels. Matched asymptotic approximation45,46, can be used to provide small parameter \(\epsilon \equiv {k}^{\star }{R}_{0}^{\star }\) corrections to this Kirchoff long-wavelength continuity condition.

Here a long-wavelength approximation is derived in order to find the effective FSI continuity equations to connect long-wavelength waves inside deformable finite elastic walls connected at some junction exemplified in Fig. 2. The derivation is performed in any dimension (2D or 3D) for pipes having any transverse section. The result also generalizes to non rectilinear vessels/pipes provided that their local curvature remains small everywhere. More precisely as far as the ratio between the mean longitudinal radius of curvature \({R}_{c}^{\star }\) and the transverse one \({R}_{0}^{\star }\) is comparable with the long-wavelength approximation parameter \({R}_{0}^{\star }/{R}_{c}^{\star } \sim \epsilon \equiv {k}^{\star }{R}_{0}^{\star }\), neglecting both being of same order is self-consistent. Fast curvature changes such as those analyzed in ref. 47 are thus not considered here. The constitutive equations for the acoustic wave inside the fluid are mass and momentum conservation

$$\frac{1}{{\rho }_{f}^{\star }{c}_{p}^{\star 2}}{\partial }_{{t}^{\star }}{p}^{\star }=-{\nabla }^{\star }\cdot {{\bf{v}}}^{\star }$$
(25)
$${\rho }_{f}^{\star }{\partial }_{{t}^{\star }}{{\bf{v}}}^{\star }=-{\nabla }^{\star }{p}^{\star },$$
(26)

whereas those for the acoustic waves within the solid are Lamé-Clapeyron equations associated with the displacement vector field ζ and the stress tensor σ⋆ read

$${\rho }_{s}^{\star }{\partial }_{{t}^{\star }}^{2}{\zeta }^{\star }={\nabla }^{\star }\cdot {\sigma }^{\star }.$$
(27)

The boundary conditions between the fluid and the solid regions are normal stress continuity and kinematic continuity

$${p}^{\star }{| }_{\partial {\Omega }_{fs}}={\bf{n}}\cdot {\sigma }^{\star }\cdot {\bf{n}}{| }_{\partial {\Omega }_{fs}}$$
(28)
$${{\bf{v}}}^{\star }\cdot {\bf{n}}{| }_{\partial {\Omega }_{fs}}={\partial }_{{t}^{\star }}{\zeta }^{\star }\cdot {\bf{n}}{| }_{\partial {\Omega }_{fs}}.$$
(29)

Integrating the time derivative of (25) inside the small connecting region \({\Omega }_{f}^{\epsilon }\) (Cf Fig. 2c), whilst using the divergence theorem, leads to

$${\partial }_{{t}^{\star }}{\int}_{{\Omega }_{f}^{\epsilon }}{\nabla }^{\star }\cdot {{\bf{v}}}^{\star }=\sum _{{e}_{k}\in {\mathcal{N}}({v}_{i})}{\partial }_{{t}^{\star }}{q}_{k}^{\star }+{\int}_{\partial {\Omega }_{fs}^{\epsilon }}{\partial }_{{t}^{\star }}{{\bf{v}}}^{\star }\cdot {\bf{n}}$$
(30)
$$=\frac{1}{{\rho }_{f}^{\star }{c}_{p}^{\star 2}}{\partial }_{{t}^{\star }}^{2}{\int}_{{\Omega }_{f}^{\epsilon }}{p}^{\star }$$
(31)

where \({\mathcal{N}}({v}_{i})\) refers to the set of edge’s neighbors of metric graph vertex vi, as the one represented in Fig. 2B, \({q}_{k}^{\star }={\int}_{\partial {\Omega }_{f}^{k}}{{\bf{v}}}^{\star }\cdot {\bf{n}}\) is the acoustic wave flux at ek edge’s border \(\partial {\Omega }_{f}^{k}\). Using boundary condition (29) in (30) then leads to

$$\sum _{{e}_{k}\in {\mathcal{N}}({v}_{i})}{\partial }_{{t}^{\star }}{q}_{k}^{\star }={\partial }_{{t}^{\star }}^{2}\left[-{\int}_{\partial {\Omega }_{fs}^{\epsilon }}{\zeta }^{\star }\cdot {\bf{n}}+\frac{1}{{\rho }_{f}^{\star }{c}_{p}^{\star 2}}{\int}_{{\Omega }_{f}^{\epsilon }}{p}^{\star }\right]$$
(32)

Since both \({\Omega }_{f}^{\epsilon }\) volume and \({\Omega }_{s}^{\epsilon }\) surface tends to zero in the limit ϵ → 0, it is easy to find that the r.h.s of (32) cancels out so that the long-wavelength edge continuity condition is related to the acoustic flux conservation (e.g., ref. 45)

$$\sum _{{e}_{k}\in {\mathcal{N}}({v}_{i})}{\partial }_{{t}^{\star }}{q}_{k}^{\star }({{\bf{x}}}_{{v}_{i}}^{\star })=0.$$
(33)

Furthermore, since in the long-wavelength approximation the pressure p is uniform along any transverse direction (radial direction if the vessel’s section is circular) (33) then, using (26) in (33) leads to the well-known23 Kirchoff condition for the fluid pressure

$$\sum _{{e}_{k}\in {\mathcal{N}}({v}_{i})}{S}_{k}^{\star l}\frac{d{p}_{{e}_{k}}^{\star }}{d{x}^{\star }}({{\bf{x}}}_{{v}_{i}}^{\star })=0,$$
(34)

where \({S}_{k}^{\star l}\) is the transverse lumen section of vessel/pipe edge ek. Integrating now (27) in \({\Omega }_{s}^{\epsilon }\), using the divergence theorem, using a stress-free external condition on ∂Ωse, i.e., \({\sigma }^{\star }\cdot {\bf{n}}{| }_{\partial {\Omega }_{se}}=0\), leads to

$${\int}_{\partial {\Omega }_{s}^{\epsilon }}{\sigma }^{\star }\cdot {\bf{n}}={\int}_{\partial {\Omega }_{fs}^{\epsilon }}{p}^{\star }+\sum _{{e}_{k}\in {\mathcal{N}}({v}_{i})}{\phi }_{\sigma }^{\star k}={\rho }_{s}^{\star }{\partial }_{{t}^{\star }}^{2}{\int}_{{\Omega }_{s}^{\epsilon }}{\zeta }^{\star }$$
(35)

where \({\phi }_{\sigma }^{\star k}={\int}_{\partial {\Omega }_{s}^{k}}{\sigma }^{\star }\cdot {\bf{n}}\) is the stress flux at ek edge’s border \(\partial {\Omega }_{s}^{k}\). Using again the fact that both \(\partial {\Omega }_{fs}^{\epsilon }\) surface and \({\Omega }_{s}^{\epsilon }\) volume zeros in the ϵ → 0 limit, leads to the long-wavelength continuity condition

$$\sum _{{e}_{k}\in {\mathcal{N}}({v}_{i})}{\phi }_{\sigma }^{\star k}({{\bf{x}}}_{{v}_{i}}^{\star })=0$$
(36)

Finally, using the asymptotic result that in the long-wavelength approximation, the longitudinal stress \({\sigma }_{zz}^{\star }\equiv {\sigma }^{\star }\) for pipe breathing axisymmetric water-hammer wave propagation is uniform along any transverse direction (radial direction if the vessel’s section is circular) within each vessel segment5,19,48, (36) leads to

$$\sum _{{e}_{k}\in {\mathcal{N}}({v}_{i})}{S}_{k}^{\star s}{\sigma }^{\star }({{\bf{x}}}_{{v}_{i}}^{\star })=0,$$
(37)

where \({S}_{k}^{\star s}\) is the solid transverse section within vessel/pipe edge ek. Finally since at each bifurcation the wave continuity conditions are associated with normal stress continuity (34) and (37) are complemented by the trivial limit of continuity conditions for the fluid pressure and the solid longitudinal stress \(\forall {e}_{k},{e}_{k{\prime} }\in {\mathcal{N}}({v}_{i})\)

$${p}_{{e}_{k}}^{\star }({{\bf{x}}}_{{v}_{i}})={p}_{{e}_{k{\prime} }}^{\star }({{\bf{x}}}_{{v}_{i}})\,\,\,\,\,\,{\sigma }_{{e}_{k}}^{\star }({{\bf{x}}}_{{v}_{i}}^{\star })={\sigma }_{{e}_{k{\prime} }}^{\star }({{\bf{x}}}_{{v}_{i}}^{\star })$$
(38)
Fig. 2: Normal stress and acoustic pressure flux continuity at bifurcations.
figure 2

A Sketch of a bifurcation between three vessel/pipe’s segments. B One-dimensional long-wavelength metric graph representation of 3D vessels bifurcation A, in the ϵ → 0 limit describing the connection between the three vessels segments represented by three metric edges for the FSI two-wave. C Inside each bifurcation, definition of the fluid and solid domains as well as their interfaces indexed by small parameter ϵ.

Full size image

Results

Quantum graph FSI Waves solutions

The vascular network is described as a compact metric graph \({\mathcal{G}}({\mathcal{V}},{\mathcal{E}})\), having vertex set \({\mathcal{V}}\), edge set \({\mathcal{E}}\) and adjacency symmetric matrix A. Here, the analysis of ref. 23 is adapted and generalized to FSI waves, but more information about operators on compact metric graph called “quantum graphs” can be found in refs. 24,25. Along each edge ij linking vertex i and j, following the curvilinear coordinate z along center-line, one defines the eigenmode solutions of the one-dimensional Laplacian \({\Psi }_{ij}^{\pm }(z)\) as

$$-{c}_{\pm }^{2}{\partial }_{z}^{2}{\Psi }_{ij}^{\pm }(z)={\lambda }^{2}{\Psi }_{ij}^{\pm }(z),$$
(39)

Defining function \({\varphi }^{\pm }(z)=\sin (\lambda z/{c}_{\pm })\), given (39), \({\Psi }_{ij}^{\pm }(z)\) can be adequately chosen as

$${\Psi }_{ij}^{\pm }(z)=\frac{{A}_{ij}}{{\varphi }^{\pm }({\ell }_{ij})}\left({\phi }_{i}^{\pm }{\varphi }^{\pm }({\ell }_{ij}-z)+{\phi }_{j}^{\pm }{\varphi }^{\pm }(z)\right).$$
(40)

so that each value of \({\Psi }_{ij}^{\pm }(z)\) at any vertex i is \({\phi }_{i}^{\pm }\). Now considering the FSI diagonal wave operator (18) as

$$\quad -{{\mathcal{C}}}_{{\mathcal{P}}}^{2}{\partial }_{z}^{2}{\Psi }_{ij}=-\left(\begin{array}{rc}{c}_{-}^{2}&0\\ 0&{c}_{+}^{2}\end{array}\right){\partial }_{z}^{2}{\Psi }_{ij}={\lambda }^{2}{\Psi }_{ij},$$
(41)

so that the 2-vector field \({\Psi }_{\lambda }^{ij}=\left(\begin{array}{r}{\Psi }_{ij}^{-}\\ {\Psi }_{ij}^{+}\end{array}\right)\) is the local eigenfunction of the diagonal 1D Laplacian (41), and the collection of continuous functions \({\Psi }_{\lambda }^{ij}\) along each metric graph edges denoted Ψλ is the metric graph Laplacian eigenfunction. It is important to stress here that every edge-specific component of the metric graph eigenfunction Ψij can equally be defined for a distinct, possibly edge-specific, velocity pair \({c}_{\pm }^{ij}\). Nevertheless, albeit, the theory can be developed within this more general framework, since precise, i.e., vessel-specific/patient-specific, data are not available, a simplified homogeneous wave velocity is chosen \({c}_{\pm }^{ij}\equiv {c}_{\pm }\) in the following (more details about the possible origin of velocities heterogeneity are given in SI). However, every result given below can easily be generalized for non-homogeneous wave velocities. The solution for the coupled FSI wave (13) can be found from an eigenmode decomposition in time domain19 using decomposition (19) to any edge the pressure/stress solution Pij inside each edge is decomposed over the base \({\Psi }_{ij}^{\pm }(z)\) such as

$${{\bf{P}}}_{ij}=\sum _{{\lambda }_{k}\in \Lambda }{a}_{{\lambda }_{k}}{\mathbf{\Pi }}{\Psi }_{{\lambda }_{k}}^{ij}{e}^{{\lambda }_{k}t},$$
(42)

where Λ is the spectrum of admissible eigenvalues the details of which will be later discussed. It is interesting to note that (42) is written, as (39) with the underlying hypothesis that the solid and fluid wave propagation parameters (i.e., wave velocities, solid to fluid density ratio, Poisson parameter) are uniform along the entire metric graph. This hypothesis is not needed for the subsequent approach to handle, but it significantly simplifies notations as well as computations. Since this hypothesis is reasonable for a broad range of contexts, it is kept to avoid unnecessary complexity. By construction since \({\Psi }_{ij}^{\pm }(z)\) define continuous fields along each edge of the metric graph, so does Ψij and Pij from linearity of (42). Hence, the continuity relations (38) are automatically satisfied (dimensionless version of relations (38) are the same without the ⋆ symbol). Furthermore, relations (34)–(37) — again, dimensionless version of those, being similar without ⋆ — can be rewritten using (42) in a general form at each vertex vi, as generalized Kirchoff–Robin conditions

$$-\sum _{j < i}{A}_{ij}^{S}{{\boldsymbol{\Pi }}}_{ij}^{{\prime} }{\boldsymbol{\Pi }}\frac{d}{dz}{\Psi }_{ij}({{\bf{x}}}_{{v}_{i}})+\sum _{j > i}{A}_{ij}^{S}{{\boldsymbol{\Pi }}}_{ij}^{{\prime} }{\boldsymbol{\Pi }}\frac{d}{dz}{\Psi }_{ij}({{\bf{x}}}_{{v}_{i}})={h}_{i}{\boldsymbol{\Pi }}{\Psi }_{ij}({{\bf{x}}}_{{v}_{i}}),$$
(43)

where we have defined \({A}_{ij}^{S}={S}_{ij}^{l}{A}_{ij}\) the non-homogeneous adjacency matrix taking into account a possible edge-dependent section \({S}_{ij}^{l}\), and also consider hi being a possible forcing chosen at vertex vi. The diagonal matrix

$${{\mathbf{\Pi }}}_{ij}^{{\prime} }=\left(\begin{array}{rc}1&0\\ 0&\frac{{S}_{ij}^{s}}{{S}_{ij}^{l}}\end{array}\right)$$
(44)

has also been introduced, where \({S}_{ij}^{s}\) & \({S}_{ij}^{l}\) are the (possibly edge dependent) liquid and solid section of eij ≡ ek’s edge. When hi = 0, (43) expresses a generalized Kirchoff condition for wave’s energy flux conservation at each bifurcation. hi ≠ 0 is used for external triggering nodes, i.e inlet node associated with the cardiac closure law or Windkessel outlet conditions as discussed in ref. 23. The inlet triggering nodes are where a time-dependent input flux is imposed into the network, thus modifying the FSI propagation within it. More precisely, they describe the heart’s coupling with aorta as a time-dependent periodic imposed flux, which, from the proportionality between the pressure-time derivative and the flux —through compliance, e.g., ref. 26— (also provided by (30)–(31)) leads to a forced pressure gradient condition. This time-dependent periodic flux is decomposed into Fourier modes and “coded” into the node’s flux boundary condition as a non-zero hi ≠ 0 parameter as detailed in ref. 23. A similar feature also apply for output triggering nodes, but instead imposing a time-dependent flux, a Windkessel law is applied herein. Eq. (43) can be re-written in order to simplify the homogeneous l.h.s (homogeneous refers to homogeneous Kirchoff’s condition at each vertex) as

$$-\sum _{j < i}{A}_{ij}^{S}\frac{d}{dz}{\Psi }_{ij}({{\bf{x}}}_{{v}_{i}})+\sum _{j > i}{A}_{ij}^{S}\frac{d}{dz}{\Psi }_{ij}({{\bf{x}}}_{{v}_{i}})=h_i {\boldsymbol \pi}_{ij}^{\prime-1} {\bf \Psi}_{ij}({{\bf x}_{v_i}}),$$
(45)

where the (non-diagonal) 2 × 2 matrix \({\boldsymbol{\pi }}^{{\prime}{-1} }_{ij}\) is

$${\mathbf{\pi }}^{{\prime}{-1} }_{ij}={\Pi }^{-1}{\boldsymbol{\Pi }}^{{\prime}{-1} }_{ij}{\boldsymbol{\Pi }},$$
(46)

being the diagonal-base change of diagonal matrix \({\boldsymbol{\Pi }}^{{\prime}{-1} }_{ij}\) which is

$${\boldsymbol{\Pi }}^{{\prime}{-1} }_{ij}=\left(\begin{array}{rc}1&0\\ 0&\frac{{S}_{ij}^{l}}{{S}_{ij}^{s}}\end{array}\right).$$
(47)

(46) shows that, when hi = 0, for generalized Kirchoff conditions, there is no supplementary coupling between the pressure and the stress from the conditions at the edges (34)–(37). Furthermore, even in the non-homogeneous case hi ≠ 0, if \({S}_{ij}^{l}/{S}_{ij}^{s}=1\) then \({\mathbf{\pi }}^{{\prime}{-1} }_{ij}={\mathbb{I}}\), no supplementary coupling arises since, the r.h.s of (45) keeps diagonal so that each \({\Psi }_{ij}^{\pm }(z)\) independently fulfill the boundary condition. Finally, defining the edge boundary condition heterogeneity parameter \(1+{s}_{ij}={S}_{ij}^{l}/{S}_{ij}^{s}\), one can evaluate

$${{\mathbf{\pi }}}_{ij}^{{\prime} -1}={\mathbb{I}}+{s}_{ij}{{\mathbf{\pi }}}^{{\prime} c},$$
(48)

with coupling matrix \({{\mathbf{\pi }}}^{{\prime} c}\)

$${{\mathbf{\pi }}}^{{\prime} c}=\frac{1}{{c}_{+}^{2}-{c}_{-}^{2}}\left(\begin{array}{rc}1-{c}_{-}^{2}&1-{c}_{-}^{2}\\ {c}_{+}^{2}-1&{c}_{+}^{2}-1\end{array}\right),$$
(49)

so that one finds a supplementary FSI coupling arising from edges conditions (34)–(37) only at the triggering sites where hi ≠ 0 and for non-homogeneous boundary condition parameter sij ≠ 0. Note that all the \({{\boldsymbol{\pi }}}^{{\prime} c}\) coefficients are all positive since c− < 1 and c+ > 1. From using (40) in (45) one can built the matrices \({{\mathfrak{A}}}^{\pm }\)

$${{\mathfrak{A}}}_{ij}^{\pm }(\lambda )=-{\delta }_{ij}\left(\sum _{m\in {\mathfrak{N}}(i)}\frac{{A}_{im}^{S}}{{c}_{\pm }}\cot (\frac{\lambda {\ell }_{im}}{{c}_{\pm }})+\frac{{h}_{i}}{\lambda }\right)+\frac{{A}_{ij}^{S}}{{c}_{\pm }\sin (\frac{\lambda {\ell }_{ij}}{{c}_{\pm }})},$$
(50)

so that relation (45) is expressed as a linear system upon vertex’s degree of freedom \({\phi }_{i}^{\pm }\)

$$\left[\left(\begin{array}{rc}{{\mathfrak{A}}}_{ij}^{-}&0\\ 0&{{\mathfrak{A}}}_{ij}^{+}\end{array}\right)-\frac{{h}_{i}}{\lambda }{s}_{ij}{{\mathbf{\pi }}}^{{\prime} c}\right]{\left(\begin{array}{r}{\phi }^{-}\\ {\phi }^{+}\end{array}\right)}_{i}\equiv {{\mathfrak{A}}}_{ij}{\left(\begin{array}{r}{\phi }^{-}\\ {\phi }^{+}\end{array}\right)}_{i}=0$$
(51)

(51) permits one to define the FSI secular matrix coefficients \({{\mathfrak{A}}}_{ij}\) as

$${{\mathfrak{A}}}_{ij}=\left(\begin{array}{rc}{{\mathfrak{A}}}_{ij}^{-}&0\\ 0&{{\mathfrak{A}}}_{ij}^{+}\end{array}\right)-\frac{{h}_{i}}{\lambda }{s}_{ij}{{\mathbf{\pi }}}^{{\prime} c}$$
(52)

Condition \(\det {\mathfrak{A}}(\lambda )=0\) sets the secular condition which permits to find the discrete spectrum of admissible λ. Following23 the spectrum Λ = ΛH ∪ Λp can be decomposed into homogeneous or internal modes (discrete, infinite, countable) set ΛH for which hi = 0 and external one Λp associated with a given triggering vertex (node) v0 where \({h}_{{v}_{O}}\ne 0\). For homogeneous modes, since \(\det {\mathfrak{A}}(\lambda )=\det {\mathfrak{{A}^{+}}}(\lambda )\det {\mathfrak{{A}^{-}}}(\lambda )\) the spectrum is the union of those of matrices \({{\mathfrak{A}}}^{\pm }\). This means that the homogeneous spectrum results from the independent superposition of each c± wave modes. The general quantum graph FSI solutions (42) can now be more precisely expressed as

$${{\bf{P}}}_{{\mathcal{G}}}=\sum _{{\lambda }_{k}\in {\Lambda }^{H}}\left({a}_{k}{e}^{i{\lambda }_{k}\tau }+{\overline{a}}_{k}{e}^{-i{\lambda }_{k}\tau }\right){\boldsymbol{\Pi }}{\Psi }_{{\lambda }_{k}}+{{\bf{P}}}_{p},$$
(53)

where \({\overline{a}}_{k}\) stands for the conjugated complex of ak, and Pp is the particular solution which decomposes into spectrum Λp which is the discrete Fourier decomposition of the finite-time triggering imposed at triggering node23. Here, the triggering node is located at cardiac valve location. Hence the total pressure/stress long-wavelength solution decomposes into a particular triggering dependent spectrum (given by the Fourier series decomposition of the heart closure valve) and a network related spectrum called homogeneous. Here we will not develop over more explicit expression for the particular solution Pp— which can be found as a direct extension of the results provided in23—, since it depends on the specific cardiac valve closing, which is not easy to measure and might be very variable among individuals. We will more closely consider the homogeneous solutions related to anatomical and physiological specificity of each individual’s arterial network.

Application to FSI waves within cerebral arteries

FSI-waves parameters within vascular networks

Using the c± parameters found in Section Method first discusses within matrix Π(17), needed in (53) and matrix (49) needed to build the secular matrix (52), one gets

$${{\mathbf{\pi }}}^{{\prime} c}=\left(\begin{array}{rc}0.04&0.04\\ 0.91&0.91\end{array}\right),\,\,\,\,\,\,{\mathbf{\Pi }}=\left(\begin{array}{rc}-1.4695&0.0666\\ 1&1\end{array}\right),$$
(54)

The last issue for explicit computations, is to set the outlet boundary conditions. Since the number of vessels involved at micro-circulatory scale is prohibitive, e.g., ref. 27, also considering that wave damping rapidly wipe-out its relevance at small scale, effective outlet boundary conditions are widely used in the literature, i.e., refs. 15,17 among others. Outlet conditions are then modeled by a Windkessel-like RCR boundary condition15,17 prescribing a local relation between the pressure p⋆ and the flux Q⋆ at the outlet

$$\frac{d{p}^{\star }}{d{t}^{\star }}+\frac{{p}^{\star }}{{C}^{\star }{R}_{2}^{\star }}=\frac{{R}_{1}^{\star }+{R}_{2}^{\star }}{{C}^{\star }{R}_{2}^{\star }}{Q}^{\star }+{R}_{1}{S}_{i}^{\star }\frac{d{Q}^{\star }}{d{t}^{\star }}$$
(55)

where, C⋆ is the micro-vascular compliance, \({R}_{1}^{\star }\) and \({R}_{2}^{\star }\) the proximal and distal resistances. In frequency domain, using the pressure-flux relation \((i\omega {\omega }_{0}^{\star }{\rho }_{f}^{\star }){Q}^{\star }=-{S}_{j}^{\star }d{p}^{\star }/d{x}^{\star }\), and introducing time relaxation parameter \({\tau }_{2}^{\star }=C{R}_{2}\), its dimensionless counterpart \({\tau }_{2}={\tau }_{2}^{\star }{\omega }_{0}^{\star }\), —\({\omega }_{0}^{\star }=2\pi {f}_{0}^{\star }\) — (55) reads as a complex Kirchoff–Robin relation applied at node j having section \({S}_{j}^{\star }\)

$$\frac{dp}{dx}=-i\frac{1+i{\tau }_{2}\omega }{1+i\frac{{R}_{1}^{\star }}{{R}_{1}^{\star }+{R}_{2}^{\star }}{\tau }_{2}\omega }\frac{K}{{\tau }_{2}{S}_{j}}p\equiv {h}_{j}(\omega )p,$$
(56)

where dimensionless parameter K given by

$$K=\frac{{L}^{\star }{\omega }_{0}^{\star }{\rho }_{f}^{\star }}{\pi {R}_{0}^{2\star }({R}_{1}^{\star }+{R}_{2}^{\star })}.$$
(57)

The parameters used for the dimensionless RBC condition (57) are taken from refs. 15 and 17 and matrix coefficients \({A}_{ij}^{S}\) used in secular matrix (43) are provided in SI.

Comparison with physiological spectrum

Figure 1A illustrates the anatomical context of Willis circle cerebral vascular network with the vessels numbering convention of ref. 28. Its vectorization using’s28 model is also shown in Fig. 1C, D where the vertex numbering is also provided. From these structural parameters, as well as physiological ones discussed in previous sections, one can evaluate quasi-analytically the wave propagation within the network from (53). The discussion and usefulness of such computation is postponed until next section, but, a specific practical application is illustrated here. The discrete spectrum obtained from the secularity matrix (52) zero determinant condition, i.e \(\det {\mathfrak{A}}(\lambda )=0\) permits to find the intrinsic vibrating modes of the FSI pressure/stress wave systems, which does not depend on the heart beating pulse. Given matrix \({\mathfrak{A}}\), and computing its determinant, this result in solving a transcendental equation numerically to find the first discrete element λH in the infinite set ΛH (more details are given in SI). Mathematically, these λH eigenvalues are invariant quantities reflecting topological (linked with graph adjacency matrix), metric (related to vessels lengths and radius), and physiological (wave velocities depend on arterial walls elastic properties) properties of the arterial network. These modes, albeit externally triggered by cardiac pulse, have intrinsic frequencies, as found in ref. 23. The external triggering only affects their amplitudes23. Hence, one should be able to decipher, from a high-frequency pressure recording the signal signature of these modes, so that the arterial networks pulse response could be specifically followed-up, retrieving the influence of the cardiac input.

Figure 3 illustrates this concept from the Fourier-spectrum recording of the Intra-Cerebral Pressure (ICP). Figure 3A provides one Fourier spectrum measured in ref. 29, into which the cardiac fundamental pulsation, and its subsequent 12th Fourier harmonics are represented (green dots). Figure 3B shows the signal’s cardiac signature in green, modeled with Gaussian-shaped peaks. Filtering these cardiac contributions—as typically done in waveform analysis30— to the original signal permits to more clearly visualize the non-cardiac features of the spectral content, illustrated in Fig. 3C & D. On the top of this retrieved spectrum, Fig. 3C, D illustrate with dotted vertical lines, the secular FSI frequencies predictions \({f}_{k}^{\star }={\lambda }_{k}{f}_{0}^{\star }\) for \({f}_{k}^{\star } < 14Hz\) and λk ∈ ΛH for two choices of outlet boundary conditions parameters. In both cases, some of the predicted FSI frequencies match quite well with measured ones (black arrows). This matching between prediction and observations has been obtained without any parameter fitting (all parameters are given in SI i.e., matrix AS, L and Table S1. and taken from the literature). Surprisingly, the first predicted eigenvalues λk ∈ ΛH, with k ∈ [1, 10] very weakly depend on the applied outlet boundary conditions, as can be qualitatively observed in Fig. 3C, D (quantitatively confirmed in SI Table S2). More precisely a sensitivity analysis on the four Windkessel RCR boundary conditions parameters changing their value by ± 10% for both models of refs. 15,17 (whom parameters are given in SI) reached a maximum relative difference of 10−8 over the first twenty eigenvalues.

Fig. 3: Retrieving cardiac’s signature out of cerebral pressure Fourier Spectrum shows intrinsinc vascular frequencies predicted by our approach.
figure 3

A Spectrum of Cerebral pressure signal of a patient from CENTER-TBI study from ref. 29. The fundamental frequency f0 = 1.06 Hz of the cardiac pulse and its various Fourier harmonics fn = (n + 1)f0 (n ∈ [1,11]) are represented with Green dots. B Same as A where Gaussian profile centered in each f0 frequencies (with standard deviation σ = 0.07) for each Cardiac mode has been represented in green. C The black signal is A’s spectrum filtered from its cardiac contribution represented in green in B. Dotted blue lines are the discrete spectrum found for the pressure pulse propagation within Willis circle using the RBC boundary condition parameters of ref. 17. D. Same conventions as C with Dotted red lines RBC boundary condition parameters of ref. 15.

Full size image

Discussion

Cardiac-filtered ICP signal has been chosen to illustrate how FSI-wave signature can be used to decipher arterial network intrinsic features. Following-up these specific features during patient’s life could result in intrinsic frequency shifts and/or intrinsic amplitude variations providing direct information about the circle of Willis response. Such clinical investigations are indeed possible using non-invasive high-frequency ICP measurements developed over the last decade—e.g., optic nerve sonography31, trans-cranial Doppler30—. Since, the circle of Willis structure and mechanical properties might play a role in the development of cerebral aneurysms32, as well as other diseases this illustration already reaches a significant clinical interest. But this is only one possible application. The proposed theory and modeling has been developed for any arterial network, and can be applied to the entire macro-vascular arterial circuit. Furthermore, more elaborated use of the proposed theory could be developed in the future. Patient-specific arterial networks can be obtained much more precisely from anatomical vessel skeletonization using medical images — e.g., computed tomography scans, magnetic resonance imaging —. Given the resulting patient-specific metric graphs, the FSI wave predictions could be matched with observations using data assimilation methods to estimate patient-specific FSI-wave parameters (e.g., vessel specific wave velocities \({c}_{\pm }^{ij}\)). This is an interesting avenue to build physically-informed parameter estimation rather than statistical indexes from clinical measurements. Hence, the hereby proposed theory can lead to real-time patient-specific monitoring of pulse wave propagation, long-term follow-up of individual arterial response, patient-specific data assimilation strategies, which could permit to anticipate and/or detect some diseases.

Let us now comment about how visco-elastic response of arterial walls might be taken into account in the proposed framework, an important issue, not yet discussed. As known in single vessel configuration6,20,33, as well as in networks,13 visco-elastic effects are affecting the pressure wave’s damping. However, even if the visco-elastic response of arterial walls both affects the pressure signal and its damping18 they neither modify the elastic mode’s peak20 (justifying Fig. 3 analysis which ignores visco-elastic effect influence) nor the coupled linear FSI wave system (13)–(14). The additional complexity provided by visco-elastic effects can be re-casted into dispersive wave-speeds—i.e frequency dependent wave speed c±(ω)— of the FSI system which can even be explicitly deduced from vessel’s rheological properties18,20. Furthermore, these dispersive wave-speeds can indeed be taken into account in the provided formalism, leading to the very same explicit solution in frequency domain. Obviously, the modeling of viscoelastic effects need more parameters to be estimated and/or known. This is where the discussion of the previous paragraph is even more relevant, for future investigations and applications. Finally, it is also interesting to mention that, although generally overlooked1,2,3,34 the influence of blood’s viscoelasticity has recently been numerically investigated in ref. 35, where it is shown that both pressure and velocity amplitude of the 1D blood hammer wave are weakly affected by non-Newtonian effects, as opposed to shear-stress and near wall profiles. To our knowledge, the influence of non-Newtonian effects on the spectral content of the pressure signal has not been investigated yet.

Conclusion

A new theory for FSI-wave propagation within arterial networks has been proposed using long-wavelength approximation to establish continuity conditions to be applied at bifurcations. From solving the coupled FSI-wave problem onto the velocity matrix diagonal base, coupling conditions at network’s nodes have been derived, providing the quantum graph boundary conditions to be satisfied at each vertex. The wave propagation solution is then decomposed onto the quantum graph’s eigenfunctions base. Vertex’s boundary conditions permits to build the secular matrix associated with FSI-wave’s propagation, from which the wave spectrum’s condition is deduced. Explicit quasi-analytical solutions for the pressure/stress have been obtained, for any network (into a compact metric graph). These solutions display a homogeneous intrinsic part related to the arterial network properties, and an external one, triggered by cardiac pulse, generalizing23 when including FSI effects. One illustration of the theory has been provided for the circle of Willis network. The predicted spectrum has been confronted with real cardiac filtered signal, showing consistent intrinsic arterial pulsating modes inside the recorded spectrum. This work open new perspectives for real-time patient-specific monitoring of pulse wave propagation, long-term follow-up of individual arterial response and patient-specific data assimilation strategies.