A physiologically realistic virtual patient database for the study of arterial haemodynamics
AA physiologically realistic virtual patient database for the study ofarterial haemodynamics
Gareth Jones , Jim Parr , Perumal Nithiarasu , and Sanjay Pant † College of Engineering, Swansea University, Swansea, United Kingdom. McLaren Technology Centre, Woking, United Kingdom. † Corresponding author: [email protected]
Abstract
This study creates a physiologically realistic virtual patient database (VPD), representing the human arterialsystem, for the primary purpose of studying the affects of arterial disease on haemodynamics. A low dimen-sional representation of an anatomically detailed arterial network is outlined, and a physiologically realisticposterior distribution for its parameters constructed through the use of a Bayesian approach. This approachcombines both physiological/geometrical constraints and the available measurements reported in the litera-ture. A key contribution of this work is to present a framework for including all such available information forthe creation of virtual patients (VPs). The Markov Chain Monte Carlo (MCMC) method is used to samplerandom VPs from this posterior distribution, and the pressure and flow-rate profiles associated with each VPcomputed through a physics based model of pulse wave propagation. This combination of the arterial networkparameters (representing a virtual patient) and the haemodynamics waveforms of pressure and flow-rates atvarious locations (representing functional response and potential measurements that can be acquired in thevirtual patient) makes up the VPD. While 75,000 VPs are sampled from the posterior distribution, 10,000 arediscarded as the initial burn-in period of the MCMC sampler. A further 12,857 VPs are subsequently removeddue to the presence of negative average flow-rate, reducing the VPD to 52,143. Due to undesirable behaviourobserved in some VPs—asymmetric under- and over-damped pressure and flow-rate profiles in left and rightsides of the arterial system—a filter is proposed to remove VPs showing such behaviour. Post application ofthe filter, the VPD has 28,868 subjects. It is shown that the methodology is appropriate by comparing theVPD statistics to those reported in literature across real populations. Generally, a good agreement between thetwo is found while respecting physiological/geometrical constraints. The pre-filter database is made availableat https://doi.org/10.5281/zenodo.4549764 . Keywords— virtual patients, stenosis, aneurysm, pulse wave haemodynamics, MCMC, screening, virtual pa-tient database
Arterial disease refers to any disease affecting the arterial system. Two of the most common forms of arterialdisease are stenosis and aneurysm. These two forms of disease are estimated to affect between 1% and 20% ofthe population [1–4]. Ruptured abdominal aortic aneurysms alone are estimated to be the cause of between 6,000and 8,000 deaths a year within the United Kingdom [5]. Current methods for the detection of arterial diseaseare primarily based on imaging techniques which can be expensive and thus unsuitable for mass screening. Apotential solution is offered by the recent emergence of data-driven Machine Learning (ML) algorithms whichcan be trained on easily acquirable pressure and flow-rate measurements. The key idea is that the presence ofdisease—both stenosis (a reduction in vessel area) and aneurysm (an increase in vessel area)—will create changesin such measurements, which may be detectable by the ML algorithms. To explore this possibility, however, alarge database of both healthy and diseased subjects along with their respective measurements is necessary. Thisstudy presents a computational method for the creation of such a physiologically realistic database.As opposed to acquiring measurements from a real population, which may be impractical, a virtual patientdatabase (VPD) based on a validated physics-based model is proposed. This idea of exploiting virtual patients(VPs) to gain meaningful information has been previously explored in the context of human arterial networks [6–8],and in the context of path clamp physiology [9]. In [6], a VPD of healthy subjects is created by varying sevenparameters of a large 55-artery network to assess the relation between foot-to-foot pulse wave velocities and aorticstiffness. In [7], 1D arterial network parameters are scaled through empirical relations based on weight, height, andage to create a virtual patient database. In [8], nine parameters of a small network of abdominal aortic bifurcationare varied to create the VPD and ML algorithms subsequently applied to detect the presence of stenosis. Thelimitations of these studies are that either not all the parameters are varied in both healthy and diseased states, or1 a r X i v : . [ phy s i c s . m e d - ph ] F e b priori distributions are assumed for the parameters. Both of these are important considerations for the creationof a general purpose realistic VPD which captures the variability as close as possible to the real population.Furthermore, the VPD must account for known variation of the measurements reported in the literature alongwith the physical constraints, while making minimal assumptions. This study presents a probabilistic frameworkfor creation of such a VPD based on a previously proposed, anatomically correct, 56-vessel and 77-segment arterialnetwork [10, 11]. The proposed framework may also be useful for the novel and active research area of in-silicoclinical trials [12], for which generation of realistic virtual patients is a key consideration.This article is organised as follows. First, a brief description of the physics-based pulse wave propagationmodel is presented. This is followed by a description of the arterial network topology—vessels, locations of diseaseand potential measurements—and its reduction. Next, parameterisation of the arterial network—geometry, vesselproperties, and boundary conditions—is presented, followed by the Bayesian framework for the VPD creation.Finally, the solution to the Bayesian problem through Markov Chain Monte Carlo method is presented, beforediscussing the results and statistics of the VPD. A validated one dimensional model of pulse wave propagation [13, 14] is adopted in this study. The model describesthe arterial vessels as one-dimensional compliant tubes, i.e. all the properties are described by a single coordinate x along the length of the vessel. Blood flow is assumed to be Newtonian, laminar, incompressible, and axisymmetric.With the arterial cross-sectional area represented as A , and considering the average lumen pressure P and theaverage velocity U to be the primary variables, the mass and momentum conservation equations at any coordinate x can be, respectively, written as ∂A∂t + ∂ ( U A ) ∂x = 0 , (1) ∂U∂t + U ∂U∂x + 1 ρ ∂P∂x = fρA , (2)where ρ represents blood density and f represents the frictional force per unit length. Depending on the assumptionof velocity profile across the arterial cross section, f can be written as f ( x, t ) = − ζ + 2) µπU, (3)where ζ is a velocity-profile dependant constant, and µ represents dynamic viscosity of blood. The system ofequations is closed by specification of the mechanical behaviour of the wall, commonly referred to as the “tubelaw”. A commonly used form [14] is P − P ext = P d + βA d ( √ A − (cid:112) A d ) , (4)with β = Eh √ π, (5)where P ext represents the external pressure, P d represents the diastolic blood pressure, A d represents the diastolicarea of the vessel, E represents the vessel wall’s Young’s modules, and h represents the vessel wall’s thickness.This system of equations has been widely employed, tested extensively, and validated against more complex three-dimensional models [13–18].Since it is computationally impractical to model every single artery in the circulatory network, it is commonto truncate the 1D network by only considering relatively large arteries which are of interest. The effect of allthe vessels that come before and after the network are then incorporated as appropriate boundary conditions.Typically, at the inlet of the network a volumetric flow-rate is specified, and at the outlets Windkessel models [19]are employed. A three-element Windkessel model, shown in Figure 1, comprises of two resistors, R and R , whichreplicate the viscous resistance of large arteries and the micro-vascular system, respectively, and a capacitor C ,which replicate the compliance of large arteries.The 1D model is solved numerically by a subdomain collocation (SDC) method [20]. This scheme is chosen for itscomputational efficiency. The SDC implementation is validated against a discontinuous Galerkin (DCG) scheme[13], which in turn has been successfully validated against a 3D model of blood-flow through stenosed arterialvessels [21]. The arterial network used as a basis for the creation of VPs is based on a version of the ADAN network [14].This network contains 56 vessels within the human arterial system, described by 77 arterial vessel segments. The2 P P out R R C Figure 1: Configuration of a three-element Windkessel model. Q and P represent the volumetric flow-rateand pressure at the terminal boundary of the 1D system respectively.topology of this network is shown in Figure 2, and each vessel within the network is identified in Appendix A. Thisnetwork is described by: •
38 vessel segments with constant reference radii ( r ) and vessel wall mechanical properties ( β ) along theirlengths. The description of each of these segments requires specification of three parameters: the length L ,reference radius r , and mechanical property β of the segment. •
39 vessel segments with linearly varying r and β along their lengths. The description of each of thesesegments requires specification of five parameters: the length L ; the reference radii r at the proximal and thedistal ends of the segment; and β at the proximal and distal ends. •
31 Windkessel models.
Every Windkessel model at the outlets requires specification of 2 resistances and acompliance. • One inlet flow-rate profile.
In this study, the inlet flow-rate is parameterised by a 5 th order Fourier Series(see Section 4.3.1 and [8]). This requires specification of 11 coefficients.Thus, to describe a VP through direct specification of all the above parameters results in the total dimensionalityof 413 per VP. This high dimensionality is problematic leading to increased complexity in creating VPs [22]. Inwhat follows, methods to reduce the dimensionality of VP description are presented. These are primarily related toeither reducing the network or employing a parsimonious parameterisation. Since the primary purpose of the VPDis deployment of ML methods for screening of aneurysm and stenosis through easily acquirable measurements,it is important that the reduction in dimensionality does not compromise i) the locations where disease andmeasurements are possible, and ii) the precision and variability of measurements at such locations. A review of literature is carried out to understand both where disease is likely to occur and where measurementscan be obtained. The latter are restricted to locations at which continuous profiles can be recorded non-invasively.While arterial disease can occur in a large number of vessels, for this study vessels with only high prevalence ofstenosis or aneurysm are considered.
Based on literature, the following locations where measurements can be taken are identified:
Pressure in the radial and common carotid arteries:
Continuous non-invasive arterial blood pressure pro-files can be obtained in the radial and common carotid arteries using applanation tonometry [23, 24]. Applanationtonometry is already clinically used [25]. The right and left radial and common carotid arteries are identified inFigure 2 as vessels 8 and 22; and 5 and 14, respectively.
Pressure in the brachial arteries:
It is possible to estimate continuous blood pressure at the brachial arteriesthrough reconstruction of finger arterial pressure [26]. Although the use of a model to estimate brachial pressurewill introduce errors, these recreated brachial pressure profiles have been shown to meet the requirements set bythe Association for the Advancement of Medical Instrumentation [26, 27]. The right and left brachial arteries areidentified in Figure 2 as vessels 7c and 21c respectively.
Flow-rate in brachial, carotid, and femoral arteries:
Using Doppler ultrasound techniques it is possibleto measure blood velocity in both the upper and lower extremities. Doppler ultrasound has been shown to workon the brachial [28], common carotid [29], and femoral [30] arteries. The first and second, right and left femoralarterial segments are identified in Figure 2 as vessels 50b and 53a; and 56b and 59a respectively.3
324 567a7b7c8 910a10b 11 1213 1415 1617 1819a19b 2021a21b21c 222324a24b2526 27 282930 31 3233a33b34 3536373839 4041 424344 45 464748 4950a50b 5152 53a53b54 55a55b55c 56a56b57 5859a59b 6061a61b61c PressureFlow-rateDisease
Figure 2: The connectivity of the reference arterial network, taken from [14]. At the inlet (free end) of vessel 1,a volumetric flow-rate is specified and at all the outlets (free ends of the terminal vessels), a Windkessel modelis specified. Locations at which pressure and flow-rates can be measured; and disease is likely to occur are alsohighlighted, see Section 3.2. 4 .2.2 Locations of disease
Four of the most common forms of arterial disease are subclavian artery stenosis (SAS), carotid artery stenosis(CAS), peripheral arterial disease (PAD), and abdominal aortic aneurysm (AAA).
SAS locations:
The ADAN network splits the subclavian arteries into two segments. The right and left instancesof the first and second subclavian artery segments are identified in Figure 2 as vessels 4 and 7a; and 18 and 21arespectively.
CAS locations:
The carotid arteries consist of the common carotid, external carotid, and internal carotid seg-ments. While the narrowing of an arterial vessel can occur within any of the three carotid segments [31], it ischosen to limit its occurrence to the common carotid arteries. The right and left instances of the common carotidarteries are identified in Figure 2 as vessels 5, and 14 respectively.
PAD locations:
The third frequent form of stenosis is PAD, referring to the stenosis of any peripheral vessel.Isolating arterial vessels that are likely to experience stenosis at a high prevalence in patients with PAD is moredifficult relative to SAS or CAS. Both SS and carotid artery stenosis have short and relatively easy to define lists ofpossible vessels they can affect. PAD, on the other hand, can cover a large range of different vessels depending onthe definition prescribed. Studies into PAD primarily focus on the lower extremities [32–34] and thus it is assumedthat VPs with PAD experience change in these vessels. Based on [34] and [33], PAD in VPs is restricted to thecommon iliacs, external iliacs, first and second femoral segments, and the first popliteal segments. The right andleft instances of each of these vessels are identified in Figure 2 as vessels 48 and 49; 50a and 56a; 50b and 56b; 53aand 59a; and 53b and 59b, respectively.
AAA locations:
The ADAN network splits the abdominal aorta into five segments. These five segments areidentified in Figure 2 as vessels 35, 41, 43, 45, and 47 respectively.
In this section, the reduction of the network by removing vessels while preserving the disease and measurementlocations is presented. It is important, however, to ensure that the removal of any vessels does not have a significantimpact on the upstream pressure and flow-rate profiles. Network reduction can be performed by removing terminalvessels and merging them into the lumped parameter (Windkessel) boundary conditions [35]. This process isadopted in this study and is summarised in Figure 3. The following three groups of vessels eligible for removal areidentified: • Group 1:
The first and second splenic segments; the left gastric; and the common hepatic (identified in Figure2 as vessels 37 and 39; 38; and 36). • Group 2:
The common interosseous, the posterior interosseous, and the second ulnar segment (identified withinFigure 2 as vessels 10a and 24a; 10b and 24b; and 11 and 25). • Group 3:
The second popliteal segments; the anterior tibial; the tibiofibular trunk; and the posterior tibial(identified within Figure 2 as vessels 55a; 54; 55b; and 55c on the right side and 61a; 60; 61b; and 61c on theleft side).Three possible reduced networks are proposed, each with the removal of a single group of vessels (group 1, 2, or 3in isolation). The pressure and flow-rate profiles produced by each reduced network are compared against the fullADAN network at all the measurement locations (see Section 3.2). The error metrics proposed in [35] are used.If the maximum error is less than 2%, the full group of vessels are omitted from the arterial network. Otherwise,the vessel segment at the proximal end of the group is re-introduced into the reduced network. This process isiteratively repeated until the maximum number of vessels that can be removed from each group individually arefound. Finally, these vessels from the three groups are removed simultaneously.It is found that the first and second splenic segments; left gastric; common hepatic; and the left and rightposterior interosseous vessels (identified within Figure 2 as vessels vessels 37 and 39; 38; 36; and 10b and 24b)can all be removed from the arterial network without introducing errors larger than 2% relative to the full ADANnetwork. The pressure and flow-rate profiles produced by the full and reduced networks are shown in Figure 4 andthe errors are reported in Table 1. The reduced network is shown in Figure 5 and has a dimensionality of 389, areduction of 24 from the full network. This network and its associated parameters—either directly taken from orcomputed based off of the ADAN network—forms the reference network for this study.5 omputation of lumped parametersCompaction into a singleWindkessel modelTransform to three elements( Q , P ) ( Q , P ) R R C ( Q ) ( P ) ( Q , P ) ( Q ) R v R R CC v ( P )( P ) ( Q )( Q , P ) C new R new ( Q ) ( P ) ( Q , P ) ( P ) R
1, new R
2, new C new ( Q ) Figure 3: The process of incorporating peripheral 1D vessels into the 0D terminal boundary Windkessel modelparameters [35]. Q and P represent the flow-rate and pressure at the proximal end of the vessel that is beingremoved; Q and P represent the flow-rate and pressure at the distal end of the vessel that is being removed; Q and P represent the the flow-rate and pressure at the outlet of the system; R , R and C represent theresistances and compliance of the original terminal Windkessel model; R v , and C v represent the viscous resistanceand compliance of the vessel being removed; R new and C new represent the resistance and compliance of the new2 element Windkessel model after the incorporation of the 1D vessel; and R , R , and C new represent theresistances and compliance of the final 3 element Windkessel model. Details on computation of R v and C v andthen R new and C new can be found in [35]. R
1, new is assumed to be equal to the characteristic impedance of theconnecting vessel (see [35]), and R = R new − R . A parsimonious representation is sought for the aforementioned reference (reduced) network so that its dimension-ality can be further reduced. The reference network has 389 parameters. The vessels in this network can be splitinto the following three categories, which are also indicated in Figure 5: • Category 1:
33 vessel segments with varying β and r , and continuous variation with respect to the prior andsubsequent vessels. • Category 2: six vessel segments with varying β and r , and discontinuous variation with respect to the priorand subsequent vessels. • Category 3:
32 vessel where β and r are constant along their lengths.For Categories 1 and 2 with varying intra-vessel property profile, the reference network assumes a linear variation[14]. Thus any property g ( x ) (reference radius r or β ) varies along length x as: g ( x ) = g − ( g − g L ) xL , (6)where g and g L represent the properties at the proximal end x = 0 and the distal end x = L , respectively. Are-parameterisation of these properties to reduce the total dimensionality is presented next. For all the aforemen-tioned three categories of vessels, a unified parameterisation is proposed, which also results in reduction in the6 P r e ss u r e ( mm H g ) (a) Right comm. carotid 0 1Time (s)6080100 P r e ss u r e ( mm H g ) (b) Right brachial 0 1Time (s)050 F l o w - r a t e ( m l s − ) (c) Right femoral I 0 1Time (s)6080100 P r e ss u r e ( mm H g ) (d) Right radial0 1Time (s)80100 P r e ss u r e ( mm H g ) (e) Left comm. carotid 0 1Time (s)6080100 P r e ss u r e ( mm H g ) (f) Left brachial 0 1Time (s)050 F l o w - r a t e ( m l s − ) (g) Left femoral I 0 1Time (s)6080100 P r e ss u r e ( mm H g ) (h) Left radial0 1Time (s)020 F l o w - r a t e ( m l s − ) (i) Right comm. carotid 0 1Time (s)510 F l o w - r a t e ( m l s − ) (j) Right brachial 0 1Time (s)020 F l o w - r a t e ( m l s − ) (k) Right femoral II0 1Time (s)020 F l o w - r a t e ( m l s − ) (l) Left comm. carotid 0 1Time (s)510 F l o w - r a t e ( m l s − ) (m) Left brachial 0 1Time (s)020 F l o w - r a t e ( m l s − ) (n) Left femoral II FullReduced Figure 4: Comparison of pressure or flow-rate profiles at all the measurement locations taken from the full ADANnetwork and the reduced network produced by the removal of the first and second splenic segments, left gastric,common hepatic, and posterior interosseous. 7 b c de fg hijkl mnopqrs p q r st tu uv vw wx xArm chains (cat. 1)Variable, not chain (cat. 2)Aortic chain (cat. 1) Leg chains (cat. 1)Constant properties (cat. 3)
Figure 5: The reduced network and grouping of vessels. The aortic chain is shown in red and consists of: the firstto forth aortic arch segments denoted by ‘a’ through to ‘d’, the first to sixth thoracic aorta segments denoted by‘e’ through to ‘j’, and the first to fifth abdominal aorta segments denoted by ‘k’ through to ‘o’. The right andleft arm chains are shown in blue and consist of: the first and second subclavian segments denoted by ‘p’ and‘q’ respectively, the axillary artery denoted by ‘r’, and the brachial artery denoted by ‘s’. The right and left legchains are shown in purple and consist of: the external iliac denoted by ‘t’, the first and second femoral segmentsdenoted by ‘u’ and ‘v’ respectively, and the first and second popliteal segments denoted by ‘w’ and ‘x’ respectively.Vessels which have variable properties but are not part of a chain are shown in green. Finally, vessels with constantproperties are shown in yellow. 8 ocation Measurement Average error [%] Systolic error [%] Diastolic error [%]
Right common carotid Pressure 0.4380 0.114 0.316Flow-rate 0.250 -0.085 -0.271Left common carotid Pressure 0.438 0.157 0.318Flow-rate 0.246 -0.111 -0.186Right brachial Pressure 0.643 1.108 0.378Flow-rate 0.875 -1.223 -1.833Left brachial Pressure 0.644 1.179 0.380Flow-rate 0.889 -1.131 -1.788Right radial Pressure 0.654 1.236 0.383Left radial Pressure 0.657 1.231 0.386Right femoral I Flow-rate 0.256 -0.174 -0.114Right femoral II Flow-rate 0.365 0.515 -0.463Left femoral I Flow-rate 0.254 -0.166 -0.115Left femoral II Flow-rate 0.363 0.514 -0.463Table 1: Percentage discrepancy between the measurements taken from the full arterial network and the reducednetwork produced by the removal of the first and second splenic segments; left gastric; common hepatic; andposterior interosseous. Percentage errors have been computed using the error metrics presented in [35].dimensionality associated with the description of Category 1 vessels. This parameterisation includes descriptionof the mechanical properties (Section 4.1), geometric properties (Section 4.2), and boundary conditions (Section4.3).
As opposed to describing the property variations linearly and individually for every vessel segment, see Equation(6), it proposed to use exponential variations, which allows for a simultaneous description of many successive ves-sels, thus reducing the total dimensionality. Such re-parameterisation of the mechanical property β for the threecategories of vessels is as follows: Category 1:
These vessels have continuous variation in properties between successive vessels. Thus, successivevessels can be lumped together into chains and a single parameterisation be adopted. The 33 vessel segments(Category 1 vessels described in Section 4 and labelled in Figure 5 as cat. 1), that meet this description can belumped into the following five chains: • The aortic chain.
This chain of vessel segments includes the first to forth aortic arch segments; the first tosixth thoracic aorta segments; and the first to fifth abdominal aorta segments. • The right and left arm chains.
These chains of vessels include the first and second subclavian segments; theaxillary artery; and the brachial artery. • The right and left leg chains.
These chains of vessels include the external iliac; the first and second femoralsegments; and the first and second popliteal segments.The β profiles for an entire vessel chain can now be described using a single function as opposed to separatefunctions (and hence separate parameters) for each individual vessel. An appropriate choice for such intra-vesselvariation is β ( x ) = β exp( − Ω d x ) , (7)where β ( x ) and β represent the material properties of the chain at the spatial coordinate x and the proximal endof the chain respectively; and Ω d ≥ Category 2:
These vessels do not form part of any chain but show an intra-vessel property variation. For con-sistency and uniformity, their variation is also described by the exponential re-parameterisation of Equation (7).These vessels are labelled in Figure 5 as cat. 2.
Category 3:
These vessels are also not part of a chain and have constant intra-vessel property variation. Again,for uniformity, the re-parameterisation of Equation (7) is used, with the added explicit condition that Ω d = 0.These vessels are labelled in Figure 5 as cat. 3. 9ith the above unified description for the three vessel categories, β properties can be hierarchically assigned. Proximal-distal coherence dictates that the β properties at the distal end of a vessel segment must be greaterthan or equal to the corresponding property at the proximal end of any subsequent daughter segments [11]. Toensure this a hierarchical procedure is adopted. The β properties of the arterial network are initialised by explicitlyassigning a value at the inlet of the first aortic arch segment. The property profile of the aortic chain, and soconsequently the vessel wall mechanical property at the proximal and distal ends of each segment within the chain,is computed using Equation (7). The properties at the proximal end of any vessel segments branching from theaorta are computed by applying a scaling term to the corresponding property at the distal end of the parent aorticsegment: β d = β pL p Ω s , (8)where β d represents the property at the proximal end of the daughter segment, β pL p represents the property at thedistal end of the parent segment, and Ω s represents a daughter-parent scaling term, where:0 < Ω s ≤ . (9)The intra-vessel or intra-chain property profiles of all daughter branches bifurcating from the aorta are computedusing Equation (7). This process is sequentially repeated through each generation of the arterial network until allterminal boundaries are reached. Here, the re-parameterisation of geometric properties—reference radius and length—is presented.
Since the variation of reference radius needs the same requirements as those for the variation of β , the variation of r is described exactly in the same manner as β , as presented in Section 4.1. Lumping arterial vessels into chains,where appropriate, and applying exponential functions to describe the properties β and r of all the vessels reducesthe number of dimensions from 389 in the reference network to 269. Each arterial vessel requires specification of its length. To individually assign a length to every vessel requires71 parameters, accounting for a large proportion of the 269 remaining dimensions. It is thus proposed to reducethe dimensionality of the network by applying a singular scaling term to the lengths of all the vessels. It isempirically found that behaviours, patterns, and overall variability observed in the pressure and flow-rate profileswhen allowing for maximum freedom to the length of vessels, i.e. assigning independent lengths to each arterialvessel, are not lost when applying a singular vessel length scaling term to all the vessels. This analysis justifies theuse of a single scaling term and can be found in Appendix B. With this assumption, the dimensionality reducesfrom 269 to 199.
The boundary conditions consist of the inlet flow rate and the terminal lumped parameters. These are parame-terised as follows:
The volumetric inlet flow-rate to the network, vessel ‘a’ in Figure 5, is described using a Fourier series (FS) [8]: Q inflow ( t ) = N (cid:88) n =0 a n sin( nωt ) + b n cos( nωt ) , (10)where and a n and b n represent the n th sine and cosine FS coefficients, respectively; N = 5 represents the truncationorder; and ω = 2 π/T , with T as the time period of the cardiac cycle. Thus, the following 11 parameters are requiredto specify the inlet flow-rate: { a = 0 , b , a , b , ..., a , b } .10 .3.2 Terminal lumped models There are 29 terminal boundaries in the network. Each of these terminal boundaries is coupled to a Windkesselmodel. Each Windkessel model requires three parameters: two resistances, and a compliance. These Windkesselmodels therefore require specification of a total of 87 parameters.The above description (presented in Sections 4.1, 4.2, and 4.3) has 199 parameters. These represent 22 decayparameters Ω d , 78 daughter-parent scaling parameters Ω s , 87 Windkessel parameters at the outlets, 11 FS coef-ficients describing the inlet flow-rate, and 1 length scaling term. These 199 parameters represent a reduction ofapproximately 48% with respect to the original description with 389 parameters. These parameters are denotedindividually by Ω i , for the i th parameter, and collectively with the vector Ω = [Ω , Ω , · · · , Ω ]. The final step inparameterisation is to describe all these 199 parameters relative to the reference network, and is described next. Instead of describing the networks by directly stating values of the 199 parameters, they are specified relative tothe reference network. Ω = θ T Ω ref , (11)where θ = [ θ , θ , · · · , θ ] T represents the scaling vector in relation to the reference network and Ω ref representsthe vector of reference network parameters. Random sampling from the distribution of θ will result in the VPD.This is described in the next section. The pulse-wave propagation model, denoted as M , takes the parameters θ as inputs and outputs the pressure andflow-rates, collectively denoted as the vector y , at all the locations: y = M ( θ ) , (12)The specific measurements are denoted by τ , which are essentially transformations of y corrupted by measurementnoise E : τ = H ( y ) + E , (13)where H represents the observation operator. Commonly, the observation operator H represents an identity trans-form of selected components of y . The measurement error E is typically assumed to be a zero-mean multivariatenormal: E ∼ N ( , Σ error ) , (14)where is a zero-vector and Σ error represents the error covariance matrix (diagonal when errors are independent).To create the VPD, random samples of θ from their joint probability distribution p( θ ) are required. Thisdistribution p( θ ) should satisfy two broad pieces of information: first, it should satisfy known physiological andgeometric constraints, and second, it should be consistent with what is known about the distributions of measure-ments τ from literature. The latter ensures that the when the physics-based model is run with the parametersdistributed as p( θ ), the resulting distributions of the measurements are close to those reported in real populations. A Bayesian formulation is adopted to find the posterior distribution of the parameters θ . This is suitable becausethe geometric and physiological constraints, in the form a prior, can be combined with the literature reportedmeasurements, in the form of a likelihood, to result in a posterior distribution that considers both the pieces ofavailable information. The Bayes’ theorem allows such a combination naturally:p ( θ | τ ) (cid:124) (cid:123)(cid:122) (cid:125) posterior = likeihood (cid:122) (cid:125)(cid:124) (cid:123) p ( τ | θ ) prior (cid:122) (cid:125)(cid:124) (cid:123) p ( θ )p ( τ ) (cid:124) (cid:123)(cid:122) (cid:125) evidence . (15)The prior distributions p( θ ) are described in Section 5.2 and are chosen to be either uninformative or weaklyinformative. This avoids imposition of strong prior beliefs which may influence the posterior. Their primarypurpose is to impose appropriate bounds on the supports of the distributions based on geometric and physiologicalconstraints. For the computation of the likelihood, the literature reported measurements are first representedas multivariate Normal distribution with mean µ lit and covariance Σ lit . Then, in the context of the statisticalmodel, Equations (12)–(15), it is assumed that τ is measured to be equal to τ = µ lit and the error covariance is Σ error = Σ lit . The likelihood for individual measurements is described in Section 5.3. Finally, in Equation (15),the evidence term is a normalising constant, which is not of interest and need not be evaluated.11 .2 Prior distributions Note that the network parameterisation is relative to the reference network, see Equation (11), and hence the dis-tributions are for the scaling terms θ . Since scarce information is available about the variations of the parameters,the prior distributions are primarily constructed based on the bounds that should be observed in the parameters.In future studies, another approach which may be employed to assign tighter prior distributions could be the useof empirical relationships between the parameters and weight, age, height, etc., see for example [7]. The followingthree types of uninformative or weakly informative priors [36] are employed:
1. Bounded parameters:
A uniform prior distribution is chosen for all the scaling parameters which aresupported on a bounded interval. The daughter-parent ratio, see Equations (8), is the only example of such aparameter. For the i th pair of vessels with reference parent-daughter ratio Ω i, ref , it can be seen from Equations (9)and (11) that the corresponding scaling term θ i is bounded between 0 and 1 / Ω i, ref . The prior probability densityfunction (PDF) is thus: p( θ i ) = (cid:40) Ω i, ref if θ i ∈ (0 , / Ω i, ref ]0 otherwise , (16)and shown in the left plot of Figure 6.
2. Semi-bounded parameters:
A log-normal prior distribution is chosen for all the scaling parameters withsemi-infinite support. These parameters with a lower bound of zero and no upper bound are: • initial values of radii r and vessel wall mechanical properties β at the inlet of the aorta (see Sections 4.1, 4.1.1,and 4.2); • the decay terms used within all exponential property profiles (see Sections 4.1, and 4.2); • the length of arterial vessels (see Section 4.2.2); and • the terminal boundary Windkessel model parameters (see Section 4.3.2)The prior PDF of the i th such parameter is thus described by the following log-normal distribution:p( θ i ) = 1 θ i σ i √ π exp (cid:32) − (ln ( θ i ) − µ i ) σ i (cid:33) , (17)where µ i and σ i represent the mean and standard deviation of the underlying Normal distribution, i.e. ln ( θ i ) ∼N ( µ i , σ i ). This log-normal distribution is created with large variance, resulting in a weakly informative prior. Forall the scaling parameters µ i is set to 0.5 and and σ i is set to 0.8, respectively. This PDF is shown in the middleplot of Figure 6.
3. Unbounded parameters:
A normal prior distribution is chosen for all the scaling parameters with infinitesupport. The only such parameters are the inlet flow-rate FS coefficients (see Section 4.3.1). The prior PDF ofthe i th FS coefficient with mean µ i and standard deviation σ i is thus described as:p( θ i ) = 1 σ i √ π exp (cid:32) − (cid:18) θ i − µ i σ i (cid:19) (cid:33) . (18)While the FS scaling parameters do have infinite support, they are expected to be within physiological range andhence within several multiplications of the reference FS coefficients. Thus, µ i = 1 and a standard deviation σ i = 2is set. This PDF is shown in the right plot of Figure 6. It is assumed that all the priors are independent of eachother, thus the joint prior PDF is: p( θ ) = N (cid:89) i =1 p( θ i ) , (19)where N = 199 is the total number of scaling parameters. The likelihood corresponds to the literature reported measurements. These measurements can be categorised asfollows 12 . . . . . P r ob a b ilit y Ω ref Daughter/parent ratioscaling term0 Ω ref P r ob a b ilit y . . . . . P r ob a b ilit y Figure 6: The three type of prior distributions used for all arterial network scaling parameters.
Measurement Measurement No. of Measured Measured Sourceslocation type patients mean ( µ dis ) std. ( σ dis )Radial artery Diastolic pressure (mmHg) 71 65 12.84 [38], [39]Systolic pressure (mmHg) 123 22.75Ascending aorta Diastolic pressure (mmHg) 69 65 7.14 [38], [40]Systolic pressure (mmHg) 103 8.24Common carotid Diastolic pressure (mmHg) 134 75.58 6.01 [41], [42],artery Systolic pressure (mmHg) 122.78 12.01 [43]Femoral artery Average flow-rate ( ml/ sec) 63 5.84 2.11 [44] [45](left and right)Table 2: The measurement of the discrete pressure and flow-rate taken from literature. When multiple sourcemeasurements are available for the same quantity, they are pooled together into a single mean µ dis and variance σ dis [37]. Pressure (diastolic and systolic) and flow-rate (average) measurements at locations (radial artery, ascending aorta,common carotid, and femoral artery) reported in literature are incorporated into the likelihood. Their statisticsand sources are shown in Table 2. When multiple measurements are available for the same quantity, they arepooled together into a single mean and variance [37]. For such a measurement τ i that was measured (or pooledtogether) to be µ i, dis with standard deviation σ i, dis , the likelihood is:p ( τ i = µ i, dis | θ ) = 1 σ i, dis √ π exp (cid:32) − (cid:18) µ i, dis − H i ( M ( θ )) σ i, dis (cid:19) (cid:33) , (20)where H i represents the observation operator (see Section 5) that extracts the τ i component of the model out-put. Denoting all such discrete measurements collectively with the vector τ dis , the combined likelihood with theassumption of independence is p ( τ dis = µ dis | θ ) = (cid:89) i p ( τ i = µ i, dis | θ ) , (21) Some time-varying information about the behaviour of pressure and flow-rates over the cardiac cycle is necessaryto obtain a physiologically realistic posterior. Statistics on any such time-varying behaviour is not reported inliterature to the authors’ knowledge. Thus, a time-varying pseudo-measurement is constructed by combining thereference (ADAN network [14]) time varying inlet flow-rate Q ref ( t ) and a measurement of average cardiac output13 . . . . . . t (s) − F l o w - r a t e ( m l s − ) Mean profile sd envelope Evaluation pointFigure 7: The measured time domain inlet flow-rate profile taken from the ADAN network is shown in solid black.The standard deviation (sd) envelope computed using Equation (23) are shown with dashed lines and the discretetime points at which the likelihood is evaluated is shown in red circles.[46, 47]. This pseudo-measurement Q meas ( t ) is described by a Gaussian Random Field (GRF, also referred to asa Gaussian Process) [48, 49] with the following mean µ pseudo ( t ) and standard deviation σ meas ( t ) µ meas ( t ) = Q ref ( t ) , (22) σ meas ( t ) = (cid:18) Q ref ( t ) − min( Q ref )max( Q ref ) − min( Q ref ) (cid:19) δ, (23)where δ is a positive scaling parameter. The scaling of σ meas ( t ) is to ensure that i) the standard deviation ispositive, and ii) the variance is proportional to the magnitude of the flow rate (i.e. the percentage variance withrespect to the magnitude remains fixed). The covariance between the pseudo-measurements at any two times t i and t j is described through a periodic kernel [48]cov ( Q meas ( t i ) , Q meas ( t j )) = σ meas ( t ) σ meas ( t ) exp (cid:18) − ( π | t i − t j | /T )( υT ) (cid:19) , (24)where T represents the cardiac period, and υ represents the ratio of the correlation length to the cardiac period.Through a trial-and-error analysis, this ratio is set to υ = 0 . δ is tuned to the average cardiac output (CO) measurement: a meanof 98 ml/sec and standard deviation of 35.55 ml/sec [46, 47]. Through a sweep of δ , an optimal value of δ = 59 . Q meas ( t ) is shown in Figure 7 with the mean and ± one standard deviation profiles.The GRF is evaluated at k = 12 time points, shown in Figure 7, for the computation of the likelihood. Denotingthese evaluation time points as { t , t , · · · , t k } , the vector of measurements at these times with τ inflow , the meanGRF vector as µ inflow = [ µ meas ( t ) , µ meas ( t ) , · · · , µ meas ( t k )] computed through Equation (22), and the vector of in-let flow rates produced by the network parameters through Equation (10) as Q inflow = [ Q inflow ( t ) , Q inflow ( t ) , · · · , Q inflow ( t k )],the likelihood for flow-rate can be written as:p ( τ inflow = µ inflow | θ ) = (2 π ) − k det( Σ ) − exp (cid:18) − ξ T Σ − ξ (cid:19) , (25)where ξ = µ inflow − Q inflow , (26)and Σ is the GRF covariance matrix whose i th row and j th column element Σ i,j = cov ( Q meas ( t i ) , Q meas ( t j )) canbe computed from Equation (24). As mentioned in Section 4.2.2, the length of the vessels are parameterised by a single scaling term relative to thereference network. Since statistics of any direct vessel length measurement are not available, it is assumed that14he lengths of arterial vessels are directly proportional to the height of a subject. A study of 25,945 twins fromeight countries reported the mean and standard deviation of the height of the full cohort to be 172.0 cm and 9.308cm respectively [50]. Since the reference arterial network has a patient of height 170 cm, the measurement datacorresponds to a mean of µ len = 1 . σ len = 0 . τ len , the likelihood for the scaling parameters θ len is thusp ( τ len = µ len | θ ) = 1 σ len √ π exp (cid:32) − (cid:18) µ len − θ len σ len (cid:19) (cid:33) , (27)Assuming all the measurements to be independent, the combined likelihood in Equation (15) is:p ( τ | θ ) = p ( τ dis = µ dis | θ ) × p ( τ inflow = µ inflow | θ ) × p ( τ len = µ len | θ ) , (28)where the RHS terms can be computed from Equations (21), (25), and (27).In summary, the prior distributions for the 199 network parameters are specified in Section 5.2. These priorswill be modified through the likelihood—specified by 7 scalar measurements of pressure and flow-rate, 1 time-varying flow-rate profile evaluated at 12 time points in the cardiac cycle, and a measurement of vessel lengths—toyield the posterior distribution of the parameters through Equation (15). Random sampling from this posteriorwill result in the virtual patient database. The sampling procedure is presented next. With the prior and the likelihood specified, the posterior distribution is given by Equation (15), and can beevaluated at any given θ up to a normalising constant (the evidence term in the denominator of the equation).Sampling from this analytically intractable posterior is achieved through the Markov-chain Monte Carlo method[51, 52], which is a widely used method to sample from unnormalised distributions. The particular MCMC method employed is the Metropolis-Hastings algorithm [53]. Starting from an initialsample, a chain of samples is built sequentially. Denoting the k th sample as θ ( k ) , a candidate for the next sampleis generated by a proposal distribution which depends only on θ ( k ) . This candidate sample θ ∗ is proposed bysampling from the following multivariate Normal distribution: θ ∗ ∼ N ( θ ( k ) , Σ step ) , (29)where Σ step is the covariance matrix, whose i th row and j th column element Σ step ,ij is given byΣ step ,ij = (cid:40) i = jσ ,i i = j , (30)where σ ,i represents the marginal variance for the i th element of θ ∗ . This proposed candidate θ ∗ is eitheraccepted or rejected based on the following ratio α of posterior probabilities at θ ∗ and θ ( k ) α = p ( θ = θ ∗ | τ )p( θ = θ ( k ) | τ ) (31)which, through the Bayes’ rule of Equation (15) applied to both the numerator and the denominator, can bewritten as α = p ( τ | θ = θ ∗ ) p ( θ = θ ∗ )p (cid:0) τ | θ = θ ( k ) (cid:1) p (cid:0) θ = θ ( k ) (cid:1) . (32)The first term and the second term in both the numerator and the denominator are the likelihoods and priors atthe candidate and current point, respectively, and can be computed by Equations (28) and (19). To accept orreject the proposed candidate α is compared to a random number γ drawn from a uniform distribution between 0and 1. If α ≥ γ , the proposed candidate, θ ∗ is accepted and the chain progress with θ ( k +1) = θ ∗ . Otherwise, theproposed candidate is rejected and a new candidate is proposed through Equation (29). This process is repeateduntil the chain is sufficiently long. After discarding the first few thousands of samples from this chain, referred asthe burn-in period [54], the remaining samples are independent samples from the posterior, and form the virtualpatient database.The acceptance rate, i.e. the proportion of candidate sample points that are accepted, is tuned by varying σ ,i . A high variance results in the Markov chain exploring the posterior density quickly, with a high proportionof the candidate points being rejected. A low variance results in slow movement of the chain moving but with a15 θθ ( ) θθθ ( )( ) θθθ ( )( ) = θθθ ( ) θθθ ( )( ) θθθ ( )( ) = θθθ ( )( ) θθθ ( )( ) θθθ ( )( ) = θθθ ( )( ) = θθθ ( ) AcceptReject
Figure 8: An example of an MCMC decision tree across two iterations of a chain. θ ( i )( j ) represents the j th possibleset (candidates) of the arterial network parameter scaling terms at the i th step of the chain. Green branches denotethe path of the chain taken if a candidate is accepted, while red branches show the path of the chain if a candidateis rejected.high proportion of candidate points being accepted. An optimum acceptance rate to balance these two apposingbehaviours is stated in literature to be 0.234 [55]. It is empirically found that using a standard deviation of 0.0375,equivalent to 3.75% of the reference values, for the scaling terms applied to the inlet flow-rate FS coefficients andthat of 0.025, equivalent to 2.5% of the reference values, for all other scaling terms results in an acceptance rateof approximately 0.2. The MCMC algorithm outlined above is inherently sequential—each subsequent sample depends on the previoussample and hence the chain grows only one sample at a time. Generating a long chain this way thus leads tovery high computational run times, since computation of the likelihood requires running a pulse-wave propagationsimulation at each step. To alleviate this issue and achieve some level of parallelisation, the concept of pre-fetching[56, 57] is used. The central idea is to think of MCMC accept-reject decisions in the form a of a decision tree.Figure 8 shows such a decision tree with a depth of η = 2. This results in 2 η leaf-nodes at the end of the tree. Acareful consideration of the decision tree (see Figure 8) shows that the only unique parameters in the entire decisiontree correspond to those at the leaf-nodes. Thus, computational simulations for all the leaf-node parameters canbe run simultaneously in-parallel. Then, η steps of the MCMC algorithm can be taken together without anycomputational overhead by walking the decision tree. A tree depth of η = 4 is used in this study. The VPD is created by generating samples through the MCMC algorithm. The results for the MCMC methodand the VPD are presented and discussed here.
It is common for the initial iterations of an MCMC chain to be considered as the “burn-in” period, during whichthe chain converges to an equilibrium distribution. This burn-in period is, therefore, discarded from the finalsampled posterior distribution. To estimate this burn-in, trace plots for every parameter at all the iterations ofthe MCMC are plotted. An example of a trace plot is shown for the scaling terms applied to the length of VPsarterial vessels in Figure 9, and all other trace plots are shown in Appendix C. To aid visualisation, these plotsare thinned out by plotting every 100 th iteration of the chain. These trace plots are visually assessed to determineburn-in. Essentially, if there is an initial period where a clear migration of the chain is seen away from the sampleused to start the chain, then the burn-in is selected to be until this initial migration is complete. These trace plotsalso show good movement of the chain around the distribution, referred to as the “mixing” of the chain. Figuresin Appendix C show that most parameters are initialised in a region of high density and hence no net migrationis observed—the samples oscillate around a central mean value. Some interesting behaviours in the trace plots,including features that deviate from such a general/desired behaviour, are:16 . . . . . . . . S ca li ng t e r m Figure 9: MCMC trace plot of the scaling term applied to the length of the arterial vessels at every 100 th iteration.In this case, not net migration of the chain is observed. • The scaling term applied to β at the inlet of the aortic chain are seen to centre around 5. The maximum scalingterm applied is approximately 7.5. • The scaling terms applied to the reference decay term of the left leg chain β profile initially oscillates around avalue of approximately 1 for 25,000 iterations, before migrating and oscillating around a value of approximately7 (see Figure 18 in Appendix). This may suggest that this distribution is multimodal. • The scaling terms applied to the reference decay term of the brachiocephalic trunk β profile oscillates around avalue of 1 for the first 40,000 iterations, before migrating and oscillating around a value of 9 (see Figure 19). Thescaling terms applied to the reference decay term of the brachiocephalic trunk r profile shows no complimentarybehaviour and remains within the region of 0–4 throughout. This may again suggest a multimodal distribution. • A spike is seen in the scaling terms applied to the compliance of the right external carotid Windkessel modelat approximately 55,000 iterations (see Figure 33 in Appendix). This behaviour is not seen within the scalingterms applied to the compliance of the left external carotid Windkessel model which remains within the regionof 0–3 throughout.Even though most of the parameters do not show a clear migration, an initial burn-in period of 10,000 iterationsis chosen as a precautionary measure to minimise any effect the initial sample position may have. Once this burn-inperiod has been removed, the VPD contains 65,000 VPs. Along with the burn-in period, all VPs with negativeaverage flow-rate within any arterial vessel are removed from the VPD. These VPs are removed from the VPDas it is physiologically unlikely for a patient to have negative average flow-rate in any vessel. Of the 65,000 postburn-in VPs, 12,857 are removed due to the presence of negative average flow-rate, reducing the VPD to 52,143.
The posterior distribution is a combination of the prior distribution and the likelihood term within Bayes’ theorem.As the prior distributions are weakly informative or uninformative, the distributions of measurements producedby the VPD should be close to the literature based measurements incorporated through the likelihood. These arenot necessarily identical though, as the prior distribution does correct the posterior to account for geometricaland physiological constraints. Furthermore, since the incorporated measurements are taken from several differentsources, the posterior distribution resolves inconsistencies between such measurements, the physics of pulse-wavepropagation, and the constraints of the prior. A comparison of the VPD pressure and flow-rate distributions to theliterature based measurements (see Table 2 and Section 5.3.1) are shown in Figure 10. A similar comparison forthe vessel length scaling term (see Section 5.3.3) is shown in Figure 12. Finally, the statistics of the time varyinginlet flow-rate profiles in the VPD are compared to the GRF (see Section 5.3.2) in Figure 13.Generally, a good agreement between the scalar pressure and flow-rate measurements are seen in Figure 10.This agreement enforces confidence in the overall approach. However, for the average flow-rates in the left andright femoral arteries, larger than expected differences are observed. This is likely due to the following reasons: • A large inconsistency between the femoral flow-rates measurement and the other pressure measurements. Notethat they are taken from different sources and hence not from the same population. Furthermore, there may be17 . . P r o b a b ili t y (a) Diastolic right radial 50 100 150 200Pressure (mmHg)0 . . P r o b a b ili t y (b) Systolic right radial50 100Pressure (mmHg)0 . . P r o b a b ili t y (c) Diastolic left radial 50 100 150 200Pressure (mmHg)0 . . P r o b a b ili t y (d) Systolic left radial40 60 80Pressure (mmHg)0 . . P r o b a b ili t y (e) Diastolic ascending aorta 80 100 120Pressure (mmHg)0 . . P r o b a b ili t y (f) Systolic ascending aorta60 80Pressure (mmHg)0 . . P r o b a b ili t y (g) Diastolic right comm. carotid 75 100 125 150Pressure (mmHg)0 . . P r o b a b ili t y (h) Systolic right comm. carotid60 80Pressure (mmHg)0 . . P r o b a b ili t y (i) Diastolic left comm. carotid 75 100 125 150Pressure (mmHg)0 . . P r o b a b ili t y (j) Systolic left comm. carotid0 5 10Flow-rate ( mls − )0 . . P r o b a b ili t y (k) Average right femoral II 0 5 10Flow-rate ( mls − )0 . . P r o b a b ili t y (l) Average left femoral IILiterature measurement MCMC distribution Figure 10: Histograms of the MCMC distributions of the pressure and flow-rate measurements at all measurementlocations. The literature based measurements and associated error distribution at each location are overlaid inblack. Diastolic and systolic pressure in the right radial artery shown in (a) and (b) respectively; the diastolic andsystolic pressure in the left radial artery are shown in (c) and (d) respectively; the diastolic and systolic pressure inthe ascending aorta are shown in (e) and (f) respectively; the diastolic and systolic pressure in the right commoncarotid artery are shown in (g) and (h) respectively; the diastolic and systolic pressure in the left common carotidartery are shown in (i) and (j) respectively; and the average flow-rate in the second segments of the right and leftfemoral artery are shown in (k) and (l) respectively. 18 5 10
Flow-rate (ml/sec) . . . P r o b a b ili t y (a) Average right femoral I Flow-rate (ml/sec) . . . P r o b a b ili t y (b) Average left femoral I Literature measurement MCMC distribution
Figure 11: Histograms of the MCMC distributions of the average flow-rate in the first segment of the right femoralartery (c), and the first segment of the left femoral artery (d). The literature based measurement of averageflow-rate in the femoral arteries and the associated error distribution is overlaid in black. . . . P r ob a b ilit y Literature measurementMCMC distribution
Figure 12: Histogram of the vessel length scaling terms assigned to VPs across the created VPD . The literaturebased measurement of vessel length scaling terms is overlaid in black.an inconsistency between the cardiac output measurement used for generating the pseudo-measurement for thetime varying inlet flow-rate measurement (see Section 5.3.2). • The presence of more measurements for pressure as opposed to flow-rates. Since the likelihoods for all themeasurements are weighted equally, it may be possible that the chain is influenced weakly by the few flow-rate measurements. Such an issue can be resolved in future studies by assigning variable weightings to themeasurements. • In the network the femoral arteries are split into the first (I) and second segments (II), shown in Figure 5 by‘u’ and ‘v’, respectively. The precise location at which the literature based femoral flow-rate measurement hasbeen taken is unknown. In this study, it is assumed that this measurement was acquired at the centre of thesecond segment ‘v’. Since the first segment ‘u’ bifurcates into ‘v’ and another segment profunda femoris (Figure5), the flow-rate in ‘v’ is smaller than that in ‘u’. Thus, it is possible that the measured flow-rate was in ‘u’. Acomparison of VPD flow-rate in this vessel ‘u’ against the measurement is presented in Figure 11, and shows asignificantly better agreement.Figure 12 shows a good agreement between the vessel length scaling term in the VPD and the literature basedmeasurement. The mean and standard deviation of this term in the VPD are 0.0551 and 0.9819, respectively,which compare well to those reported in literature with values of 0.0548 and 1.0118, respectively.Figure 13 shows a good agreement between the statistics of the time varying inlet flow-rate profile in the VPDand the pseudo-measurement constructed through the GRF (Section 5.3.2). The mean profile closely follows thatof the GRF throughout the cardiac cycle, with maximum difference of approximately at 10% peak systole. Figure13 also shows ± one standard profiles for both the VPD flow-rate and the GRF pseudo-measurement. Throughoutthe cardiac cycle, the agreement and the trend is well captured in the VPD, with largest errors of approximately10% magnitudes at peak systole. 19 . . . . . . Time (s) − F l o w - r a t e ( m l s − ) Evaluation pointLiterature meanLiterature standard deviationMCMC meanMCMC standard deviation
Figure 13: Comparison of the empirical distribution of the time varying inlet flow-rate profiles between the MCMCsamples and the pseudo-measurement created through a GRF. All standard deviation curves in dashed lines depictmean ± Random samples from the VPD are assessed to gain further insights on the VPs and the behaviour of pressure andflow-rate profiles. Pressure profiles are examined at the ascending aorta; right and left radial arteries; and rightand left common carotid arteries. Flow-rate profiles are examined at the right and left second femoral segments. 15VPs are randomly drawn from the VPD (excluding the burn-in), and the pressure and flow-rate profiles associatedwith each are shown at all examination locations in Appendix D. Five VPs of interest are extracted from AppendixD and shown in Figure 14.Figures in Appendix D show that pressure profiles in proximal vessels, i.e. the ascending aorta and two commoniliacs, show more consistency and similarity to the reference profiles compared to the pressure and flow-rate profilesin distal vessels, i.e. the radial artery pressure and the femoral artery flow-rate. This suggests that the pseudo-measurement of the time varying inlet flow-rate profile is sufficient to impose sufficient control over the shape ofpressure and flow-rate profiles in the proximal vessels. As the spatial distance from the inlet of the aorta increases,an increases in the variability of pressure and flow-rate profiles is observed. This is expected as the measurementsthat guide the posterior are only available at few locations. Thus, access to more measurements, distributed acrossthe network in both space and time, will result in an even more realistic database as the likelihood at several suchlocations will guide the posterior.Figure 14 shows a selection of undesirable behaviour in the VPD profiles. In the case of Patient-E, oscillatorybehaviour is seen within all the pressure and flow-rate profiles. This is caused by the shape of the inlet flow-rate profile prescribed for this patient. While this behaviour is not very common in the database, if needed itsoccurrence can be further reduced by imposing stronger correlations between inlet flow-rates at two time pointsthrough the parameter υ in Equation (24) in the GRF.Highly oscillatory pressure profiles are observed in the left radial artery, shown in Figure 14(c), of Patient-Aand Patient-E. In the case of VP Patient-A there is no significant oscillatory behaviour in any other profile, butthe pressure waveform in the right radial artery is featureless with no clear systolic and diastolic points. Thisapparent over- and under-damping of the pressure profiles in the radial arteries may suggest an imbalance in thecompliances or resistances of the right and left arms. This hypothesis of left-right imbalance is further supportedby Figures 14(f) and (g) when the femoral flow-rates are observed for Patient-I: the flow-rate in the right femoralartery shows a high-mean high-oscillation behaviour while in the left femoral artery it is low-mean with significantlylower oscillations. Similar behaviour is seen within Patient-J and Patient-O. In future studies, a symmetry metricthat balances, while still allowing for some variability, the left and right side parameters for symmetric left andright side vessels may be implemented. In this study, it is chosen to filter out vessels that show a high left-rightimbalance with an analysis that is described next. 20 .
00 0 .
25 0 .
50 0 .
75 1 . P r e ss u r e ( mm H g ) (a) Ascending aorta Reference networkPatient-APatient-E Patient-IPatient-JPatient-O0 .
00 0 .
25 0 .
50 0 .
75 1 . P r e ss u r e ( mm H g ) (b) Right radial .
00 0 .
25 0 .
50 0 .
75 1 . P r e ss u r e ( mm H g ) (c) Left radial .
00 0 .
25 0 .
50 0 .
75 1 . P r e ss u r e ( mm H g ) (d) Right common carotid .
00 0 .
25 0 .
50 0 .
75 1 . P r e ss u r e ( mm H g ) (e) Left common carotid .
00 0 .
25 0 .
50 0 .
75 1 . F l o w - r a t e ( m l s − ) (f) Right femoral .
00 0 .
25 0 .
50 0 .
75 1 . F l o w - r a t e ( m l s − ) (g) Left femoral Figure 14: Pressure profiles in the ascending aorta (a), pressure profiles in the right radial artery (b), pressureprofiles in the left radial artery (c), pressure profiles in the right common carotid artery (d), pressure profiles inthe left common carotid artery (e), flow-rate profiles in the right second femoral artery (f), and flow-rate profiles inthe left second femoral artery (g). In each figure the profiles taken from the reference network are shown in black;and the literature reported measurements and associated variances are shown by the solid (mean) and dashed greylines (mean ± The left-right lower extremity imbalance is assessed by reducing the network (see Section 3.3 and Figure 3) up tothe second bifurcation in the leg vessels, i.e up to the end of vessels labelled ‘u’ in Figure 5. To assess the left-rightimbalance, the ratio of femoral pulse (maximum flow-rate minus the minimum flow-rate) on the left and right sidesis considered. It is found that this ratio shows highest correlation against the ratio of reduced network complianceson the left and right sides. A plot of the femoral pulse ratio against the compliance ratio is shown in Figure 15.The left-right imbalance is apparent in this figure which shows that the ratio of left-to-right femoral pulses variesfrom 10 − and 10 across the VPD. To limit this, a filter on the left-to-right compliance ratio between 0.2–5 isintroduced. These limits are shown as vertical lines in Figure 15 and constrains the the femoral ratio between 10 − and 10 across the dataset. With this filter, approximately 45% (23,275 patients) of the patients are discarded,leaving with 28,868 physiologically realistic patients in the VPD.This VPD reflects both prior beliefs/constraints and measurements from the population. Thus, it is suitableto study population characteristics and behaviour. More importantly, disease can be virtually created in thesepatients and subsequent assessment of machine learning be made on the database. The pre-filter VPD is publiclymade available at Zenodo repository . Link to Zenodo repository: https://doi.org/10.5281/zenodo.4549764 − Post second bifurcationleg compliance ratio10 − − R a t i oo f t h e f e m o r a l fl o w - r a t e pu l s e Figure 15: Left-to-right femoral pulse ratio against the left-to-right compliance ratio of reduced network (up tovessel ‘u’ in Figure 5). Filter range on the compliance ratio is shown with vertical lines.
A physiologically realistic virtual patient database is presented for the human arterial network. A methodologyto create virtual patients guided by prior beliefs, geometrical/physiological constraints, and literature reportedmeasurements is presented. Starting from a reference network describing the arterial network, the methodologyincludes i) network reduction without compromising relevant behaviour; ii) re-parameterisation to reduce dimen-sionality; iii) incorporation of geometrical and physiological constraints in the form of a prior; iv) incorporationof literature reported clinical measurements in the form of the likelihood; v) combination of the prior and likeli-hood to generate the posterior; and vi) sampling from the posterior with MCMC to create the VPD. This genericmethodology, given a mathematical description of a biological system, can be adopted to create virtual patientsfor any biological system while accounting for all available information.The key conclusion of this study is that despite computational challenges and dimensional complexity, theBayesian approach to creating virtual patients is viable and yields a realistic database of patients consistent withthe expected behaviour in terms of the range of outputs expected in a population. The study also shows thatappropriate incorporation of priors and measurements (distributed both spatially and in time) is necessary forthe generation of physiologically realistic virtual patients. In the absence of these, the virtual patients may bephysiologically unrealistic while still statistically reproducing the literature measurements. In this study thisphenomenon was observed as no likelihood or prior term incorporating left-right symmetry was included. Whilesuch issues can be fixed by incorporation of a filter, the end result is wasted effort in the generation of virtualpatients which are subsequently discarded. Nevertheless, the Bayesian approach and the framework presented isan attractive framework that is capable of incorporating such expectations in the methodology so that such issuesare avoided.A realistic database of 28,868 patients representing human arterial network is presented and made available.This database is created with a view that it is useful in subsequent studies which exploit the advances in machinelearning methods for the detection of arterial disease. This forms the primary area of future investigation withthis database.
Funding
This work is supported by an EPSRC studentship ref. EP/N509553/1 and an EPSRC grant ref. EP/R010811/1.22 eferences [1] Fowkes, F. G. R., Rudan, D., Rudan, I., et al. “Comparison of global estimates of prevalence and risk factors for peripheralartery disease in 2000 and 2010: a systematic review and analysis”. In:
The Lancet
Journal of the American College of Cardiology
Cerebrovascular diseases
PLoS One
Journal of vascular surgery
American Journal of Physiology-Heart and Circulatory Physiology
Biomechanics and Modeling in Mechanobiology (2020), pp. 1–17.[8] Jones, G., Parr, J., Nithiarasu, P., et al.
A proof of concept study for machine learning application to stenosis detection . 2021.arXiv: .[9] Celik, N., O’Brien, F., Brennan, S., et al. “Deep-Channel uses deep neural networks to detect single-molecule events frompatch-clamp data”. In:
Communications Biology
Biomechanics and modeling in mechanobiology
IEEE Transactions on biomedical engineering
Briefings inBioinformatics . Virtual PiE Led t/a BHR Group: Lisbon, Portugal. 2012, pp. 401–442.[14] Boileau, E., Nithiarasu, P., Blanco, P. J., et al. “A benchmark study of numerical schemes for one-dimensional arterial bloodflow modelling”. In:
International journal for numerical methods in biomedical engineering
Journal of engineeringmathematics
Annals of biomedical engineering
AmericanJournal of Physiology-Heart and Circulatory Physiology
Journal of biomechanics
Medical & biological engineering & computing
International journal for numerical methods in biomedical engineering
Notices of the AMS
Blood pressure monitoring
American journal of hypertension
Clinical Research in Cardiology
Journal of hypertension
Blood pressure monitoring
Scandinavian journal of clinical and laboratory investigation
Journal of medical ultrasound
Journal of Applied Physiology
Archives of Neurology
New England Journal of Medicine
33] Aboyans, V., Desormais, I., Lacroix, P., et al. “The general prognosis of patients with peripheral arterial disease differs accordingto the disease localization”. In:
Journal of the American College of Cardiology
Journal of the American Heart Association
American Journal of Physiology-Heart and Circulatory Physiology
Bayesian theory . Vol. 405. John Wiley & Sons, 2009.[37] Rudmin, J. W. “Calculating the exact pooled variance”. In: arXiv preprint arXiv:1007.1012 (2010).[38] Pauca, A. L., Wallenhaupt, S. L., Kon, N. D., et al. “Does radial artery pressure accurately reflect aortic pressure?” In:
Chest
Circulation
Circulation
Heart and vessels
Pfl¨uger’s Archiv f¨ur die gesamte Physiologie des Menschen und der Tiere
Ultrasound in medicine & biology
Hepatology
Ultrasound in medicine & biology
Heart
The American journal of cardiology
Advanced lectures on machine learning . Springer, 2004, pp. 63–71.[49] Ibragimov, I. A. and Rozanov, Y. A.
Gaussian random processes . Vol. 9. Springer Science & Business Media, 2012.[50] Silventoinen, K., Sammalisto, S., Perola, M., et al. “Heritability of adult body height: a comparative study of twin cohorts ineight countries”. In:
Twin Research and Human Genetics
Statistical science (1992), pp. 473–483.[52] Gilks, W. R. “M arkov Chain M onte C arlo”. In:
Encyclopedia of biostatistics
The american statistician
Markov chain Monte Carlo in practice (1996), pp. 115–130.[55] Roberts, G. O., Gelman, A., Gilks, W. R., et al. “Weak convergence and optimal scaling of random walk Metropolis algorithms”.In:
The annals of applied probability
Journal of Computational and GraphicalStatistics
Computational Statistics& Data Analysis List of vessels within the reference network
The 77 arterial segments in the reference network are listed in Table 3.
Reference number Vessel name Reference number Vessel name
Table 3: The 56 arterial vessels, described by 77 segments, within the reference arterial network, taken from [14],are outlined above. The numbers assigned to each vessel within this table align with those in Figure 2.25
Arterial vessel length
To assess the importance of vessel length on the variability of pressure and flow-rate profiles, two groups of VPsare created using the reference network (see Section 3.3). For both groups, all properties of each VPs arterialnetwork are set to their reference values, except for the length of each vessel. The lengths of VPs arterial vesselsare randomised for each group by: • Group 1: applying an individual independent scaling term to the reference length of each vessel. Thesescaling terms are sampled from normal distributions with mean of 1 and a standard deviation of 0.2. • Group 2: applying a single scaling term to the reference length of all vessels. As with the first case, thisscaling term is sampled from a normal distribution with mean of 1 and standard deviation of 0.2.For each of the two groups outlined above 30,000 VPs are sampled, and the pressure and flow-rate profilesassociated with each computed. Pressure or flow rate profiles are taken from each VP at all measurement locations,highlighted within Section 3.2, and the average, maximum, and minimum of each profile recorded. The average,maximum, and minimum pressure is also recorded at the inlet of the aorta, to compare the affect of vessel length onpressure profiles at a location with known flow rate. The mean and standard deviation of the average, maximum,and minimum pressure at each appropriate examination location is shown for patients with a constant lengthscaling term in Figure 16, and individual scaling terms in Figure 17.Comparing Figure 16 and Figure 17 it can be seen that as is expected, the mean of all pressure measurementsare relatively consistent for both groups of VPs. The standard deviation of all pressure measurements is seen tobe larger for group 2, relative to group 1. This increase in standard deviation is most noticeable when comparingthe maximum and minimum pressure values. The standard deviation of the average pressure at each location isrelatively low for both groups of VPs. While there is an increase in the standard deviation of pressure measurementswhen using a constant scaling term, relative to individual scaling terms, the extent of this increase is reasonablewhen weighed up against the reduction in dimensionality achieved by using a single scaling term.Further analysis is carried out by comparing the mean and standard deviations of each flow rate measurement;and the correlation between all pressure and flow rate measurement for each of the two groups of VPs. Similarbehaviour is seen to the comparison of pressure measurements highlighted above.26 o r ti c B r ac h i a l L B r ac h i a l R C a r o ti d L C a r o ti d R R a d i a l L R a d i a l R P r e ss u r e ( mm H g ) Average ± sdMax ± sdMin ± sd Figure 16: The mean and standard deviation of the average, maximum, and minimum pressure at all appropriateexamination locations for the 30,000 VPs created using a constant length scaling term. A o r ti c B r ac h i a l L B r ac h i a l R C a r o ti d L C a r o ti d R R a d i a l L R a d i a l R P r e ss u r e ( mm H g ) Average ± sdMax ± sdMin ± sd Figure 17: The mean and standard deviation of the average, maximum, and minimum pressure at all appropriateexamination locations for the 30,000 VPs created using individual length scaling terms applied to each vessel withinthe network. 27
MCMC trace plots
The trace plots of all parameters within each VPs arterial network are split into the following figures: • The scaling terms applied to the reference parameters of the vessel wall mechanical property profiles of vesselswith varying properties along their lengths are shown in Figures 18, 19, and 20. • The scaling terms applied to the reference parameters of the vessel radius property profiles of vessels withvarying properties along their lengths are shown in Figures 21, 22, and 23. • The scaling terms applied to the reference daughter-parent ratios of the vessel wall mechanical propertiesand the radii of all vessels bifurcating off of the aorta with constant properties are shown in Figures 24 and28 respectively. • The scaling terms applied to the reference daughter-parent ratios of the vessel wall mechanical propertiesand the radii are shown for vessels with constant properties in the right upper extremities in Figures 25 and29 respectively, and left upper extremities in Figures 26 and 30 respectively. • The scaling terms applied to the reference daughter-parent ratios of the vessel wall mechanical propertiesand the radii of vessels with constant properties within the lower extremities are shown in Figures 27 and 31respectively. • The scaling terms applied to the reference parameters of the Windkessel models in the aortic region, rightupper extremities, left upper extremities, and the legs are shown in Figures 32, 33, 34, and 35 respectively.To aid in visualisation of the results shown in the above listed figures, all trace plots are thinned out by plottingonly each 100 th iteration of the chain. This is done to filter out high frequency noise, clarifying the low frequencybehaviour of the chain. D Pressure and flow-rate profiles from random VPs
Within the subsequent three figures (Figures 36, 37, and 38) pressure and flow-rate profiles are shown from 15 VPsrandomly sampled from the VPD. The pressure or flow-rate profile at each location measured within the referencenetwork is also included, as well as the literature based measurements and associated error incorporated into theposterior distribution. 28 . . . S ca li ng t e r m (a) S ca li ng t e r m (b) . . . S ca li ng t e r m (c) S ca li ng t e r m (d) S ca li ng t e r m (e) . . S ca li ng t e r m (f) S ca li ng t e r m (g) . . S ca li ng t e r m (h) . . . S ca li ng t e r m (i) S ca li ng t e r m (j) Figure 18: The parameter scaling terms at each 100 th iteration of the Markov chain applied to the: aorta chaininitialising value (a), aorta chain decay term (b), right arm chain daughter/parent ratio (c), right arm chain decayterm (d), left arm chain daughter/parent ratio (e), left arm chain decay term (f), right leg chain daughter/parentratio (g), right leg chain decay term (h), left leg chain daughter/parent ratio (i), and the left leg chain decay term(j) for the β properties of the network. 29 . . . . S ca li ng t e r m (a) S ca li ng t e r m (b) . . . S ca li ng t e r m (c) . . . S ca li ng t e r m (d) S ca li ng t e r m (e) S ca li ng t e r m (f) Figure 19: The parameter scaling terms at each 100 th iteration of the Markov chain applied to the: brachiocephalictrunk daughter/parent ratio (a), brachiocephalic trunk decay term (b), right common carotid daughter/parent ratio(c), right common carotid decay term (d), left common carotid daughter/parent ratio (e), and the left commoncarotid decay term (f) for the β properties of the network.30 S ca li ng t e r m (g) S ca li ng t e r m (h) . . . S ca li ng t e r m (i) S ca li ng t e r m (j) . . . S ca li ng t e r m (k) S ca li ng t e r m (l) Figure 20: The parameter scaling terms at each 100 th iteration of the Markov chain applied to the: celiac trunkdaughter/parent ratio (g), celiac trunk decay term (h), right common iliac daughter/parent ratio (i), right commoniliac decay term (j), left common iliac daughter/parent ratio (k), and the left common iliac decay term (l) for the β properties of the network. 31 S ca li ng t e r m (a) . . S ca li ng t e r m (b) . . S ca li ng t e r m (c) S ca li ng t e r m (d) S ca li ng t e r m (e) S ca li ng t e r m (f) . . S ca li ng t e r m (g) . . S ca li ng t e r m (h) . . S ca li ng t e r m (i) . . S ca li ng t e r m (j) Figure 21: The parameter scaling terms at each 100 th iteration of the Markov chain applied to the: aorta chaininitialising value (a), aorta chain decay term (b), right arm chain daughter/parent ratio (c), right arm chain decayterm (d), left arm chain daughter/parent ratio (e), left arm chain decay term (f), right leg chain daughter/parentratio (g), right leg chain decay term (h), left leg chain daughter/parent ratio (i), and the left leg chain decay term(j) for the r properties of the network. 32 S ca li ng t e r m (a) S ca li ng t e r m (b) . . . S ca li ng t e r m (c) S ca li ng t e r m (d) S ca li ng t e r m (e) . . . S ca li ng t e r m (f) Figure 22: The parameter scaling terms at each 100 th iteration of the Markov chain applied to the: brachiocephalictrunk daughter/parent ratio (a), brachiocephalic trunk decay term (b), right common carotid daughter/parent ratio(c), right common carotid decay term (d), left common carotid daughter˙parent ratio (e), the left common carotiddecay term (f) for the r properties of the network. 33 S ca li ng t e r m (g) S ca li ng t e r m (h) . . S ca li ng t e r m (i) S ca li ng t e r m (j) . . . S ca li ng t e r m (k) S ca li ng t e r m (l) Figure 23: The parameter scaling terms at each 100 th iteration of the Markov chain applied to the: celiac trunkdaughter/parent ratio (g), celiac trunk decay term (h), right common iliac daughter/parent ratio (i), right commoniliac decay term (j), left common iliac daughter/parent ratio (k), and the left common iliac decay term (l) for the r properties of the network. 34 S ca li ng t e r m (a) S ca li ng t e r m (b) S ca li ng t e r m (c) S ca li ng t e r m (d) S ca li ng t e r m (e) S ca li ng t e r m (f) S ca li ng t e r m (g) S ca li ng t e r m (h) Figure 24: The parameter scaling terms at each 100 th iteration of the Markov chain applied to the daughter/parentratio of the: inferior mesenteric (a), right renal (b), left renal (c), superior mesenteric (d), second left posteriorintercostal (e), second right posterior intercostal (f), first left posterior intercostal (g), first right posterior intercostal(h) for the β properties of the network. 35 . . . S ca li ng t e r m (a) . . . S ca li ng t e r m (b) . . . S ca li ng t e r m (c) S ca li ng t e r m (d) . . . S ca li ng t e r m (e) . . . S ca li ng t e r m (f) Figure 25: The parameter scaling terms at each 100 th iteration of the Markov chain applied to the daughter/parentratio of the: first right ulnar (a), right common interosseous (b), right radial (c), right vertebral (d), right externalcarotid (e), right internal carotid (f) for the β properties of the network.36 . . . S ca li ng t e r m (a) . . . S ca li ng t e r m (b) . . . S ca li ng t e r m (c) S ca li ng t e r m (d) . . . S ca li ng t e r m (e) . . . S ca li ng t e r m (f) Figure 26: The parameter scaling terms at each 100 th iteration of the Markov chain applied to the daughter/parentratio of the: first left ulnar (a), left common interosseous (b), left radial (c), left vertebral (d), left external carotid(e), left internal carotid (f) for the β properties of the network.37 S ca li ng t e r m (a) S ca li ng t e r m (b) S ca li ng t e r m (c) S ca li ng t e r m (d) S ca li ng t e r m (e) S ca li ng t e r m (f) S ca li ng t e r m (g) S ca li ng t e r m (h) Figure 27: The parameter scaling terms at each 100 th iteration of the Markov chain applied to the daughter/parentratio of the: right anterior tibial (a), left anterior tibial (b), right posterior tibial (c), left posterior tibial (d), rightprofunda femoris (e), left profunda femoris (f), right internal carotid (g), left internal carotid (h) for the β propertiesof the network. 38 S ca li ng t e r m (a) S ca li ng t e r m (b) S ca li ng t e r m (c) S ca li ng t e r m (d) S ca li ng t e r m (e) S ca li ng t e r m (f) S ca li ng t e r m (g) S ca li ng t e r m (h) Figure 28: The parameter scaling terms at each 100 th iteration of the Markov chain applied to the daughter/parentratio of the: inferior mesenteric (a), right renal (b), left renal (c), superior mesenteric (d), second left posteriorintercostal (e), second right posterior intercostal (f), first left posterior intercostal (g), first right posterior intercostal(h) for the r properties of the network. 39 . . . S ca li ng t e r m (a) . . . . S ca li ng t e r m (b) . . S ca li ng t e r m (c) S ca li ng t e r m (d) S ca li ng t e r m (e) . . . S ca li ng t e r m (f) Figure 29: The parameter scaling terms at each 100 th iteration of the Markov chain applied to the daughter/parentratio of the: first right ulnar (a), right common interosseous (b), right radial (c), right vertebral (d), right externalcarotid (e), right internal carotid (f) for the r properties of the network.40 . . . S ca li ng t e r m (a) S ca li ng t e r m (b) . . S ca li ng t e r m (c) S ca li ng t e r m (d) S ca li ng t e r m (e) . . . S ca li ng t e r m (f) Figure 30: The parameter scaling terms at each 100 th iteration of the Markov chain applied to the daughter/parentratio of the: first left ulnar (a), left common interosseous (b), left radial (c), left vertebral (d), left external carotid(e), left internal carotid (f) for the r properties of the network.41 S ca li ng t e r m (a) S ca li ng t e r m (b) S ca li ng t e r m (c) S ca li ng t e r m (d) S ca li ng t e r m (e) S ca li ng t e r m (f) S ca li ng t e r m (g) S ca li ng t e r m (h) Figure 31: The parameter scaling terms at each 100 th iteration of the Markov chain applied to the daughter/parentratio of the: right anterior tibial (a), left anterior tibial (b), right posterior tibial (c), left posterior tibial (d), rightprofunda femoris (e), left profunda femoris (f), right internal carotid (g), left internal carotid (h) for the r properties of the network. 42 S ca li ng t e r m (a) S ca li ng t e r m (b) S ca li ng t e r m (c) S ca li ng t e r m (d) S ca li ng t e r m (e) S ca li ng t e r m (f) S ca li ng t e r m (g) S ca li ng t e r m (h) R R C Figure 32: The parameter scaling terms at each 100 th iteration of the Markov chain applied to the Windkesselmodel parameters at the terminal boundary of the: inferior mesenteric (a), right renal (b), left renal (c), supe-rior mesenteric (d), second left posterior intercostal (e), second right posterior intercostal (f), first left posteriorintercostal (g), and the first right posterior intercostal (h).43 S ca li ng t e r m (a) S ca li ng t e r m (b) S ca li ng t e r m (c) S ca li ng t e r m (d) S ca li ng t e r m (e) S ca li ng t e r m (f) R R C Figure 33: The parameter scaling terms at each 100 th iteration of the Markov chain applied to the Windkesselmodel parameters at the terminal boundary of the: second right ulnar (a), right common interosseous (b), rightradial (c), right vertebral (d), right external carotid (e), and the right internal carotid (f).44 S ca li ng t e r m (a) S ca li ng t e r m (b) . . . S ca li ng t e r m (c) S ca li ng t e r m (d) . . . S ca li ng t e r m (e) S ca li ng t e r m (f) R R C Figure 34: The parameter scaling terms at each 100 th iteration of the Markov chain applied to the Windkesselmodel parameters at the terminal boundary of the: second left ulnar (a), left common interosseous (b), left radial(c), left vertebral (d), left external carotid (e), and the left internal carotid (f).45 S ca li ng t e r m (a) S ca li ng t e r m (b) S ca li ng t e r m (c) S ca li ng t e r m (d) S ca li ng t e r m (e) S ca li ng t e r m (f) S ca li ng t e r m (g) . . S ca li ng t e r m (h) R R C Figure 35: The parameter scaling terms at each 100 th iteration of the Markov chain applied to the Windkesselmodel parameters at the terminal boundary of the: right anterior tibial (a), left anterior tibial (b), right posteriortibial (c), left posterior tibial (d), right profunda femoris (e), left profunda femoris (f), right internal iliac (g), leftinternal iliac (h). 46 .
00 0 .
25 0 .
50 0 .
75 1 . P r e ss u r e ( mm H g ) (a) Ascending aorta Reference networkPatient-APatient-B Patient-CPatient-DPatient-E0 .
00 0 .
25 0 .
50 0 .
75 1 . P r e ss u r e ( mm H g ) (b) Right radial .
00 0 .
25 0 .
50 0 .
75 1 . P r e ss u r e ( mm H g ) (c) Left radial .
00 0 .
25 0 .
50 0 .
75 1 . P r e ss u r e ( mm H g ) (d) Right common carotid .
00 0 .
25 0 .
50 0 .
75 1 . P r e ss u r e ( mm H g ) (e) Left common carotid .
00 0 .
25 0 .
50 0 .
75 1 . F l o w - r a t e ( m l s − ) (f) Right femoral .
00 0 .
25 0 .
50 0 .
75 1 . F l o w - r a t e ( m l s − ) (g) Left femoral(g) Left femoral
75 1 . F l o w - r a t e ( m l s − ) (g) Left femoral(g) Left femoral Figure 36: Examples of pressure and flow-rate profiles taken from five VPs randomly drawn from the VPD. Thesubplots show the: pressure profiles in the ascending aorta (a), pressure profiles in the right radial artery (b),pressure profiles in the left radial artery (c), pressure profiles in the right common carotid artery (d), pressureprofiles in the left common carotid artery (e), flow-rate profiles in the right second femoral artery (f), and flow-rateprofiles in the left second femoral artery (g). In each figure the profiles taken from the reference network are shownin black; and the literature reported measurements and associated error are shown by the solid and dashed greylines respectively. 47 .
00 0 .
25 0 .
50 0 .
75 1 . P r e ss u r e ( mm H g ) (a) Ascending aorta Reference networkPatient-FPatient-G Patient-HPatient-IPatient-J0 .
00 0 .
25 0 .
50 0 .
75 1 . P r e ss u r e ( mm H g ) (b) Right radial .
00 0 .
25 0 .
50 0 .
75 1 . P r e ss u r e ( mm H g ) (c) Left radial .
00 0 .
25 0 .
50 0 .
75 1 . P r e ss u r e ( mm H g ) (d) Right common carotid .
00 0 .
25 0 .
50 0 .
75 1 . P r e ss u r e ( mm H g ) (e) Left common carotid .
00 0 .
25 0 .
50 0 .
75 1 . F l o w - r a t e ( m l s − ) (f) Right femoral .
00 0 .
25 0 .
50 0 .
75 1 . F l o w - r a t e ( m l s − ) (g) Left femoral(g) Left femoral
75 1 . F l o w - r a t e ( m l s − ) (g) Left femoral(g) Left femoral Figure 37: Examples of pressure and flow-rate profiles taken from five VPs randomly drawn from the VPD. Thesubplots show the: pressure profiles in the ascending aorta (a), pressure profiles in the right radial artery (b),pressure profiles in the left radial artery (c), pressure profiles in the right common carotid artery (d), pressureprofiles in the left common carotid artery (e), flow-rate profiles in the right second femoral artery (f), and flow-rateprofiles in the left second femoral artery (g). In each figure the profiles taken from the reference network are shownin black; and the literature reported measurements and associated error are shown by the solid and dashed greylines respectively. 48 .
00 0 .
25 0 .
50 0 .
75 1 . P r e ss u r e ( mm H g ) (a) Ascending aorta Reference networkPatient-KPatient-L Patient-MPatient-NPatient-O0 .
00 0 .
25 0 .
50 0 .
75 1 . P r e ss u r e ( mm H g ) (b) Right radial .
00 0 .
25 0 .
50 0 .
75 1 . P r e ss u r e ( mm H g ) (c) Left radial .
00 0 .
25 0 .
50 0 .
75 1 . P r e ss u r e ( mm H g ) (d) Right common carotid .
00 0 .
25 0 .
50 0 .
75 1 . P r e ss u r e ( mm H g ) (e) Left common carotid .
00 0 .
25 0 .
50 0 .
75 1 . F l o w - r a t e ( m l s − ) (f) Right femoral .
00 0 .
25 0 .
50 0 .
75 1 . F l o w - r a t e ( m l s − ) (g) Left femoral(g) Left femoral