In silico coronary wave intensity analysis: application of an integrated one-dimensional and poromechanical model of cardiac perfusion.

Coronary wave intensity analysis (cWIA) is a diagnostic technique based on invasive measurement of coronary pressure and velocity waveforms. The theory of WIA allows the forward- and backward-propagating coronary waves to be separated and attributed to their origin and timing, thus serving as a sensitive and specific cardiac functional indicator. In recent years, an increasing number of clinical studies have begun to establish associations between changes in specific waves and various diseases of myocardium and perfusion. These studies are, however, currently confined to a trial-and-error approach and are subject to technological limitations which may confound accurate interpretations. In this work, we have developed a biophysically based cardiac perfusion model which incorporates full ventricular-aortic-coronary coupling. This was achieved by integrating our previous work on one-dimensional modelling of vascular flow and poroelastic perfusion within an active myocardial mechanics framework. Extensive parameterisation was performed, yielding a close agreement with physiological levels of global coronary and myocardial function as well as experimentally observed cumulative wave intensity magnitudes. Results indicate a strong dependence of the backward suction wave on QRS duration and vascular resistance, the forward pushing wave on the rate of myocyte tension development, and the late forward pushing wave on the aortic valve dynamics. These findings are not only consistent with experimental observations, but offer a greater specificity to the wave-originating mechanisms, thus demonstrating the value of the integrated model as a tool for clinical investigation.

Introduction

Impaired myocardial blood flow underlies a wide range of cardiac diseases. Notably, there are 3M coronary angiographies undertaken every year to diagnose coronary artery disease in Europe alone (Cook 2011). The ischaemic cascade, which begins with a perfusion deficit and progresses through mechanical and electrical dysfunction to ultimately culminate in myocardial infarction, has been well documented (Nesto and Kowalchuk 1987). In addition, myocardial ischaemia has received attention as a secondary contributor to several other cardiac diseases including hypertrophic cardiomyopathy (Maron et al. 2009), aortic stenosis (Rajappan et al. 2003) and heart failure (Nakanishi et al. 2012). Perfusion deficits can be brought on by disparate causes, including stenotic lesions in large vessels or microvascular dysfunction, as well as the increase in extravascular compression which is known to heighten subendocardial vulnerability to ischaemia (Hittinger et al. 1995; Heusch 2008). Due to this physiological complexity, the differential diagnosis and the therapeutic trajectories must be pieced together from additional clinical evidences. Accordingly, there has been much interest in a diagnostic approach which would unite the currently overlooked assessments of coronary microcirculation and perfusion–contraction crosstalk into the existing routine clinical workflow.

One technique which has demonstrated potential as an integrated measure is wave intensity analysis (WIA) (Parker and Jones 1990; Parker 2009). The theoretical development of this approach was based on a linearised one-dimensional approximation of wave propagation in elastic arteries, and the characteristic analysis of the resulting hyperbolic system. From simultaneous pressure and velocity waveforms measured via coronary catheterisation, this technique enables the assessment of the combined myocardial and haemodynamic functions within a single diagnosis (Davies et al. 2011; Kyriacou et al. 2012). Recently, a landmark predictive application of WIA was demonstrated for the first time in de Silva et al. (2013), which showed that a real-time WIA-derived index can be used to predict the functional recovery following myocardial infarction. Nonetheless, this study, as well as all other coronary WIA investigations, is currently encumbered by the lack of a method with which to relate WIA observations to the mechanistic foundations underlying the pathophysiology. Without the ability to attribute the altered coronary flow and myocardial dynamics to the chain of events that generate the coronary waveforms, the reverse process of mapping these signals to a specific underlying disease process remains in the domain of qualitative trial and error.

As an alternative to this approach, the convergence of medical imaging with computational modelling allows subject-specific coronary WIA to be studied from a causal mechanistic basis. This is demonstrated in the current work using an animal subject. We require such a model to consist of coronary flow and myocardial mechanistic components coupled together, as the coronary waves are produced by their interaction. To tractably model the coronary circulation necessitates a multiscale strategy, whereby the inertia-dominated flow in large conducting vessels and the predominantly viscous flow in the distal circulation employ alternative modelling paradigms (Lee and Smith 2012). The one-dimensional vascular flow modelling framework has been well established in the literature in the past decade and has been shown to reproduce wave propagation behaviour with high accuracy (Vosse and Stergiopulos 2011). In this work, we combine and extend our previous studies in coupled one-dimensional flow–myocardium (Smith et al. 2002) and poroelastic perfusion modelling (Cookson et al. 2012).

The theory of poromechanics provides an apt framework for perfusion modelling, as it forgoes the need for explicit microvascular geometries and innately addresses crosstalk effects with spatial variability. In the standard poroelasticity, the macro-scale elasticity and fluid conservation laws are derived based on general assumptions without knowledge of the detailed microstructure (Bowen 1980). This is in contrast to the alternative approach of mathematical homogenisation which rigorously determines the governing macroscopic laws and effective constants directly from a specific micro-scale geometry (Cioranescu and Donato 1999). However, although several studies have employed homogenisation to vascularised tissues (Chapman et al. 2008; Shipley and Chapman 2010; Rohan and Cimrman 2010) including those in the large-deformation regime (Rohan 2006; Rohan and Lukes 2011), the requirements of microstructural periodicity and the necessity of numerical solution of the cell problem at every Gauss point currently limit the utility of this approach in large domains (Rohan and Lukes 2011). Accordingly, studies aiming to investigate the perfusion dynamics of the whole heart (Rossi 2007; Chapelle et al. 2010; Cookson et al. 2012) have applied large deformation poroelasticity Biot (1972); Coussy 1995, 2004); Boer 2005). It should be mentioned also that although the theory deals exclusively with macroscopic quantities, the micro–macro averaging technique (Whitaker 1986) has retrospectively provided theoretical rigour to the fundamental governing laws and remains an indispensable tool for characterising the macroscopic properties of a specific physical medium with known microstructure.

It is noteworthy that, to date, the most advanced application of the porous perfusion models has yet to achieve a clinically relevant milestone. Key remaining challenges include the development and validation of a suitable cardiac poroelastic constitutive law, characterisation of hydraulic permeability tensors through coronary microvascular analysis, and an effective strategy for posing boundary conditions on the porous flow domain. While the former two tasks represent an ongoing programme of research beyond the current scope, in this article, we will approach the boundary condition issue in a data-driven manner, i.e. through explicit image-based modelling of the flow in proximal vessels, thus capturing the complex spatial distributions of the feeding sites. Consequently, the presence of the major coronary vessels in the model allows WIA to be performed in silico, which will in turn benefit from the increased realism of the distal flow conditions and perfusion–contraction coupling.

The principal objective of this work can thus be summarised as the establishment of an integrated multiscale computational model of myocardial blood flow that incorporates the effects of perfusion–contraction crosstalk. Its utility as a physiological research platform is illustrated through an in silico application of wave intensity analysis under various perturbed conditions, thus allowing the direct wave dependence on individual myocardial and haemodynamic parameters to be established.

Materials and methods

The integrated computational model described herein consists of four major components (see Fig. 1), encompassing passive and active myocardial mechanics, coronary flow and myocardial perfusion.

Fig. 1

A schematic overview of the integrated perfusion model

Poromechanics

Porous material In this work, cardiac tissue is idealised as a saturated porous medium consisting of a solid phase (myocardial tissue) and a fluid phase (blood). To proceed, the microstructural particularities are disregarded, and the two constituents are conceptualised to exist co-spatially with varying volume fractions, while retaining their individual properties. Thus, in the following discussions, we refer to the macroscopic solid continuum (referred to as the ‘skeleton’) and the fluid continuum. This concept had been firmly consolidated by the time the classical work of Biot (1941) and mixture theory Truesdell (1957); Truesdell and Noll 1965) arrived, and allows us to consider a solely macroscopic framework.

Denoting a representative elementary volume (REV) as , we can define the respective initial and current Eulerian porosities as the fluid volume fraction, i.e.where the subscript o and superscript f correspond to reference configuration and fluid, respectively. For convenience, we also denote the solid fraction as . Note that in describing the deformation of the total porous medium, it is convenient to exploit the skeleton deformation since it is directly observable (Coussy 1995). Therefore, the term Lagrangian implies with respect to the skeleton.1 Following the standard definitions of large deformation kinematics, the material initially located at X in the reference configuration is identified with the deformed position x, uniquely mapped bywhere U denotes the displacement vector. The deformation gradient tensor is defined asleading to the right Cauchy–Green strain tensor and Jacobian . The fluid flow through the porous medium can thus be described in terms of the Darcy velocity ,which is the velocity of fluid relative with respect to that of the skeleton weighted by the porosity. Darcy velocity can be expressed in the Lagrangian frame by introducing the following identity relating the flow through an infinitesimal oriented area element which, in combination with Nanson’s formula, gives Particle and material derivatives To examine the conservation of physical quantities, two types of derivatives particular to the porous medium can be defined. We focus on the Lagrangian formulation in the following. A particle derivative with respect to a solid or a fluid particle is the time derivative that an observer attached to the particle would derive (Coussy 2004). Let be the Eulerian volume density of a physical quantity in the current configuration. We can writefor the corresponding Lagrangian particle density . Now denoting the volume integral byit follows from Eq. (7), the particle derivative is expressedThe material derivative accounts for the incoming flux of the considered quantity carried by the fluid, additionally to the solid particle derivative. Hence,where denotes the Lagrangian density function of the fluid particle, which is located by the skeleton position vector . By applying (9) and divergence formula to the above, we obtainEquations (9) and (11) state that the particle derivative can be considered as the time derivative of the function , but the same cannot be assumed of the material derivative in general since the skeleton and the fluid may undergo different motions.

Conservation of mass The conservation of fluid mass can be considered by posing (7) in terms of mass density where m is defined as the current additional fluid mass content per unit reference volume , above the initial level. If a distributed source is present,Applying (11) to the above leads to the Lagrangian fluid continuity equationIn the absence of solid phase mass creation, the balance of the skeleton mass is simply Equation of motion The conservation of momentum in the porous medium is often examined for the whole mixture, since the balance within an individual phase gives rise to a distributed interaction force term which must be additionally defined. For brevity, we refer the reader to a derivation elsewhere (Coussy 1995) and focus instead on describing the results. The Lagrangian equation of motion readswhere and represent the second Piola–Kirchoff stress, body force density and acceleration, respectively. The two density-like quantities and denote, respectively, the skeletal and fluid mass content per unit reference volume , such thatnoting that remains constant and equal to the reference by virtue of (15).

Stress partition Consistent with the standard theory, the Cauchy stress is related to and first Piola–Kirchoff stress viawhere the total stress in the macroscopic medium is the sum of individual phase stresses , which in turn, is posited to be the volume-weighted average of the actual stress in the microscopic constituent. Consequent investigations with micro–macro averaging technique has confirmed the correspondence between the integral of microscopic free energy potential with the macroscopic stresses (Buhan et al. 1998), and derived macroscopic equation of motion from the microscopic momentum balance (Coussy 2004).

The fluid stress tensor can be expressed as , where p is termed pore pressure and represents the fluid shear stress. In practice, the shear term is often disregarded. In the coronary circulation, we assume that such a contribution is secondary to the drag at the internal walls (the main mechanism behind the Darcy flow; see below), and adopt in the following, including also for the one-dimensional flow formulation. Lastly, we remark that this hydrostatic stress state inherently embeds a constitutive assumption, and a specific formulation should be posed consistently.

Darcy’s law The pore fluid motion is modelled by Darcy’s law, which is a linear conduction law in the formNote that here the effects of gravity, fluid inertia or shear (Brinkman correction) are ignored, since in the coronary perfusion, we regard such terms to be of secondary importance. In general, the permeability K is a tensor. To ensure that a pressure gradient accelerates the flow in a consistent direction, K must be positive definite. The phenomenological foundations underlying Darcy’s law has been consolidated through subsequent micro–macro averaging derivation, which arrived at the identical expression by assuming Stoke’s law and fluid incompressibility in a general porous medium without discontinuity (Whitaker 1986).

As the pores deform, the permeability of the medium will also undergo changes. By considering the Poiseuille conductance C of a representative cylindrical vessel, we writewhere l, r, a and v correspond to the length, radius, area and volume of the segment. Upon changes to the vessel volume,If we assume that the change in segment length is negligible, (21) becomesIf, on the other hand, isotropic permeability is assumed, in three dimensions such thatEquations (22) and (23) represent the upper and lower bounds on scaling.

Constitutive law

This section details the development of a specific constitutive law applied in numerical simulations (Sect. 3). The formulation is intended for cardiac tissue, for which both the solid and fluid phases can reasonably be treated as incompressible.

Thermodynamic perspectives on constitutive law In a continuum, the second law of thermodynamics takes the form of Clausius–Duhem inequality. The general expression for a porous medium has been previously derived (Coussy 1989; Dormieux and Stolz 1992). In the isothermal regime, the intrinsic dissipation D iswhere is the (Lagrangian) free energy density of the skeleton. Note that represents the strain work performed by action of the pore pressure on the skeleton via the internal walls. It is observed that in the absence of this term, we would recover the dissipation of a standard solid involving only the strain work rate . The term is known as the Lagrangian porosity, which measures the current fluid content per unit reference volume. While on heuristic grounds, one may expect to be associated with the fluid action, because fluid inflow will, in general, lead to an associated change in offers a more appropriate account of the fluid work rate. Thus (24) shows that p is the thermodynamic force driving the change in .

Poroelasticity is characterised by zero intrinsic dissipation in the system. Hence, putting leads to the state equationsNote that inherent in (25) is an assumption of normality between the state variables E and , that is, variation of either particular state can occur independently of the other variable. This property is demonstrated by materials with compressible microscopic solid constituent, whereby the pore volume can change in the absence of macroscopic deformation.

Matrix incompressibility

The altered state of the skeletal stress under solid matrix incompressibility can be examined by adapting the dissipation expression appropriately. By putting , Eq. (15) can be written assuch that . Now using the volume change rate identity (Bonet and Wood 2008, Chpt4)Equation (24) can be re-expressed asThe new expression shows that the physically pertinent stress in the medium is no longer the total stress S, but rather the component excluding pore pressurewhich is classically referred to as the effective stress. It also shows that is no longer a state variable, thus reducing (25) to simplyIn other words, the knowledge of p or is no longer required to determine the skeletal free energy, as the actual work performed by the fluid on the skeleton via internal walls will be self-evident from the observed deformation F by virtue of the constraint (26). Note however, in contrast to the hyperelastic incompressibility, the new constitutive equation (30) does not have to be posed in terms of the strictly distortional component of the strain tensor. In fact, it would be undesirable to do so in the current context, as it would imply zero resistance against volumetric dilation. The difference arises from the fact that matrix incompressibility does not imply skeletal incompressibility—through fluid inflow, macroscopic volume change is still possible, thus the macroscopic medium must be regarded as compressible.

Method of Lagrange multiplierAlthough the definition of effective stress reveals physical insights, (30) leaves pore pressure indeterminate, and the solution process must devise a way to calculate consistent fluid stress states. For materials with nearly incompressible matrix, a procedure described in Chapelle and Moireau (2014) may be employed. In this work, the incompressibility of the material is addressed by the method of Lagrange multiplier. First, the general free energy expression admitting matrix compressibility effects is augmented with the constraint (26)where the superscript denotes a constrained quantity. It follows from (25) that where (33) was expanded using (32). Upon comparison with (29), the effective stress is identified aswhich can be shown to be consistent with (30) if we consider as a function of E such that , and the equivalence in rates of J and due to (26) upon fulfilment of the constraint.

Compaction limitFor a physically sensible behaviour, the constitutive law must appropriately address the limit cases, which occur when the porosity reaches a value of 0 or 1. As pointed out in Chapelle et al. (2010), Eq. (26) guarantees for ; therefore, no explicit measures are required to ensure . Against the compaction limit , a barrier potential previously proposed by Federico and Grillo (2012) is considered here. While the purpose of the original formulation was to control the bulk modulus in the solid skeleton, unnecessary in the current work due to the Lagrange multiplier approach, it has several desirable properties which are exploited here to modify the pore pressure characteristics. In particular, the adapted potential remains inactive until compaction is approached (in that it contributes zero pressure and zero elastance at non-negative volumetric strain), but when active, tends towards infinite pressure and elastance, which can be used to prevent further fluid extraction. These conditions can be expressed as and the barrier function readswhere denotes the Heaviside step function which ensures that (38) becomes active only when , where . The parameter r is selected within (0, 1], and q is a positive integer. Given these choices, (38) will satisfy (35)–(37).

Free energy compositionBased on the above developments, we now compose a specific constitutive law. The general form of the free energy has been built up to this point aswhere fully characterises the behaviour of the porous material above the compaction limit. To facilitate parameter decoupling, we further propose a splitting of this potential such thatSuch a decomposition permits the direct incorporation of previous (separate) models of cardiac constitutive and coronary pressure–volume relations.

On the other hand, the response of the macroscopic skeleton depends on the specific composition of its constituent phases. In the simplest case, this can be represented by the initial volume fractions and in the free energy expression, to capture correct proportional contribution of the skeleton to the total energy density. Considering the extreme cases, when , we would expect to vanish. Likewise, will recover a purely solid material. In between these bounds, it can be expected that the lower the initial porosity, the greater the pressure required will be in order to raise the fluid content of the medium by a given amount. A simple modification to (40) can satisfy these conditionsDue to the free energy being a linear summation, these new scaling terms can be absorbed into the constitutive parameters of the form (40). However, (41) serves to demonstrate that these parameters reflect the volume fractions of the phases at the reference state.

For the specific form of , an existing cardiac constitutive law can be employed. Here, we select the structurally based law of Holzapfel and Ogden (2009), which encompasses the preceding orthotropic-type laws (and their various simplifications). For , we adopt the experimentally characterised pressure–volume relation of Bruinsma et al. (1988),leading towhere the constant exists to satisfy the condition at rest volume . In (43), the exponential term dominates during large net mass inflow, while the log term dominates during small or negative mass inflow.

Active stressFor active stress generation we adopt a modified form of a previously proposed model (Kerckhoffs et al. 2003) of the formwhere and . Here, and denote the peak stress scaling parameter, activation duration and fibre stretch ratio, respectively. is assigned as an additional component in the total stress tensor.

One-dimensional vascular flow

Network flow model

The one-dimensional model of vascular flow has been well established in the literature (see Lee and Smith (2012) and Vosse and Stergiopulos (2011) for review). We follow the formulation employed in previous work (Lee and Smith 2008) with conservation equationswhere A and Q denote the cross-sectional average area and flow rate, K represents the friction coefficient at the lumen wall, which can be written as (Smith et al. 2002)The parameter controls the shape of the flow profile across the cross section of the vessel.

The derivation of the above system involves posing a pressure–area relationwhich leads to the distensibility D and wave speed c showing that vascular wave speed is proportional to the square root of wall stiffness parameter . The junction equations include the conservation of mass and total pressure and the compatibility relations as presented in Lee and Smith (2008).

Coupling with porous domain

The geometrical interface between the explicit vascular and porous flow regimes within the modelling framework is represented by the meso-scale vessels which progressively bifurcate into smaller segments forming a distributed network. A simplified treatment of this transition zone as a point source in the porous domain is undesirable, since localised outflow will lead to the development of unphysiological pressure and velocity peaks. Therefore, we assume that the exchange of fluid between a terminal vessel and the porous tissue occurs within a volume surrounding the distal end of the vessel , such thatwhere, for clarity, the variables associated with the one-dimensional vascular domain are denoted by the subscript 1D and S denotes a distributed source in the porous domain. Furthermore, we express S via a distribution function f which, together with (51) impliesIn the absence of detailed anatomical information, we approximate f with a Gaussian function. The pressure–flow relationship of the coupling interface can be established by regarding the meso-scale vessels as a predominantly resistive element such thatwhere, as a first approximation, we define the average pressure as

Systemic haemodynamics

As is common in the literature, we employ a windkessel-type systemic boundary conditions to the ventricular model. The aortic valve dynamics has been previously characterised by Korakianitis and Shi (2006), based on an orifice model accounting for leaflet angular position. However, the original formulation allowed instantaneous changes in the aortic valve flow, which resulted in rapid fluctuations of pressure following the isovolumic phases. To address this issue, the valve flow dynamics has been altered to be a first-order ODE such that where and denote the valve angular position, aortic pressure and LV pressure, respectively. The systemic boundary condition is modelled with a three-element windkessel, which had been shown to accurately reproduce physiological systemic impedance (Westerhof et al. 2009) where represents the lumped systemic pressure. The flow through the mitral valve is modelled asThe coupling between this system of equations and the ventricular model is achieved through the cavity pressure , which is induced by the myocardial contraction and set in (59), and which is updated by (56) and imposed on the deformation through application of the following constraint to the finite element problem using a Lagrange multiplierwhere denotes the endocardial surface.

Solution procedures

The numerical solution of the coupled model is accomplished using CHeart, an mpi-based multi-physics finite element solver developed at King’s College London (Lee et al. 2016). Due to the disparate parallel computational requirements of each subproblem, we employ a sequential approach whereby the solution of each coupled system is iterated until a nonlinear system convergence is achieved at each time step. The poroelastic system (14), (16) and (19) was discretised using a Galerkin finite element method with quadratic/linear/linear elements for the displacement/pressure/mass mixed formulation, while the one-dimensional vascular system (47) was discretised using a fifth-order spectral element method with Crank–Nicolson time stepping. Details of mesh construction and refinement are presented in Sect. 3.1. The simulations were executed on the in-house HPC resource (SGI Altix UV 1000). With 128 cores per simulation using a time step size of , the typical solution wall time was around 25 h.

Further details of the numerical formulation of individual physics can be found elsewhere (coronary flow (Lee and Smith 2008), elasticity (Hadjicharalambous et al. 2014) and porous flow in large deformation (Cookson et al. 2012)).

Results

Geometrical model construction

The image-to-model generation process is briefly outlined here and summarised in Fig. 2. The heart and coronary geometries were obtained from a previous experimental study performed on a 40-kg Danish Landrace pig (Schuster et al. 2010). After injection with a fluorescent vascular casting polymer in the left main stem, the heart was frozen in an unloaded state and imaged using a cryomicrotome multichannel acquisition at m resolution (Spaan et al. 2005). Subsequently, the myocardial geometry was segmented, and a custom template mesh fitting procedure was employed. The right ventricular wall and papillary muscles were excluded, as is customary in clinical perfusion analysis.

Fig. 2

The computational simulation meshes were created from 50-m-resolution 3D image stacks of porcine myocardium and coronary vasculature. A maximum intensity projection (MIP) of the vascular images is shown in a. The myocardial mesh was truncated near the valve plane, and vascular network was modified to achieve an even distribution of terminal segments throughout the myocardium (see text for details). The distribution of lengths from vascular root to each distal terminal node  cm) is shown in c

The vascular images were reconstructed as described in Goyal et al. (2013) to full detail. The resulting mesh was then truncated to restrict its extent to the upper arterial network, by dividing the myocardium into approximately 0.02 mL volumes and seeking a supplying vessel segment within each, with a target diameter of m. These seeded segments were traced back up to the root and the set of all traversed segments comprised the simulation domain. This ensured an even distribution of vascular termini throughout the myocardium. Due to there being no polymer injected in the right coronary vessel, the inferoseptal and parts of the anteroseptal wall lack vasculature. In addition, vessel segments were assumed to be of uniform unstressed area, as the resolution of the imaging system did not permit a reliable characterisation of the tapering behaviour in all segments. The resulting truncated mesh has an LV volume of 40.6 mL, unloaded LV cavity volume 48.5 mL, and 3910 vascular segments with 1990 distal termini totalling 2.4 mL excluding the microcirculation. The distribution of path lengths from the root to each terminal ( cm) is shown in Fig. 2c.

The computational mesh for the poroelastic problem was generated using quadratic/linear hexahedral elements, yielding around 30,000 degrees of freedom. The vascular mesh was refined to yield around 187,000 degrees of freedom, which was found to be necessary to resolve the accurate wave propagation behaviour as indicated by a mesh refinement study.

Physiological baseline conditions

Parameters The physiological baseline parameters were largely adopted from the literature. The passive myocardial parameters were obtained from Holzapfel and Ogden (2009), for which a downscaling was necessary to tune the diastolic behaviour. A uniform factor of 0.25 was applied to the parameter set. In the absence of subject-specific fibre measurement, a linear endo-to-epi distribution of fibre angles from to was used (Zhang et al. 2013). For the active tension model, adjustments were made to and to achieve a sufficiently fast relaxation rate. For clarity in coronary wave interpretation, a simple linear activation wave spreading from endocardium to epicardium lasting 60 was prescribed without a base–apex gradient.

The permeability of the porous domain was assumed to be isotropic. Although no porcine-specific data was available, a recent study on the capillary network of rat myocardium has estimated the principal component of the permeability tensor to be on the order of (Smith et al. 2014). Using the ratio between reported capillary densities in pig (Koudstaal et al. 2013) and rat (Kerckhoven et al. 2004) myocardium, the permeability was adjusted to be . The porous constitutive parameters were based on the capillary compartment parameters in Bruinsma et al. (1988), but , which scales the volume response, was adjusted to include a partial influence from the larger, arteriolar microvessels. As only a single porous compartment was used in this model, the role of the porous domain should be considered to represent the combined behaviour of both arteriolar and capillary microcirculation.

The elastic properties of the vessel segments were determined by assuming a constant uniform wave speed of , thus implying a constant distensibility. Currently, there is no gold standard technique for measuring local coronary pulse wave velocity, and inter-subject variations ranging from (Aguado-Sierra et al. 2006) to above (Rivolo et al. 2014) have been reported in pigs. All other parameters were determined by manual tuning. The complete set of baseline parameters are listed in Table 1.

Table 1

Baseline parameters for the systemic, haemodynamic and contractile models

Systemic circulation parameters			Coronary haemodynamic parameters			Active stress parameters
θ_max	75	^∘	ρ	1.05×10^-6	kgmm^-3	a_4	3.2	s
T_ao	0.01	s	μ	3.36×10^-6	kPa s	a_6	2.0	Dimensionless
CQ_ao	9.59×10^5	mm^3/(skPa^0.5)	α	1.05	Dimensionless	a_7	0.7	Dimensionless
K_1ao	750	(kPas)^-1	c	15,000	mms^-1	T_0	15,000	kPa
K_fao	50	s^-1	q_1	0.022	kPa	t_r0	0.16	s
R_a	1.2×10^-5	kPas/mm^3	q_2	1.009	kPa	t_d	0.03	s
R_s	2.2×10^-4	kPas/mm^3	q_3	80	Dimensionless	t_max	0.3	s
C	7000	mm^3/kPa	ϕ_o	0.06	Dimensionless	QRSd	0.06	s
CQ_mi	2.0×10^5	mm^3/s	γ	1.0×10^-7	mm^3s^-1kPa
C_mi	37.5	(kPa)^-1	p_drain	0.5	kPa

Boundary conditions To prevent free body translation of the ventricle, the epicardial ring of nodes at the base was fixed in the longitudinal direction; however, the base plane surface was permitted to deform under a Maxwell-type constraint. This was necessary to prevent large unphysiological hydrostatic pressures developing in the basal region which could adversely affect the porous flow.

The lateral component of translation was restricted by a linear-spring-type model applied to the epicardium which penalised displacement in the surface-normal direction.

The pressure at the inlet of the vascular network was prescribed to be equal to the aortic sinus pressure determined by the systemic windkessel model. The outflow conditions from the porous domain was prescribed as a variable resistor of the formwhere was subject to the same scaling as the permeability as in (23).

Baseline function The results of the baseline simulation are shown in Fig. 3. In addition, the bull’s eye diagrams in Fig. 3e, f depict segmental perfusion (estimated by the sum of coronary outflow from the vascular mesh) and fluid mass accumulation (m) in the mid-equatorial ring with three transmural layers. The global functional measures are as follows: ejection fraction 49.5 %, perfusion 6.2, perfusion pressure 129/69, diastolic time fraction 63 %.

Fig. 3

Baseline results a coronary inlet pressure, LV cavity pressure and distal pore pressure (median and inter-quartile range over the nodes). Although not shown in this figure, the maximum pore pressure in excess of LV pressure was found during systole, consistent with experimental observations. b Inflow and total outflow of the upper arterial vascular network. Inflow exceeds the outflow in early systole as the influence of cardiac contraction is greater on the distal arterial flow. The reverse happens during late systole, as the stored flow is discharged. c LV cavity volume and aortic outflow. Transient reversal in flow is observed at end systole, enabled by the modified valve dynamics (see text for details). d Coronary flow across transmural layers shows an augmented systolic flow in the subepicardial layer and a reversed flow in the subendocardial layer. e Tissue segmental perfusion in the mid-equatorial portion of myocardium. Systolic endo-to-epi fluid shift caused by increased pressure can clearly be seen. The inferoseptal segment receives zero flow, due to a lack of vasculature in the region. Slice orientation is shown on the first diagram (anterior/inferior, lateral/septal). f Variation of blood volume in the tissue as a percentage of total reference material volume. Around 1 % maximum variation is observed in the baseline conditions. The deficit in the subendocardial layer lags perfusion rate in time due to capacitance effects

Wave intensity analysis

Wave intensity is defined as the product of the increments in pressure and velocity during a small time interval . The theory of WIA allows the total intensity to be separated into components according to their direction of travel (denoting forward as proximal-to-distal) and classified according to their pushing/compression (increasing pressure) or suction/expansion effects on pressure (Parker 2009). As mechanical disturbance on vessels can result in a propagating wave, coronary wave dynamics are sensitive to the events taking place in the myocardial environment surrounding the vasculature, as well as those internal to the vessel network.

The application of WIA in the human coronary circulation identified six major waves and the associated cardiac events from which they originate (Davies et al. 2006). Of these, the dominant backward-travelling suction wave has since received the most attention, as it signals the forward acceleration of arterial inflow in diastole, during which around 80 % of anterograde flow occurs. The source of this wave has been attributed to the ventricular relaxation which eases the compressive forces acting on the small vessels within the myocardium. In addition, the dominant forward-travelling pushing wave originating from early ventricular ejection has been shown to carry an equivalent magnitude of cumulative wave intensity (also referred to as wave energy, wave area or the time integral of wave intensity (Siebes et al. 2009).

The six major waves can be identified in the WI profile calculated from the baseline simulations, in Fig. 4. The WI shape was observed to exhibit a dependence on the measuring location. Reported results were obtained from a location approximately 6 cm along the LAD. The % contributions of the individual cumulative wave intensities compare well with those reported in Davies et al. (2006), as summarised in Table 2.

Fig. 4

Wave intensity profiles and coronary pressure and velocity calculated from baseline results. The waves associated with anterograde acceleration of flow are coloured in black. The six major waves are labelled in accordance with Table  2

Table 2

Cumulative proportional wave intensity: comparing model results with experimental measurements of Davies et al. (2006)

Wave-type	Experimental %	Model %
① Early backward pushing wave	1.9±2.1	5.1
② Dominant forward pushing wave	22.3±7.9	27.4
③ Late backward pushing wave	20.5±2.9	2.8
④ Forward suction wave	18.9±4.0	13.4
⑤ Backward suction wave	30±5.7	37.3
⑥ Late forward pushing wave	6.1±2.4	14.1

Wave sensitivity to model parameters

The dependency of wave characteristics on the underlying parameters was investigated through local perturbation around the baseline conditions. Specifically, we examined the effects of contractile (active tension), haemodynamic (viscous resistance, wave speed and outflow), porous (outflow and pressure–volume characteristics) and valve dynamic (rate of transition) parameters on coronary waves. The modified wave intensities and WI indices under the altered parameter sets are summarised in Fig. 5. The top row summarises the perturbation applied to the baseline parameters.

Fig. 5

Modified coronary waveforms under parameter perturbation. The applied changes are illustrated in the top row, and results of adjusting QRS duration (first column), active tension transient (second column), vascular resistance (third column), distal outflow pressure (fourth column) and aortic valve transition rate (fifth column) are shown below. The bottom row shows the change in % area of dominant forward pushing wave (DFPW) and backward suction wave (BSW). The grey patches show the baseline results for comparison. Refer to Sect. 4.2.1 for detailed descriptions

Discussion

Baseline parameterisation

The baseline pressure results in Fig. 3a reveal a subtle but important difference between the LV cavity pressure and the microvascular pore pressure , which is that the peak of the two pressures does not necessarily coincide in time. This is significant because, as wave-generating mechanisms, these pressures act on opposite locations of the coronary network and at times can act in an independent manner. The shape of can be modified through the systemic haemodynamic parameters and under intact aortic valve function is a dominant driving force for the coronary inlet pressure as the figure shows. On the other hand, during systole, is governed mainly by the intrinsic myocardial contraction. Therefore, in early systole, the rate of rise of dominates the forward initiated coronary waves, but in late systole, as systemic capacitance begins to become exhausted, stabilises while the intramyocardial stress continues to rise at near-constant LV pressure. Thus, the shape of in late systole follows largely that of the tension transient (e.g. see Fig. 5). It is the sharp drop in the intramyocardial stress which enables the generation of the backward suction wave. This mechanism is examined further below in Sect. 4.2.2.

In the literature, several mechanisms have been proposed to characterise the flow impeding forces exerted on coronary vessels by the myocardium, although no single mechanism has been shown to fully explain the experimental observations (Westerhof et al. 2006). In a recent model-based screening investigation, it was suggested that a combination of two forces—shortening-induced intracellular pressure (SIP) and cavity-induced extracellular pressure (CEP)—can qualitatively account for a broad range of existing data (Algranati et al. 2010). This investigation was, however, limited by the phenomenological manner in which each mechanism and their interaction were formulated. CEP was prescribed to be a linear function of the wall depth according to experimental observations; however, the transmural position of vessels was idealised instead of being based on a measured distribution. SIP was prescribed by a simple scaling of the LV volume waveform (Algranati et al. 2010) and applied without variation in local activation timing or fibre architecture. This was replaced by a scaled LV elastance waveform in later work (Mynard et al. 2014) for a more realistic coronary flow. The two pressures were regarded to be independent and linearly summed. While these assumptions assist our intuitive understanding, they can also affect the resulting wave profile in unknown ways. Both SIP and CEP originate from cellular contraction, and their resultant manifestation is subject to force equilibrium under a given set of boundary conditions. Our model accounts for these physical principles through multiscale modelling of conservation laws. With this, we showed from a biophysical basis that the maximal pressure rise from regionally activating myocytes corresponds to late systolic wave features in coronary flow such as the second dip in velocity, while the maximum cavity pressure is associated with the initial dip (Fig. 3a, b).

The total inflow and outflow from the one-dimensional arterial network shown in Fig. 3b reflect another mechanism of interaction, between coronary flow and contraction. While both traces feature a sharp drop during the isovolumic phases, at isovolumic contraction inflow is greater than the outflow, while the trend reverses at isovolumic relaxation. This dynamic is a product of the vascular compliance (both proximal and microvascular), as initially described by the intramyocardial pump model (Spaan 1985). When the total flow is decomposed by the supplied myocardial layer, the transmurally varying influence of the contraction of coronary flow is clearly seen (Fig. 3d). Opposite trends are observed during systole in the subendocardial (reversed) and subepicardial (augmented) flows, while the mid-myocardial layer features a relatively constant anterograde flow. Since the total inflow is always positive, this trend can be interpreted as a cardiac phase-dependent redistribution of flow between the layers. This mechanism has been reported by experimental studies (Rovai et al. 1989) and is also reflected by the segment-wise flow and fluid mass variation as shown in Fig. 3e, f. The time lag of fluid mass behind flow is another manifestation of the vascular compliance, showing for example that forward subendocardial flow which begins during early diastole requires until late diastole to replenish the fluid deficit inherited from the previous systole.

Further qualitative agreements are found with regards to experimental observations, namely that may exceed locally (Westerhof 1990), features a dicrotic notch, and the reversal of cardiac outflow following aortic valve closure. On the other hand, the difference between and is less than the range observed experimentally, leading to an absence of a crossover point of the two pressures at late systole. This indicates an insignificant pressure gradient across the aortic valve (and into coronary ostium) resulting from suboptimal valve model characterisation. However, the cWIA results are expected to be minimally affected by these since the time course of dP remains largely unchanged.

In the current analysis, the proportional cumulative areas of the individual waves (‘% wave area’) are used as the principal index for comparison with experimental data due to its independence on the cardiac load. Furthermore, our previous work has shown that peak magnitude of the wave is highly sensitive to the standard post-processing procedure, and thus is unsuitable for robust comparison (Rivolo et al. 2014).

The baseline WIA (Fig. 4; Table  2) shows overall an accurate model reproduction of the measured % area of individual waves. The two largest waves (⑤ backward suction wave & ② forward pushing wave) occupy a greater proportion than the mean values reported in Davies et al. (2006); however, given that their ratio is almost identical to the reported value (), the difference may be attributed to the % errors in the other, less significant waves. The largest difference was found in the ③ late backward pushing wave, which is a reflection of the dominant forward pushing wave and thus is the wave which is most sensitive to the internal wave transmission characteristics of the vascular network. The underestimation of this wave is expected to have artificially inflated the % of all other waves. Often in the literature, it has been found that the early and late backward pushing waves (① and ③) cannot be reliably distinguished, and thus only the combined value was reported (Davies et al. 2011; Rolandi et al. 2012; Silva et al. 2013). With this approach, the model currently underestimates the observation by about ( vs ). The ⑥ late forward pushing wave was overestimated, but expectedly due to the fully reflecting boundary condition imposed at the coronary ostium.

The % wave areas can be grouped and summed according to contraction (①, ②, ③) versus relaxation phase, accelerating (②, ⑤, ⑥) versus decelerating, or forward (②, ④, ⑥) versus backward-travelling tendencies. As with the individual % areas these measures exhibit a dependence on the sampling location along the vessel. Figure 6 shows the variation of these indices along the distal half of the LAD where, notably, the total forward-travelling % wave area can be seen to diminish with distance. Conversely this implies that the total % wave energy balance is shifting from the forward to the backward waves distally towards the microcirculation.

Fig. 6

Changes in % area indices along distal LAD

These results must, however, be interpreted with caution and experimentally verified, due to the lack of specificity in network parameter tuning in the current model. The vascular mesh employed in this work lacks tapering as explained, which would result in a lack of wave reflections. In addition, wave speed governs the characteristic impedance of each vessel and, under a network setting a mismatch of impedances at a junction will lead to unphysiological wave reflections. Using a linearised analysis, the reflection coefficient of a pressure perturbation can be defined as (Sherwin et al. 2003)where the superscript refers to the parent (1) or daughter segments (2,3). While it has long been suggested that coronary junctions are well matched for forward-travelling waves ( Arts et al 1979), the more recently identified prominence of the backward-travelling waves and the resultant trade-off between forward and backward matching have not been reconciled experimentally into a unified description. The forward and backward reflection coefficients under the constant wave speed regime are illustrated in Fig. 7. Coincidentally, the forward reflection coefficients are predominantly near 0 under this simple assumption while the backward-travelling waves mostly encounter negative coefficients due to the expanding area ratio. Our preliminary investigations indicate that when the wave speed is allowed to vary segment-to-segment, a reduction in the backward reflectances while simultaneously preserving the forward reflectance requires a significant progressive increase in the distal wave speed, unless one of the daughter vessels is of much larger calibre than the other.

Fig. 7

Forward (left) and backward (right) linear reflection coefficients under a constant wave speed of

Modulators of coronary waves

Contraction parameters The wave intensity profiles were sensitive to both the shape of the tension transient as well as the QRS duration (QRSd). However, the dependence on QRSd was only an indirect one through modulation of the aortic sinus pressure. The reduction of QRSd to did not cause a significant difference in the shape of the contraction waves (①-③), neither was the timing of these waves altered as the timing of ejection was preserved under new activation sequence (Fig. 5, first column), whereas the accelerated rate of relaxation led to an augmented backward suction wave. On the other hand, a faster contraction and slower relaxation led to a change in all six waves (second column), correspondingly increasing the % area of contraction waves and reducing relaxation waves. An expanded analysis is described in Sect. 4.2.2.

Haemodynamic parameters A consequence of modifying in the current model is the modulation of the friction coefficient through the factor [see Eq. (48)]. The result of perturbation to from 1.05 was an approximate halving of the friction coefficient (third column), which increased the magnitude of all waves. This change was caused mostly through the increased flow rate in the vessels throughout the whole cycle, but affected the relaxation waves most significantly since the vascular resistance was at its minimum and the most rapid inflow gradient could be produced then. The effect of increasing the venous outflow pressure was much less pronounced (fourth column), with a near-uniform downward and upward shift in coronary velocity and pressure respectively but no change in the wave magnitudes. This is due to the fact that wave intensity is solely a function of the rates of change in velocity and pressure and has no dependence on their magnitude. Modifying the wave speed of the vessels to resulted in an insignificant change to the wave intensity as a uniform scaling of wave speed preserves reflection coefficients as can be seen in (65). Also, due to the comparatively short vascular path lengths even a change of in the wave speed translated to only a few extra milliseconds in the wave propagation times.

Aortic valve parameters The modification of the aortic valve flow dynamics in Eq. (56) was a necessary element in tuning the forward wave intensities. Clearly, capturing the dynamics of the aortic sinus/valve region with a zero-dimensional model imposes limitations—however, the modification does, to an extent, alleviate an otherwise serious flaw in the model. The results of the fifth column emulate the wave dynamics in the absence of this modification, through reducing the time constant by a factor of 10. The abrupt closure of the aortic valve leads to brief, but much steeper rates of rise in both coronary pressure and velocity, leading to around a fortyfold increase in the peak of the late forward pushing wave. While the changes to other waves are much less pronounced, their % areas collectively reduce due to the greatly expanded proportion of the forward pushing wave.

Contraction, relaxation and coronary waves

The results of the parameter sensitivity analysis indicate that coronary waves are strongly dependent on the relative rates of activation and relaxation, the synchrony of contraction and the time course of aortic sinus pressure. The codependency of these factors is further examined in Fig. 8 through simultaneous variation of the tension transient and QRS duration. The active tension parameters were adjusted to shift the timing of maximal tension by of the twitch duration, while keeping its magnitude constant. Figure 8a shows the changes in coronary pressure at the WI sampling location, where it can be seen that the time course of pressure is disproportionately sensitive to an advancement of peak tension in time, as opposed to when it is delayed. The effect of QRSd was more moderate, affecting most significantly the onset of the diastolic reduction in coronary pressure. These changes were directly reflected by the forward pushing wave % area (Fig. 8b), which exhibited a monotonic dependence on the tension transient (but with a greater sensitivity at shift). It may appear counterintuitive that forward pushing wave decreases with improved synchrony—in terms of the actual magnitude, there was only a slight dependence on QRSd (possibly owing to the fact that activation gradient was imposed only transmurally in this model), and the observed trend was mostly a result of the expanded proportion of the waves that occur during relaxation.

Fig. 8

Wave sensitivity to contraction parameters a coronary pressure sampled at 6 cm along LAD. Modifying the time course of active tension resulted in a much greater effect on the pressure trace compared to changes in QRSd. Shown in dotted lines are the timing of the maximal tension, for a myocyte which is activated at  s b forward pushing wave (FPW) % area shows a relatively linear dependence on both peak tension time and activation synchrony c backward suction wave (BSW) % area is positively correlated with peak tension time under synchronous activation, but the maximal rise associated with the rapid relaxation rate is lost with increasing level of dyssynchrony

In contrast, the dependence of backward suction wave % area on QRSd (Fig. 8c) was accompanied by changes in the wave magnitude and, under all variations of the tension transient, exhibited a reduction with increasing dyssynchrony. Furthermore, increasing level of dyssynchrony was found to be capable of overcoming even the accelerated relaxation rate of individual myocytes thereby reducing the % area of the backward suction wave as shown in Fig. 8c, top right row. Such a trend was not seen in the forward pushing wave, and highlights the stronger dependence of the backward waves on the intramyocardial stresses.

Increase in ⑤ backward suction wave has been reported to be associated with optimal biventricular pacing and the associated coronary flow increase in patients with dyssynchronous heart failure (Kyriacou et al. 2012). But the link between myocardial perfusion, coronary flow, and dyssynchronous heart failure is poorly understood at present, and the study of Kyriacou et al. (2012) remains the sole application of cWIA to heart failure to date. The key challenge is presented by the difficulty in separating the pathophysiological response from normal physiological response that arises as a result of altered electromechanics (Claridge et al. 2015). In this regard, the proposed model offers a new method by which to address these unanswered questions. They will be pursued in future investigations.

Modulators of perfusion

Although WI was found to depend only weakly on the distal pressure , both perfusion and fluid accumulation in the porous domain were significantly affected. An interesting observation, when was raised as shown in Fig. 5 (fourth column), was that the backward suction wave remained unmodified even when the pores were in a state of dilation at its onset. Under the current model assumptions, this implies a dominant role of the myocardial relaxation rate in producing the early diastolic suction over factors pertaining to the intracoronary haemodynamics. That is, wave behaviour was largely independent of the venous load.

The baseline mean perfusion rate of observed is significantly higher than the 1–2 range typically found in resting porcine experimental conditions (Fedor et al. 1978). However, the baseline parameterisation of was tuned to experimentally observed venous pressures, rather than to compensate for the lumped nature of the distal circulation in the current model. Indeed, when was raised to , more in line with the experimentally observed zero flow pressure (Bellamy 1978), the mean perfusion was reduced to . In addition, the pressure drop along LAD was brought closer to a physiological level . Likewise, in the absence of anatomical details, the baseline coupling interface resistances ( in (54)) were prescribed with a simple objective of minimising the uncontrolled wave reflections. When the observation that approximately half of the total pressure drop occurs between the arteriolar and capillary compartments under normal vasomotor tone (Chilian et al. 1989) is approximated by a tenfold increase in along with raised , physiological mean perfusion and pressure drop along LAD are recovered . These changes occurred largely free of deviations from the baseline WI, except for a doubling of the ③ late backward pushing wave % area (see Table 2).

The consequence of modifying the porous compartment parameters was contrary to our expectation. The pressure–volume relationship of the porous domain was adjusted to reflect the arteriolar, rather than capillary properties as shown in Fig. 9a. As can be seen in 9b, its major effect was to raise the systolic subendocardial flow and diastolic subepicardial flow, leading to a sustained increase in the total coronary flow throughout the cardiac cycle. This was caused by the greater negative pressure generation in the compressed pores, thereby explaining the observed spatiotemporal variation in perfusion.

Fig. 9

Effects of porous medium parameters a modified pressure–volume relation, b The modifications lead to an increase in systolic subendocardial and diastolic subepicardial flows. Shown in grey are baseline results for comparison

Where accurate reproduction of regional perfusion is concerned, these observations present a case for employing a multi-compartment porous domain framework (Michler et al. 2012; Hyde et al. 2014), which can distinguish the pressure and flow across the multiscale hierarchy of distal vessels. With a single-compartment model, reproduction of wave dynamics may require a compromise in the distal characteristics. However, we maintain that for WIA applications, the single-compartment model can sufficiently reproduce a broad range of behaviours. This is because waves which originate in the deep distal vessels have a limited influence over the waveforms observed in the large arteries due to wave-trapping and reflection effects (Davies et al. 2012).

Age, wave speed and WIA

The present model-based investigations largely reaffirm the previously inferred wave-generating mechanisms in the experimental studies. In addition, further clarifying remarks can be made with respect to the age-related changes in wave characteristics. Even though age-related changes in coronary arteries may involve a myriad of mechanisms, wall stiffening had received a central emphasis. Although there is no reported evidence in the coronary arteries as yet, the increase in the pulse wave speed in the aorta with age is a well-documented phenomenon which may reasonably be expected to take place in the coronaries also. This being the case, Davies et al. (2006) found no correlation between the ⑤ backward suction wave and age, whereas a significant reduction in its % area was observed in patients with LV hypertrophy. Both of these observations are consistent with the simulation results presented here and reinforce the view that altered initiation of the suction wave by myocardial relaxation has a greater repercussion than its modifications during transmission to the proximal artery.

On the other hand, Davies et al. (2006) reported that ageing correlated positively with the ④ forward suction wave magnitude and offered two possible explanations. The ventricular–arterial hypothesis attributed the increased coronary stiffness and its ability to transmit a greater magnitude of energy into the vasculature. The aortic hypothesis ascribed the greater rates of rise, and correspondingly, the decline in the aortic pressure brought on by its own stiffening. The present work lends support to the dominant role of the second mechanism, especially given the absence of reported augmentation in the other forward-travelling waves.

Shortcomings of the model

The present work is subjected to several shortcomings which should be addressed in future investigations. Of these, the design of a suitable constitutive law will require perhaps the most attention, as anisotropic poroelastic response in large strain regime is currently a sparsely explored area even in the field of biomedical engineering applications. Formal homogenisation approaches provide an alternative avenue, through which macroscopic description can be obtained from well-defined microscopic behaviours (see Davit et al. (2013) for an overview of methodologies). Even if these methods may lead to a computationally intensive coupled micro–macro problem, they may still play an important role in the validation of the practical approach adopted in this work.

The proposed model requires tuning of a large number of parameters to govern its behaviour. While global functional indices were shown to be physiological in this work, further experimental validation is key if the model is to be extended to predictive applications. Vascular network parameterisation presents additional difficulties in this regard due to their limited accessibility, and thus, theoretical developments are currently under way to complement these efforts. Balancing the distal resistance parameters (which modulate perfusion) with wave reflection behaviour on a vessel-wise basis was a challenging task, and should be explored further in future work.

In addition to the mechanisms described in Sect. 4.1, there are additional ones by which the coronary–myocardial coupling could take place that we have not explored in this work. Two prominent examples are (1) the direct extravascular compression on the explicitly represented vessels and (2) the effect of vessel motion on the flow and waveforms. However, we do not believe the omission of these mechanisms has resulted in significant errors for the following reasons. Of the 94,000 vessel segments reconstructed from the data, only proximal segments were retained for simulation.

Of these, those situated along the epicardium can be disregarded since it is known that systolic compression is insignificant there (Toyota et al. 2005). Considering the significant distal expansion in the number of vessels and total cross-sectional area (Kassab et al. 1993), the one or two generations of intramural vessels that are currently unaccounted for representing a tiny fraction of the total surface area available in the coronary network for vessel–myocardial interaction. Temporally, these transmurally extended vessels are  cm in length, such that waves initiating at the endocardium would experience only a -ms delay, which is less than half the sampling interval at which most clinical WIA are conducted (Rivolo et al. 2014). Taken together, it is clear that in our model, the majority of the systolic compression occurs in, and is captured by, the region of vascular network represented by the poroelastic domain.

The impact of vessel motion is convoluted by numerous ways in which it may exert its influence on the vascular flow. Multiple simulation studies have revealed that the frequency and phase of the vessel motion relative to flow have the greatest impact, whereas the curvature variation of vessel segments was found to be of minimal consequence. The inertial effects were shown to be significant at 5 Hz vessel motion (Pivkin et al. 2005), whereas at 1 Hz they were negligible (Santamarina et al. 1998). The justification for applying motion at nearly five times the heart rate came from a frequency domain analysis of biplane-cineangiography images (Gross and Friedman 1998), although a later analysis of the data showed that around of the power was contained at or under 3 Hz (Moore et al. 2001), implying the modelling results are likely to have been exaggerated to a degree. In addition, the relative phase between vessel motion and flow was reported to be a major determinant of the resulting flow splitting at a bifurcation (Pivkin et al. 2005). The maximal effect was found when the flow was ahead of the outward ventricular motion, at which the flow ratio of the daughter vessels deviated around from the static myocardium / steady flow case. The deviation at phase difference was less than 5 %. We estimate the phase difference in our results to be in the range since the majority of the anterograde arterial flow acceleration occurs during early diastole, to be followed by the bulk of the ventricular inflow. When compared to the beat-to-beat variations in cathlab measurements and errors that may occur during WIA post-processing (Rivolo et al. 2014), we conclude that the existing evidence for cardiac-motion-induced effects in coronary flow indicates a secondary role that is within the order of magnitude of experimental error.

Patient-specific modelling

While the translation of the present model framework to the clinics is beset by challenges in both data acquisition and modelling at current time, continuing advances in the medical imaging field present optimistic prospects. The recent introduction of 320-multidetector CT scanner combined with nonlinear iterative reconstruction now allows low-dose/high-resolution angiography for coronary vessel characterisation (Williams et al. 2013). Likewise, perfusion MR imaging is solidifying its gold standard status through phantom validation (Chiribiri et al. 2013) and patient trials (Greenwood et al. 2012), permitting high-resolution quantification of myocardial blood flow as recently demonstrated (Villa et al. 2015). Meanwhile, the scope of validation is currently being expanded through first-in-human simultaneous measurements of coronary pressure/velocity and ventricular pressure/volume by our group, which reveal a wealth of information regarding coronary–myocardial coupling. A reliable personalisation of a model is conditional upon the conformance between obtainable data and the complexity of the model. Although the eventual form of the model will likely entail a further reduction to the clinical essentials, it is likely to benefit from the closing gap between the data and model, as well as inverse parameter estimation techniques that have been successfully employed to myocardial mechanical characterisation (Asner et al. 2015).

Conclusions

In this work, an in silico platform for coronary wave intensity analysis was developed and used to investigate the wave dependence on individual parameters underlying contraction, perfusion and systemic haemodynamic processes. To the authors’ knowledge, this is the first and only instance of such a study based on an integrated biophysical model of cardiac perfusion has been reported. While experimental coronary WIA is often afflicted by temporal waveform misalignment and motion-related issues which degrade the data quality, the present work offers a noise-free, repeatable and quantitative alternative. In addition, it overcomes the technological limitations in the bandwidth of measurements that can be simultaneously acquired, opening up new avenues in diagnostic index discovery and predictive clinical applications.

The simulation studies of wave mechanisms highlighted the direct and indirect relationships between variables underlying ventricular–aortic–coronary coupling. These findings will aid in the interpretation of clinical studies with respect to mechanisms of pathophysiology. Several milestones must be achieved prior to exploiting the full potential of the present work, which is the diagnostic application to patient data. Fine-tuning of vascular network parameters and validation against experimental data represent two immediate challenges that lie ahead. Efforts to address these issues are currently under way, through controlled in vivo experiments which build on the work of Schuster et al. (2010).