Fractional-Order Windkessel Boundary Conditions in a One-Dimensional Blood Flow Model for Fractional Flow Reserve (FFR) Estimation

Abstract

Recent studies have demonstrated the benefits of using fractional derivatives to simulate a blood pressure profile. In this work we propose to combine a one-dimensional model of coronary blood flow with fractional-order Windkessel boundary conditions. This allows us to obtain a greater variety of blood pressure profiles for better model personalization An algorithm of parameter identification is described, which is used to fit the measured mean value of arterial pressure and estimate the fractional flow reserve (FFR) for a given patient. The proposed framework is used to investigate sensitivity of mean blood pressure and fractional flow reserve to fractional order. We demonstrate that the fractional derivative order significantly affects the fractional flow reserve (FFR), which is used as an indicator of stenosis significance.

1. Introduction

Atherosclerotic diseases of coronary vessels are the main reason for myocardial ischemia frequently resulting in disability or death. These diseases are mainly caused by blockages due to an abnormal narrowing in a blood vessel—stenosis [1]. The choice of medical treatment involves evaluation of stenosis significance, which may require invasive measurements. To assess the severity of each stenosis case, clinicians use various hemodynamic indices. The most popular and well-developed index is the fractional flow reserve (FFR), which is a ratio between mean pressure distal (downstream) to stenosis and mean aortic pressure during artificially induced hyperemia [2,3]. Stenoses with values of FFR below 0.8 are considered to be significant and should be surgically treated.

Measuring FFR involves expensive pressure sensors and specialized equipment. Some patients have multiple stenoses with complicated interactions. These problems led to the development of coronary blood flow models capable of estimating FFR from coronary computed tomography angiography (CCTA) and patient’s data (age, heart rate, stroke volume, blood pressure, etc.). Some of these models are based on solving three-dimensional Navier–Stokes equations [4], but in this work, we concentrate on one-dimensional (1D) models of blood flow [5,6,7,8]. The 1D approach is less time-consuming, and it was shown that 3D and 1D FFR calculations demonstrate similar results [9].

One-dimensional models with fractional derivatives can provide a compromise between accuracy and computational cost. Fractional-order models extend the concept of differentiability and incorporate nonlocal and memory effects using a small amount of parameters. This feature can be used to describe complex flows over different space and time scales without splitting the problem into smaller subproblems. Fractional-order models have been proposed in hemodynamic applications. Some examples are the model of blood flow in viscoelastic vessels [10,11] and the model of heat and mass transfer through an arterial segment, which takes into account the interaction with a magnetic field [12]. Fractional derivatives are used to obtain a more realistic prediction of pulse wave forms, which involves the development of parameter identification procedures [13,14]. The latter task has become especially relevant in recent years with the development of computer technology and increasing interest in inverse problems.

In order to extract the patient-specific data of coronary arteries, CCTA images are used. However, there is no geometry data for arteries of systemic circles beyond coronary vessels which can be accounted for in the model. One approach to resolving this issue is to simulate the whole systemic circle with some averaged parameters [8]. This is a physiologically based approach, and it requires a lot of computational resources and includes many parameters that are difficult to estimate. Another option is to impose pressure-derived boundary conditions directly on the inlet of coronary arteries [9], which represents the impact of smaller arteries and microcirculation. This approach simplifies the model but makes it difficult to investigate the effect of various heart conditions on coronary blood flow. Another approach is to take into account the impact of smaller arteries and microcirculation using a submodel coupled with the blood flow model in coronary arteries. A popular choice is Windkessel-type models which are based on the representation of blood vessels (or the whole systemic circulation) as elastic reservoirs with resistance [5]. This leads to small amount (from 2 to 4) of parameters to describe the influence of systemic circulation. An alternative option, similar to the one we use in this work, is to include a part of the aorta in the model and impose a boundary condition on the end of ascending aorta [15]. This approach allows us to use cardiac output as an inlet boundary condition and calculate pressure in the aorta.

Boundary conditions in blood flow models usually imitate the impact of smaller arteries and microcirculation. A porous media-based approach was previously used to simulate microcirculation [16], and fractional derivatives were used to describe flow in porous media [17]. Fractional derivatives were also used to simulate blood flow in capillary vessels [18].

We describe the flow in the systemic circle and microcirculation using Windkessel-type boundary conditions that utilize fractional derivatives. Fractional derivatives have already been used in Windkessel models to simulate hypertensive and normal blood pressure profiles [14]. We propose to couple the Windkessel fractional derivative model with a 1D coronary blood flow model to obtain a greater variety of aortic pressure profiles. We demonstrate that the resulting shape of the aortic pressure profile allows for better personalization of the model and affects the calculated FFR as well as the patient’s diagnosis.

2. Materials and Methods

2.1. Coronary Blood Flow Model

We simulated coronary blood flow and calculated the FFR with a 1D hemodynamic model [19,20]. This model is based on the flow of incompressible viscous fluid through a network of one-dimensional elastic tubes. The conditions for mass and momentum conservation within the network are expressed as a system of hyperbolic equations for each tube:

$$\frac{\partial A}{\partial t}+\frac{\partial Au}{\partial x}=0,$$

(1)

$$\frac{\partial u}{\partial t}+\frac{\partial }{\partial x}\left(\frac{u^2}{2}+\frac{P}{\rho }\right)=-8\pi \mu \frac{u}{A},$$

(2)

where t is time, x is the coordinate along the vessel (tube), $A=A(x,t)$ is the cross-sectional area, $u=u(x,t)$ is the velocity averaged over the cross-section, $P=P(x,t)$ is the blood pressure, $\rho =1.06$ g/cm${}^3$ is the blood density, and $\mu =4cP$ is the blood viscosity. The right-hand side of Equation (2) represents friction force. An additional relation between the blood pressure and cross-sectional area of the vessel wall is required to close the system:

$$P\left(A\right)={\rho }_w{c}^{2}f\left(A\right),f\left(A\right)=\left\{\begin{array}{c}exp\left(\frac{A}{A_0}-1\right)-1,\frac{A}{A_0}\ge 1\hfill \\ ln\frac{A}{A_0},\frac{A}{A_0}<1,\hfill \end{array}\right.$$

(3)

where ${\rho }_w=1.1$ g/cm${}^3$ is the blood vessel wall density, $A_0$ is the cross-sectional area of the unstressed vessel, and c is elasticity index. The physiological meaning of c is the pulse wave velocity or velocity of small disturbances propagated in the vessel wall [21]. Equation (3) is an analytical approximation of the pressure–area curves obtained in experimental studies [22].

The computational domain consists of the aortic root, aorta, left coronary artery (with branches), and right coronary artery (with branches). The diameters, lengths, and topology of vessels can be extracted from CCTA scans. A simplified version of arterial network is presented in Figure 1. We simulated stenosis as a separate segment with decreased diameter.

One-dimensional vessels are connected to each other in junction points to create an arterial structure. The conditions of mass conservation and total pressure continuity are imposed at the junction points:

$$\sum _i{Q}_i=0,$$

(4)

$$\frac{\rho u_i^2}{2}+P_i=\frac{\rho u_j^2}{2}+P_j,i\ne j.$$

(5)

Equation (4) represents an algebraic sum of influxes and effluxes, where i is the index of a vessel connected to a junction. $u_i$, $u_j$ and $P_i$, and $P_j$ in (5) are the velocities and blood pressures of vessels with indices i and j near the junction point.

On the inlet of the aorta (segment 1 on Figure 1), we set cardiac the output as a periodic time function $Q\left(t\right)$ (Figure 2). The shape of the cardiac output was proposed in [23]. It can be adjusted according to the patient’s heart rate (HR), stroke volume (SV), or other data (peak-to-mean flow ratio, QT-interval).

On the outlet of coronary arteries (segments 4,7,8 on Figure 1), we impose a pressure drop condition:

$$P_k-P_{out}=R_k{Q}_k,$$

(6)

where k is the index of a segment, $P_k$ is the blood pressure at the boundary point, $Q_k$ is the blood flux at the boundary point, and $R_k$ is the hydraulic resistance, $P_{out}$ is the outflow pressure. Outflow pressure can be described as a value of blood pressure at which the microcirculation between arteries and veins stops. It ranges between 20 and 60 mm Hg [5] and can be adjusted according to the patient’s data. Resistances $R_k$ are distributed according to empiric Murray’s law through an algorithm described in [24]. Resistances increase during the systolic phase to simulate contractions of myocardium tissue that hinder coronary blood flow [19].

The boundary condition on the outlet of the aorta (segment 2 on Figure 1) differs from boundary condition on the terminal coronary arteries (6) since the former one represents the whole systemic circle as well as microcirculation. In order to describe the behavior of the microcirculation vessels, the two-element Windkessel model [25] is extensively used:

$$Q_a\left(t\right)=\frac{P_a\left(t\right)-P_{out}}{R_a}+C\frac{d{P}_a}{dt},$$

(7)

where $Q_a$ and $P_a$ are the blood flow and pressure in the aorta, and $R_a$ is the hydraulic resistance of the systemic circle and microcirculation. In (7), compliance C is introduced. which represents the ability of blood vessels to distend and store blood volume. Larger values of C correspond to greater vessel elasticity. Compliance C can be adjusted according to patient’s systolic and diastolic blood pressures, and resistance $R_a$ can be derived from the systemic vascular resistance—the ratio between mean blood pressure and cardiac output [24]. The elastic index c from (3), on the other hand, increases with the rigidity of the vessels and, thus, has a different physiological interpretation.

We solve the hyperbolic system (1) and (2) inside each vessel numerically with the help of an explicit grid-characteristic method [26], which is monotone and first-order accurate. Compatibility conditions imposed on junctions with Equations (4) and (5) and boundary points with conditions (6) and (7) form the system of nonlinear equations which is solved with the Newton method. Compatibility conditions are discretized implicitly with a first-order approximation. Discretizations and convergence studies are presented in [27].

The described model can be used to calculate FFR at any point of the coronary arteries. We calculate FFR as a ratio between mean pressure in the coronary artery distal to stenosis (${\overline{P}}_{dist}^h$) and mean aortic pressure (${\overline{P}}_{aortic}^h$) during vasodilation of the coronary vessels (hyperemia) [2]:

$$FFR=\frac{{\overline{P}}_{dist}^h}{{\overline{P}}_{aortic}^h}.$$

(8)

Hyperemia is simulated by decreasing terminal resistances $R_k$ by 70% [28]. FFR values below 0.8 are considered to be significant. This means that the coronary vessel has an abnormal narrowing (stenosis) that should be surgically treated.

By adjusting C and $P_{out}$ in (6) and (7), we can reproduce patient’s systolic and diastolic blood pressures. However, the range of possible blood profiles is quite limited. To improve the mathematical model, additional elements can be introduced into the arterial structure and the Windkessel model [8]. Instead, in this paper, we propose to use the fractional time derivative in Equation (7).

2.2. Fractional-Order Boundary Conditions

We impose boundary condition on the terminal end of the aorta using the fractional-order Windkessel model, which can be written as follows:

$$Q_a\left(t\right)=\frac{P_a\left(t\right)-P_{out}}{R_a}+C_{\alpha }{D}_t^{\alpha }{P}_a\left(t\right),$$

(9)

where $D_t^{\alpha }$ is a fractional differentiation operator; $\alpha$ is a fractional differentiation order, which is assumed to be between 0 and 2 in this work; and $C_{\alpha }$ is a pseudo-compliance (pseudo-capacitance). Fractional differentiation order $\alpha$ determines the relative degree of interaction between the capacitance of the microvasculature vessels, elastic compliance of the vessels, and the dissipation forces inside them. This, in turn, defines the physiological meaning of $C_{\alpha }$.

Windkessel models with fractional derivatives were extensively studied before [13,29]. However, combining a one-dimensional hemodynamic model with the fractional-order Windkessel boundary condition is a new approach that allows us to represent a greater variety of storage and dissipation effects that can be represented with a single additional parameter $\alpha$.

A large number of different definitions have been proposed for the fractional differentiation operator $D_t^{\alpha }$ [30]. We use the Caputo fractional derivative in this work:

$$D_t^{\alpha }P\left(t\right)=\frac{1}{\Gamma (\lceil \alpha \rceil -\alpha )}\underset{0}{\overset{t}{\int }}\frac{P^{(\lceil \alpha \rceil )}\left(t^{\prime }\right)}{{(t-t^{\prime })}^{1+\alpha -\lceil \alpha \rceil }}d{t}^{\prime },$$

(10)

where $\Gamma$ is a gamma function, $\lceil \alpha \rceil$ is a ceiling of $\alpha$ (smallest integer greater than $\alpha$), and $P^{(\lceil \alpha \rceil )}\left(t^{\prime }\right)$ is a derivative with an integer order $\lceil \alpha \rceil$. There are many other definitions of the fractional derivative: the Atangana–Baleanu fractional integral [31], Riemann–Liouville fractional derivative [32], Riesz derivative [33], etc. The choice of the Caputo fractional derivative is due to its simplicity in representation (relative to other fractional derivatives) and availability of well-studied numerical methods with approximation estimates for various problems. Another useful representation derived in [34] can be obtained using integration by parts in (10):

$$D_t^{\alpha }P\left(t\right)=\frac{1}{\Gamma (-\alpha )}\underset{0}{\overset{t}{\int }}\frac{P\left(t^{\prime }\right)}{{(t-t^{\prime })}^{1+\alpha }}d{t}^{\prime }.$$

(11)

This representation is valid for $0<\alpha <2$. We use it to approximate the integral in (11) with a trapezoidal rule. For the interval $[0,t]$ with a grid $\{t^n=n\tau :n=0,1,2,..,N\}$, assuming a constant time-step $\tau$ and $P\left(0\right)=0,P^{\prime }\left(0\right)=0$, numerical approximation of $D_{\tau }^{\alpha }{P}_N$ can be expressed as

$$D_{\tau }^{\alpha }{P}_N=\frac{{\tau }^{-\alpha }}{\Gamma (2-\alpha )}\sum _{n=0}^N{a}_{n,N}{P}_{N-n},$$

(12)

where coefficients $a_{n,N}$ are

$$a_{n,N}=\left\{\begin{array}{c}1,n=0,\hfill \\ {(n+1)}^{1-\alpha }-2n^{1-\alpha }+{(n-1)}^{1-\alpha },0<n<N,\hfill \\ (1-\alpha )N^{-\alpha }-N^{1-\alpha }+{(N-1)}^{1-\alpha },n=N.\hfill \end{array}\right.$$

The error and stability analysis of this numerical approximation was described in detail in [30,35]. It was demonstrated that the order of approximation error is $O\left({\tau }^{2-\alpha }\right)$.

We use an explicit approximation scheme for (9) with the discretization of fractional derivative described above. This allows us to determine the blood pressure $P_N$ at the terminal end of the aorta from the values of pressure and flux at the previous time steps. Then, we calculate the outflow $Q_N$ at the terminal end of the aorta using compatibility conditions (1) and (2) and wall-state Equation (3).

2.3. Model Personalization

One of the most important problems in patient-specific blood flow modeling is to identify the parameters of the model. Diameters, lengths, and overall arterial structure can be extracted from CCTA images with the help of segmentation and skeletonization algorithms [36]. Parameter c in (3) represents the pulse wave velocity and can be estimated from the patient’s age, blood pressure, heart rate, and stroke volume with the help of machine learning methods [20]. Terminal resistances $R_k,R_a$ in (6), (7) and (9) are calculated from systemic vascular resistance (the ratio between mean pressure and cardiac output) and the diameters of the terminal arteries [24].

Outflow pressure $P_{out}$ and compliance C in (7) are estimated to reproduce measured systolic and diastolic blood pressure. A number of algorithms for $P_{out}$ and C estimation are presented in [25]. The following procedure to estimate these parameters is used in this work. We calculate initial value of C as a ratio between stroke volume and pulse pressure ($PP=P_{sys}-P_{dia}$) and the initial value of $P_{out}$ as 50% of diastolic pressure. After this, C and $P_{out}$ are iteratively adjusted until the measured diastolic and systolic blood pressures match:

$$P_{out}^{i+1}=P_{out}^i\frac{P_{sys}^{true}+P_{dia}^{true}}{P_{sys}^i+P_{dia}^i},C^{i+1}=C^i\frac{P_{sys}^{true}-P_{dia}^{true}}{P_{sys}^i-P_{dia}^i}.$$

(13)

Adjustment of $C_{\alpha }$ in (9) is performed using the same procedure, but initial estimation is usually less precise since the dimension and interpretation of $C_{\alpha }$ changes with $\alpha$.

Unfortunately, relying solely on systolic and diastolic pressures may produce incorrect diagnostic outcomes. For example, the hemodynamic significance of stenosis may vary for the same values of systolic and diastolic pressures. In order to obtain more accurate estimates, additional available information about the pressure profile, such as mean pressure, must be taken into account.

We propose to estimate an additional parameter, fractional derivative order $\alpha$, based on the value of the mean pressure $P_{mean}$. To do this, we first iteratively adjust $\alpha$ to match the mean pressure and then adjust $C_{\alpha }$ and $P_{out}$ for each $\alpha$ to match the systolic and diastolic blood pressures. After this, we calculate $P_{mean}$ and compare it with the measured mean pressure. If the calculated mean pressure is higher than the measured one, the order $\alpha$ should be decreased, and vice-versa. As we will see from the results, the relationship between $\alpha$ and $P_{mean}$ is very close to linear. Therefore, in most situations, it is sufficient to perform two preliminary calculations for $\alpha =1$ and $\alpha =1.5$ or $\alpha =0.5$ to estimate $\alpha$, which provides the appropriate value for the mean pressure.

The procedure of parameter identification described above is summarized on Figure 3.

2.4. Patient Data

We tested the parameter identification procedure (Figure 3) on a publicly available dataset [5]. This dataset contains data from ten patients, including the geometry of arterial networks and the location stenoses in various coronary vessels. We kept the numeration of patients from [5] but we excluded Patient 9 from our study since the measurement of mean pressure is unavailable. As a result, our study included nine patients (

Table 1. Characteristics of the patient dataset (mean ± standard deviation). Details are presented in [5]. $\theta =\frac{P_{mean}-P_{dia}}{P_{sys}-P_{dia}}$ is a measure of blood profile thickness.

Characteristic	Value
Number of patients	9
Number of males	5
Heart rate, bpm	69 ± 14
Systolic pressure $P_{sys}$, mm Hg	141 ± 23
Diastolic pressure $P_{dia}$, mm Hg	73 ± 8
Mean pressure $P_{mean}$, mm Hg	103 ± 11
BMI, kg/m${}^2$	29 ± 4
$\theta$	0.45 ± 0.10

with 13 stenoses).

Figure 4 presents examples of two patient-specific structures for Patient 1 and Patient 4.

Table 2. Characteristics of Patient 1 and Patient 4. Stroke volume (SV) was estimated from patient age and BMI values presented in [5].

	Patient 1	Patient 4
Age, years	80	68
HR, bpm	67	88
SV, mL	82	70
$P_{sys}$, mm Hg	174	130
$P_{dia}$, mm Hg	76	66
$P_{mean}$, mm Hg	111	94
Stenosis location	LAD	LAD
Stenosis degree	70%	60%
FFR measured	0.89	0.82

presents patient metadata as well as measured FFR values. These two patients were chosen based on the value $\theta$, which is the measure of the blood profile thickness:

$$\theta =\frac{P_{mean}-P_{dia}}{P_{sys}-P_{dia}}.$$

(14)

The average value of $\theta$ across all nine patients is ${\theta }_{ave}=0.45$. Patient 1 has the lowest value ${\theta }_1=0.36$, and Patient 4 has the closest to the average value ${\theta }_4=0.44$. As a result, Patient 1 has the thinnest blood pressure profile and Patient 4 has the most typical profile.

3. Results

3.1. Blood Pressure and FFR Sensitivity to Order $\alpha$

The proposed model was applied to calculate the blood pressure profiles for various fractional differentiation orders $\alpha$ in a simplified network of coronary arteries (Figure 1). We also calculated the FFR for 66% of the stenoses for various $\alpha$.

First, we calculated aortic blood pressure for $\alpha =1.0$ and adjusted the model parameters to acquire the physiological systolic (125 mm Hg) and diastolic (75 mm Hg) blood pressures. Then, we performed calculations for other values of $\alpha$ with the same set of model parameters. The value of compliance C remains constant for various $\alpha$. This is technically incorrect since the physiological meaning and dimensional formula for C depend on $\alpha$, but it helps us to explore changes in pressure profiles with the change of $\alpha$ (Figure 5a).

Second, we adjusted C and $P_{out}$ for each order $\alpha$ to obtain the same values of systolic and diastolic blood pressures (Figure 5b). As $\alpha$ decreases, the pressure peak shifts to the left, and the pressure profile becomes thinner. This results in a drop in mean pressure. Conversely, as $\alpha$ increases, the profile peak shifts to the right, the profile becomes thicker, and the mean pressure grows.

Systolic and diastolic blood pressures are one of the most commonly available patient-specific parameters. All blood profiles on Figure 5b have the same systolic and diastolic blood pressures, but the mean pressures are different for each $\alpha$. If the mean pressure is available for a given patient, we can use it to select an appropriate $\alpha$ and calculate FFR. Figure 6 demonstrates how the calculated mean pressure and FFR depend on $\alpha$, assuming that the systolic and diastolic blood pressures are the same (125/75 mm Hg). The relationship between mean pressure and $\alpha$ is very close to linear within the considered interval (from 0.25 to 1.5). Therefore, this simplifies the process of $\alpha$ identification from patient’s mean pressure: if we calculate the mean pressures for any two values of $\alpha$, we can derive an appropriate $\alpha$ for a given mean pressure value using linear interpolation.

At the same time, FFR drops with increasing $\alpha$ (Figure 6b). This relation resembles exponential decay. For $0<\alpha <1$, FFR drops rapidly from 0.95 ($\alpha =0.25$) to 0.78 ($\alpha =1.0$). The threshold between the significant and insignificant lesions is 0.8, so the choice of $\alpha$ affects the diagnostic outcome. For $1<\alpha <2$, FFR is almost constant.

3.2. Patient-Specific Calculations

In this section, we describe applying a parameter estimation algorithm (Figure 3) to estimate FFR for nine patients (Table 1). We start with two examples: Patient 1 and Patient 4.

Patient 1 had stenosis in the left anterior descending artery (LAD) with a corresponding measured FFR value of 0.89 and a mean aortic pressure of 111 mm Hg. Blood flow calculations for $\alpha =1$ (utilizing boundary condition (7)) yielded $FF{R}^{\alpha =1}=0.83$ and $P_{mean}^{\alpha =1}=119$ mm Hg. Then, we applied the proposed model with fractional-order Windkessel boundary conditions (9) with the parameter identification procedure to fit the given systolic, diastolic, and mean pressures. As a result, the following estimates were obtained: $P_{mean}^{\alpha =0.56}=111$ mm Hg was achieved for $\alpha =0.56$, and the resulting FFR was $FF{R}^{\alpha =0.56}=0.87$. The error in FFR estimation was significantly reduced after applying optimization of $\alpha$ for mean pressure.

Patient 4 had stenosis in the LAD with a corresponding measured FFR value of 0.82 and a mean aortic pressure of 94 mm Hg. Blood flow calculations for $\alpha =1$ (utilizing boundary condition (7)) yielded $FF{R}^{\alpha =1}=0.82$ and $P_{mean}^{\alpha =1}=94$ mm Hg. No further optimization was required in this case since the calculated mean pressure matched the measured mean pressure with good accuracy. The calculated FFR value also matched the measured one.

FFR estimations for other patients are presented in

Table 3. FFR estimations for the patient dataset. Patient data: $P_{mean}$ is the measured mean pressure, mm Hg; $\theta$ is a measure of the pressure profile thickness (14); Loc. is a location of stenosis; $FFR$ is the invasively measured FFR. The original approach for FFR estimation (order $\alpha =1.0$): $FF{R}_{\alpha =1}$ is the calculated FFR value with a boundary condition (7); $P_{mean}^{est}$ is the calculated mean pressure with order $\alpha =1$. The fractional derivative approach (order $\alpha ={\alpha }_{opt}$) involves adjusting order $\alpha$ so that the calculated mean pressure matches the measured one: $FF{R}_{\alpha ={\alpha }_{opt}}$ is the calculated FFR value with the boundary condition (9), and fractional order $\alpha ={\alpha }_{opt}$; ${\alpha }_{opt}$ is the optimal fractional order. Patient 9 was excluded due to the absence of a mean pressure measurement.

	Patient Data				Order $\mathit{\alpha }=1.0$		Order $\mathit{\alpha }={\mathit{\alpha }}_{\mathit{opt}}$
№	${\mathit{P}}_{\mathit{mean}}$	$\mathit{\theta }$	Loc.	$\mathit{FFR}$	${\mathit{FFR}}_{\mathit{\alpha }=1}$	${\mathit{P}}_{\mathit{mean}}^{\mathit{est}}$	${\mathrm{FFR}}_{{\alpha }_{\mathrm{opt}}}$	${\alpha }_{\mathrm{opt}}$
1	111	0.36	LAD	0.89	0.83	118	0.87	0.56
2	83	0.46	LAD	0.86	0.87	83	0.87	1.0
3	125	0.40	RCA	0.88	0.89	125	0.89	1.0
4	94	0.44	LAD	0.82	0.82	94	0.82	1.0
5	99	0.39	LAD	0.82	0.8	102	0.82	0.82
6	99	0.40	LADp	0.9	0.98	101	0.98	0.91
			LADd	0.82	0.87		0.88
			DA	0.81	0.84		0.85
7	98	0.37	LAD	0.75	0.68	102	0.71	0.78
			LCx	0.84	0.82		0.85
8	110	0.51	LAD	0.88	0.92	107	0.92	1.3
			LCx	0.89	0.97		0.96
10	108	0.51	LAD	0.72	0.81	90	0.8	1.85

. The original approach to estimate FFR involves a boundary condition (7) without fractional derivative. We ignored the measured mean pressure and adjusted our model to achieve the measured systolic and diastolic blood pressures. The fractional-order approach involves adjusting the fractional derivative order to obtain measured mean pressure. The RMSE for the original approach was 0.05, and the RMSE for the fractional order approach was 0.04. The RMSE was mainly defined by large errors in the FFR estimations of patients 6, 8, and 10. We assumed that stenosis degree and length were not identified properly for these patients. The FFR estimation was improved for patients with a “thin” blood profile ($\theta <0.4$), including Patient 1, Patient 5, and Patient 7. Patients 2, 3, 4, and 6 had am optimal fractional order ${\alpha }_{opt}$ close to 1.0 (or equal to 1.0), so FFR estimations for both approaches were similar. The FFR estimations for patient 6 were less precise for the fractional order approach, but the difference was very small. Patients 8 and 10 had an optimal fractional order ${\alpha }_{opt}>1.0$, and the FFR estimation was similar for both approaches. This is due to the fact that FFR is almost constant for $\alpha >1.0$ (Figure 6).

Results of the FFR estimation showed that the fractional-order approach provided benefits for a certain group of patients (patients 1, 5, and 7) with a thin blood pressure profile. These patients may have a common cardiovascular pathology. According to [14], low fractional orders can be used to simulate the blood profiles of patients with hypertension. Patients 1, 5, and 7 have a BMI > 30 kg/m${}^2$ which is a good predictor of hypertension. However, their blood pressure levels are not the highest in the dataset. We need a larger sample of patients to make further conclusions.

4. Discussion

We proposed coupling the well-established one-dimensional hemodynamic model of coronary blood flow with a fractional-order boundary condition, as well as a procedure for estimating its parameters. The fractional derivative can be a useful tool for more accurate modeling of the pressure profile. The actual blood pressure profiles depend on many factors, such as age, height, weight, medical history, and artery elasticity. The most commonly available characteristics of blood pressure profiles are systolic and diastolic blood pressures. Unfortunately, relying solely on these two values alone may produce incorrect diagnostic outcomes. For example, the hemodynamic significance of stenosis may vary for the same systolic and diastolic blood pressures. This fact has motivated researchers to look for new tools to model blood pressure profiles.

We used the fractional derivative order to match calculated mean blood pressure with the measured one. Adjusting mean pressure can be performed in other ways: introducing a larger Windkessel model, expanding the arterial network, etc. The fractional-derivative Windkessel model is a good compromise: we introduced a very small amount of new parameters (1–2) and gained the ability to simulate a whole spectrum of dissipative and storage mechanisms with the help of fractional order $\alpha$. Our approach does require additional information on the blood pressure profile. These data can be obtained with simple noninvasive procedures. Unfortunately, in many cases, these data are unavailable, and all the pressure profile information is reduced to systolic and diastolic pressure values. This was the case for Patient 9 who was excluded from our study.

The proposed approach has a number of shortcomings that need to be resolved in the future. First, in some cases, the only data available regarding a patient’s pressure profile were systolic and diastolic blood pressures. Identifying the fractional derivative order $\alpha$ in this case would require a completely different approach that can be based on other data, such as patient medical history. Second, calculating fractional derivatives calls for significant computational resources. This negates one of the main advantages of one-dimensional blood flow models—computational efficiency. Instead of the basic approach presented in this work, new efficient numerical approximations can be used. Third, the proposed approach should be tested on a larger number of patients with a wide range of FFR values. The FFR is nearly constant for fractional differentiation orders $\alpha >1$, so for some patients, its adjustment will be useless.

Further research will focus on integrating more efficient approaches to identify model parameters and to approximate a fractional derivative. Other areas of work include collaborating with clinicians to find effective methods for pulse wave assessment and further validation on a larger amount of patients. The proposed approach has great potential to provide an alternative means to simulating arterial stiffness and pulse waves.