ML-ROM wall shear stress prediction in patient-specific vascular pathologies under a limited clinical training data regime

TL;DR

This paper combines machine learning and reduced order modeling to predict blood flow stress in patient-specific vascular conditions using limited clinical data.

Contribution

A novel integration of POD-based ROM and ML models for WSS prediction in vascular pathologies with limited training data.

Findings

ML-ROM models achieved significant speed-up ratios (10⁴) in predicting hemodynamic indices.

Flowrate-coefficients mapping model outperformed autoregressive models in limited data scenarios.

Accuracy was higher in Peripheral Arterial Disease compared to Aortic Dissection cases.

Abstract

High-fidelity numerical simulations such as Computational Fluid Dynamics (CFD) have been proven effective in analysing haemodynamics, offering insight into many vascular conditions. However, these methods often face challenges of high computational cost and long processing times. Data-driven approaches such as Reduced Order Modeling (ROM) and Machine Learning (ML) are increasingly being explored alongside CFD to advance biomechanical research and application. This study presents an integration of Proper Orthogonal Decomposition (POD)-based ROM with neural network-based ML models to predict Wall Shear Stress (WSS) in patient-specific vascular pathologies. CFD was used to generate WSS data, followed by POD to construct the ROM. The ML models were trained to predict the ROM coefficients from the inlet flowrate waveform, which can be routinely collected in the clinic. Two ML models were explored: a simpler flowrate-coefficients mapping model and a more advanced autoregressive model. Both models were tested against two case studies: flow in Peripheral Arterial Disease (PAD) and flow in Aortic Dissection (AD). Despite the limited training data sets (three flowrate waveforms for the PAD case and two for the AD case), the models were able to predict the haemodynamic indices, with the flowrate-coefficients mapping model outperforming the autoregressive model in both case studies. The accuracy is higher in the PAD case study, with reduced accuracy in the more complex case study of AD. Additionally, the computational cost analysis reveals a significant reduction in computational demands, with speed-up ratios in the order of 104 for both case studies. This approach shows an effective integration of ROM and ML techniques for fast and reliable evaluations of haemodynamic properties that contribute to vascular conditions, setting the stage for clinical translation.

Figures

Diagrammatic overview of the study methodology.(a) Full Order Model: Patient-specific data is used for CFD simulations to obtain WSS. (b) Reduced Order Model: SVD is applied to WSS data to generate POD mode structures and temporal coefficients. (c) Machine learning: a model is trained to predict temporal coefficients from inlet mass flowrate waveforms. The dashed line indicates the supervised training, i.e., the models see the true coefficients during the training phase. (d) Reconstruction: Predicted coefficients reconstruct WSS, enabling calculation of TAWSS and OSI.

ML models.(a) Flowrate-coefficients mapping model and (b) Autoregressive model.

Patient-specific geometry of the femoral artery with boundary conditions,showing a flow split of 33% and 67% at the outlets. The graph below presents mass flowrate waveforms for training (μ1, μ2, μ3) and testing (μ4, μ5) datasets over a cardiac cycle. (Figure modified from Ninno et al. [10] with permission.)

POD-based ROM performance in case study 1: PAD.(Left) Reconstruction error of ROMs retaining the first i modes and (right) their RIC. The selected ROM retains r = 10 modes to achieve a reconstruction error less than 5.00%.

Performance comparison of the flowrate-coefficients mapping and autoregressive models on the test dataset, case study 1: PAD.Each bar shows average NMAE and NRMSE for TAWSS and OSI, with 95% confidence intervals

Comparison of TAWSS in the PAD under μ5:ML prediction (top), original CFD (middle), and their differences (bottom). The detailed view shows a region with a relatively high magnitude of absolute differences.

Comparison of OSI in the PAD under μ5:ML prediction (top), original CFD (middle), and their differences (bottom). The detailed view shows a region with a relatively high magnitude of absolute differences.

Bland-Altman plots for TAWSS (left) and OSI (right) comparing ML predictions and CFD results in the PAD case under μ5.The mean difference and limits of agreement (±1.96 SD) are indicated. To enhance readability, the graph displays a subset of only 2,000 randomly chosen data points.

Patient-specific geometry of the AD with boundary conditions,showing the inlet velocity profile and three-element Windkessel model at each outlet. The graph on the top left corner shows mass flowrate waveforms for training (μ−25%, μ+25%) and testing (μ+0%) datasets over a cardiac cycle. TL and FL are colored blue and grey, respectively. (Figure modified from Stokes et al. [5] with permission).

POD-based ROM performance in case study 2: AD.(Left) Reconstruction error of ROMs retaining the first i modes and (right) their RIC. The selected ROM retains r = 26 modes to achieve a reconstruction error less than 5.00%.

Performance comparison of the flowrate-coefficients mapping and autoregressive models on the test dataset, case study 2: AD.Each bar shows average NMAE and NRMSE for TAWSS and OSI, with 95% confidence intervals.

Comparison of TAWSS in the AD under μ+0%:ML prediction (left), original CFD (middle), and their differences (right).

Comparison of OSI in the AD under μ+0%:ML prediction (left), original CFD (middle), and their differences (right).

Bland-Altman plots(left) and OSI (right) comparing ML predictions and CFD results in the TL and FL of AD under μ+0%. The mean difference and limits of agreement (±1.96 SD) are indicated. To enhance readability, the graph displays a subset of 1,000 randomly chosen data points each from the TL and FL (a total of 2,000 data points displayed).

Equations

Full text

Introduction

Wall Shear Stress (WSS) is a haemodynamic metric that has been found to be closely linked to many Cardiovascular Diseases (CVDs) [1–10]. Accurate evaluation of WSS is thus highly important for understanding disease progression and aiding clinical decision-making. High-fidelity simulations such as Computational Fluid Dynamics (CFD) offer powerful tools to analyse 3D blood flow and obtain WSS in cardiovascular systems, providing valuable insights into disease mechanisms and potential treatment strategies [11]. These tools have significantly contributed to our understanding of vascular flow behaviour across various medical conditions, such as aortic aneurysm [1, 2], aortic dissection (AD) [3–5, 12], peripheral arterial disease (PAD) [6–10], and coronary artery disease [13–15].

Despite its benefits, CFD demands substantial expertise and is highly reliant on proprietary software. Moreover, it often involves a trade-off between accuracy and complexity [16]. Accurate simulations demand high computational costs and time, unsuitable in clinical settings that require rapid decision-making [16, 17]. To overcome this challenge, researchers have increasingly adopted Machine Learning (ML) alongside traditional CFD to push the boundaries of biomechanical research and applications [18–34]. These methods have proven effective in various haemodynamics studies, ranging from predicting blood flow quantities and haemodynamic indices [19–26, 32–34] to enhancing flow data resolution and noise reduction [27–31]. A key advantage of ML models is their ability to leverage complex relationships within large datasets using data-driven approaches [35, 36].

The complexity in most models scales with the data dimensionality, thus reducing these dimensions can decrease computational and memory demands. Moreover, simpler models (with fewer inputs) tend to exhibit less variance against noise and outliers [36]. Dimensionality Reduction (or Model Order Reduction) is a class of data-driven techniques used to transform a high-dimensional Full Order Model (FOM) into a lower-dimensional form, known as Reduced Order Model (ROM) [18, 36, 37]. Proper Orthogonal Decomposition (POD) is among the most widely used methods for this purpose. [31, 39]. POD works by decomposing the FOM into a set of orthogonal modes, capturing the most significant features with minimal loss of information. This can be achieved through Singular Value Decomposition (SVD) [18, 37, 38, 40]. POD has been employed in numerous cardiovascular flow investigations. For instance, Chang et al. [41] used POD to construct computationally efficient ROMs to study the flow and WSS in Abdominal Aortic Aneurysms (AAA) with varied inflow angle; Di Labbio and Kadem [42] compared the use of POD and Dynamic Mode Decomposition (DMD) in identifying coherent flow structures in a left ventricle with aortic regurgitation; Buoso et al. [43] developed a computational approach utilizing a parameterised ROM based on POD to accelerate the calculation of pressure drop along stenotic blood vessels. More recently, Chatpattanasiri et al. [31] explored the use of a variation of POD, called Robust POD (RPOD), to construct computationally efficient ROMs of the velocity field inside an AD.

POD-based ROMs (or PCA-based ROM) can also be integrated with ML predictive models to help simplify the prediction of haemodynamic quantities. This entails two major steps: offline and online. In the offline step, the FOM data is collected through traditional CFD or in vitro experiments, and then processed to construct the ROM through POD. This step also involves training the ML model to predict POD coefficients that represent haemodynamic quantities of interest. In the online step, the trained model is employed to make fast and accurate predictions of those quantities in unseen cases (test cases). Notable examples of this approach include the work by Pajaziti et al. [23] who used PCA and Feed-forward Neural Networks (FFNNs) to predict velocity and pressure fields in different aorta geometries. Drakoulas et al. [25] developed a model referred to as FastSVD-ML-ROM which utilized an SVD update methodology and a Convolutional Autoencoder for dimensionality reduction. Their approach also involved FNNs and a Long short-term memory (LSTM) network for predicting the ROM coefficients (sometimes referred to in their study as latent variables or temporal scales of the reduced representations). Siena et al. [24] combined a POD-based ROM with FFNNs to predict time-dependent velocity, pressure, and WSS in coronary artery bypass grafts with varying levels of stenosis. More recently, MacRaild et al. [32] used POD and different designs of FFNNs to predict the time-dependent velocity magnitude in an internal carotid artery with aneurysm, and Gerdroodbary and Salavatidezfouli [34] used POD along with LSTM to predict time-dependent haemodynamic quantities in patient-specific carotid arteries with and without stenosis. Beyond POD, other dimensionality reduction techniques have also been integrated with ML predictive models. For instance, Liang et al. [20] predicted velocity and pressure fields for different aortic shapes using Autoencoders and FNNs.

This work focuses on developing ML models to predict WSS from input quantities that are commonly measured in the clinic such as flowrate waveforms, with the model trained on highly limited datasets typically available in such environments. The methodology involves high-fidelity CFD simulations to generate WSS data, followed by the application of POD to construct the ROM. The ML models are trained to predict the ROM coefficients from the inlet mass flowrate waveforms. The predicted coefficients can then be converted to the 3D WSS data and its related haemodynamic indices. This approach is demonstrated through two case studies: Case study 1: PAD and Case study 2: AD. The former serves as a simple case study with predominantly unidirectional and laminar flow, while the latter represents a more complex case involving flow splitting into two channels: the true lumen (TL) and false lumen (FL). This introduces more intricate flow patterns and turbulent flow regimes. Both case studies involve very limited training datasets, using only three flowrate waveforms for the PAD case and two for the AD case, considering that past studies of comparable complexity typically used 10 or more conditions in the training dataset [25, 44, 45]. It is crucial to highight that this is the reality of routinely acquired clinical datasets, which is often at odds with research requirements. Motivated by this limitation, we aim to achieve high accuracy and robustness with a simpler and more interpretable ML model that provides fast and reliable WSS predictions, enhancing the potential for clinical applications in cardiovascular disease diagnosis and treatment planning.

Methods

Fig 1 illustrates a diagrammatic overview of the study methodology, divided into four major phases. The first phase involves the high-fidelity modelling of vascular haemodynamics. CFD is employed to simulate the flow fields inside the blood vessel of interest with multiple flowrate waveforms. The time-dependent WSS field is calculated and used as the FOM. More details can be found in Computational fluid dynamics and FOM. The second phase focuses on the construction of the ROM, where POD is applied to the WSS data to extract eigenmodes and the corresponding temporal coefficients as detailed in ROM via POD. The third phase involves the development of the ML model designed to predict the ROM coefficients from the mass flowrate waveform explained in Machine learning predictive model. Lastly, the predicted coefficients are used to reconstruct the predicted WSS using the ROM, and the haemodynamic indices: Time-average WSS (TAWSS) and Oscillatory Shear Index (OSI), are calculated, as detailed in ROM via POD.

Diagrammatic overview of the study methodology.(a) Full Order Model: Patient-specific data is used for CFD simulations to obtain WSS. (b) Reduced Order Model: SVD is applied to WSS data to generate POD mode structures and temporal coefficients. (c) Machine learning: a model is trained to predict temporal coefficients from inlet mass flowrate waveforms. The dashed line indicates the supervised training, i.e., the models see the true coefficients during the training phase. (d) Reconstruction: Predicted coefficients reconstruct WSS, enabling calculation of TAWSS and OSI.

Computational fluid dynamics and FOM

In our study, the FOM was derived from CFD simulations. Since blood is an incompressible fluid, its motion can be described by the Navier-Stokes (NS) and continuity equations given below:

$$[eqn]$$ $$[eqn]$$

in the domain $[eqn]$ , where $[eqn]$ and $[eqn]$ are the unknown velocity and pressure fields, with x representing the position vector in 3D coordinates, t representing time, T is the period of the cardiac cycle, and $[eqn]$ is a set of controlled physical parameters (in this study, it is the inlet mass flowrate waveform). $[eqn]$ is the fluid density, $[eqn]$ is the shear stress tensor, and f is the body force per unit volume (e.g., gravity). Appropriate boundary conditions are applied at the domain boundaries ( $[eqn]$ ) to enforce the influence of $[eqn]$ . CFD involves solving the NS and continuity equations numerically on the flow domain that has been discretized into a mesh, which can be done with CFD solver packages such as Ansys Fluent or CFX. After u and p are obtained, WSS can be calculated by:

$$[eqn]$$

where $[eqn]$ is the unit vector normal to the vessel walls. This WSS data derived from the CFD results is used as the FOM.

Specific assumptions and configurations (including the mesh, numerical schemes, turbulence model, etc.) for PAD and AD cases are discussed separately in Case study1: PAD and Case study 2: AD, respectively.

ROM via POD

POD decomposes the data into a set of modes where structures are arranged depending on their energy content. The higher energy modes represent the coherent structures in the flow. A detailed description of POD can be found in Berkooz et al. [40] and in the textbook by Brunton and Kutz [37]. Only a brief overview is provided here.

POD is implemented using the method of snapshots. Consider a 3D WSS dataset (i.e. FOM) $[eqn]$ under $[eqn]$ , described on $[eqn]$ by a position vector x. The dataset consists of N spatial positions and Nt temporal snapshots (usually N>>Nt). The instantaneous WSS data is first separated into a time-independent reference value, commonly taken as the mean value $[eqn]$ [46], and the disturbance from the reference $[eqn]$ . In this context, the mean value $[eqn]$ refers to the overall mean across the population of $[eqn]$ . Thus, $[eqn]$ is independent of $[eqn]$ . In this study, the mean of the training dataset is used as the proxy of $[eqn]$ . Then, the disturbance part is further decomposed into a set of spatial structures $[eqn]$ multiplied by temporal coefficients ai as follows:

$$[eqn]$$ $$[eqn]$$

To compute $[eqn]$ , $[eqn]$ is arranged in a matrix format, stacking each point and each vector component into a single column ( $[eqn]$ ), and arranging all the columns together in a $[eqn]$ matrix $[eqn]$ called snapshot matrix of WSS under $[eqn]$ :

$$[eqn]$$

A large snapshot matrix representing WSS data under multiple conditions of $[eqn]$ can then be constructed by concatenating multiple $[eqn]$ together :

$$[eqn]$$

Singular Value Decomposition (SVD) is then applied directly to the large snapshot matrix $[eqn]$ :

$$[eqn]$$

where U and V are the left and the right singular vectors of $[eqn]$ respectively. Each column of U contains the POD mode structure $[eqn]$ . The POD temporal coefficients can then be computed by projecting $[eqn]$ onto $[eqn]$ : $[eqn]$ .

The singular matrix ( $[eqn]$ ) is a diagonal matrix containing the singular values ( $[eqn]$ ) of $[eqn]$ . The singular values rank in descending order, and they indicate the level of contribution of each corresponding POD mode to the overall dynamics. Many complex dynamical systems show a rapid decline in singular values [18, 47], allowing the use of low-dimensional ROMs to approximate the high-fidelity FOM with high accuracy. Additionally, $[eqn]$ can be used to compute the Relative Importance Content (RIC) which is another metric used to quantify the contribution of the retained modes to the overall system dynamics. It is calculated as $[eqn]$ .

The ROM based on POD can be obtained by slightly modifying Eq 3b:

$$[eqn]$$

where r denotes the number of modes included in the ROM. When setting r = N, Eq 7 yields the FOM. $[eqn]$ can be arranged into a reconstructed snapshot matrix $[eqn]$ or $[eqn]$ , and the reconstruction error is then defined as:

$$[eqn]$$ $$[eqn]$$

Eq 8a and Eq 8b are for a single case and multiple cases, respectively.

In the subsequent phase, the temporal coefficients ai were used to develop the ML prediction model. By training the ML model with these coefficients, we enabled it to predict the POD coefficients for unseen conditions in the parameter space. These predicted coefficients ãi can then be substituted into Eq 7 to obtain the estimated WSS and its related indices.

Machine learning predictive model

Two ML models were explored in this study: a Flowrate-coefficients mapping model and an Autoregressive model. The flowrate-coefficients mapping model is a straightforward prediction model that maps flowrate data to output coefficients. The Autoregressive model is a more advanced model that predicts future values based on past data, building on techniques used in multiple previous studies [25, 44, 48].

Flowrate-coefficients mapping model

The flowrate-coefficients mapping model (Fig 2a) takes $[eqn]$ that has been arranged into a window of w time steps $[eqn]$ as the input. The input is then processed through 2 LSTM layers with 200 neurons per layer, and another 2 dense layers with 100 neurons per layer. It then predicts the output of $[eqn]$ at the same set of time steps ( $[eqn]$ ). The window size w is set to be 8 in this study.

ML models.(a) Flowrate-coefficients mapping model and (b) Autoregressive model.

The flowrate-coefficients model was trained in a supervised manner using the Mean Squared Error (MSE) loss function and the Adam optimizer with a learning rate of 10^−3^. A randomly selected $[eqn]$ of the training dataset was reserved as a validation dataset. An early stopping technique was implemented to automatically end the training when the validation MSE stopped improving for 20 consecutive epochs.

Autoregressive model

The autoregressive model (Fig 2b) was designed to advance the prediction of ROM coefficients into future time steps based on multiple past steps. The network autoregressively predicts $[eqn]$ ) using the previous w = 8 steps of coefficients (at $[eqn]$ ). $[eqn]$ , that has been arranged into a window of w time steps $[eqn]$ , is concatenated to the input to improve the generalizability of the model. This approach was also used in Drakoulas et al. [25].

The autoregressive model consists of 2 LSTM layers (200 neurons per layer) and 2 dense layers (100 neurons per layer). The model was trained in a supervised manner with the Adam optimizer, MSE loss function, and early stopping criteria in the same configurations as the flowrate-coefficients model. Since this model relies on the previous time steps of $[eqn]$ to start the prediction, the flowrate-coefficients mapping model was used as the initialiser to predict the first set of $[eqn]$ , providing the starting point for the autoregressive model.

All the mentioned hyperparameters including the number of layers and the number of neurons per layer in both ML models were chosen empirically during the training process. Further optimization using techniques such as Bayesian optimization is possible [49, 50], although this is beyond the scope of this research.

Results

The performance of the proposed ML models was evaluated via two clinical case studies: Case study 1: PAD representing a simpler flow scenario, and Case study 2: AD which involves more intricate flow patterns. In both case studies, the primary aim was to predict $[eqn]$ from $[eqn]$ . Two key haemodynamic indices, TAWSS and OSI, were calculated to assess the accuracy of the predicted WSS using the following equations [5, 10]:

$$[eqn]$$ $$[eqn]$$

The model accuracy is assessed using two metrics: Normalized Mean Absolute Error (NMAE) and Normalized Root Mean Square Error (NRMSE):

$$[eqn]$$ $$[eqn]$$

where $[eqn]$ can be TAWSS or OSI. These error metrics are similar to those used in Liang et al. [20], but here, they are normalized differently to reflect the true magnitude of the error in our study.

3D reconstructions of TAWSS and OSI, and Bland-Altman plots for both quantities were used to evaluate the results. In addition, the 3D reconstructions of WSS at four phases of the cardiac cycle (acceleration, peak systole, deceleration, and diastole) and a plot of the mean absolute error over a cardiac cycle are provided in the Supplementary Material.

Case study 1: PAD

Problem description

PAD is a circulatory condition primarily caused by atherosclerosis where the buildup of fats, cholesterol, and other substances in the arterial walls results in a narrowing of the arterial lumen, reducing blood flow to the limbs. This leads to symptoms ranging from leg pain and numbness to gangrene and ulceration, the latter of which is prone to infection. In severe instances, these symptoms can progress to the point where amputation becomes necessary, significantly affecting quality of life and raising healthcare costs [51, 52]. It often requires interventions such as angioplasty or bypass surgery to restore blood circulation. However, restenosis in PAD may develop over time as the body’s response to the treatment leads to a gradual re-narrowing of the arteries, causing reduced blood circulation and the recurrence of the complications described previously. While the exact cause of restenosis in PAD is still unclear, researchers have identified that WSS-related indices are linked with the risk and progression of the restenosis [7, 9, 10]. Therefore, developing predictive tools for WSS may significantly improve monitoring and treatment strategies for restenosis in PAD.

This case study utilised data from a recent study by Ninno et al. [10] exploring how discrepancies in the timing between Computed Tomography (CT) scans (which facilitate the reconstruction of vessel geometry) and Doppler Ultrasound (DUS) images (which defines inlet flow boundary conditions) affect the assessment of haemodynamic indices in predicting restenosis (Fig 3). This work received ethical approval from West Haven VA Connecticut Healthcare Systems (approval number AD0009). The CFD package Ansys Fluent (Ansys Inc., PA, USA) was used to solve NS and continuity equations describing blood flow in patient-specific femoropopliteal bypasses. The fluid domain was discretised using tetrahedral elements with refined layers near the wall. Blood was modelled as a non-Newtonian fluid with Carreau viscosity and constant density. The flow was assumed as laminar. Transient simulations were conducted for each bypass using inlet velocity waveforms extracted from DUS images. A parabolic profile was imposed at the inlet, and a flow split of 33% to profunda femoral and 67% to bypass was prescribed at the outlets. The vessel wall was assumed rigid with no-slip conditions. Two cardiac cycles were simulated, and the first cycle was excluded to eliminate the influence of initialisation parameters. The WSS data was obtained using Eq 2.

Patient-specific geometry of the femoral artery with boundary conditions,showing a flow split of 33% and 67% at the outlets. The graph below presents mass flowrate waveforms for training (μ1, μ2, μ3) and testing (μ4, μ5) datasets over a cardiac cycle. (Figure modified from Ninno et al. [10] with permission.)

One patient (patient 3 or PT3) from the patient cohort studied in Nino et al. [10] was selected for this study. There were four waveforms for this patient (acquired by DUS at different dates) presented in Ninno et al. [10]. To enrich the dataset for a more comprehensive analysis, an additional simulation was performed using another waveform from the same patient, thus expanding the total number to five waveforms. The first three waveforms ( $[eqn]$ , $[eqn]$ , and $[eqn]$ ) formed the training dataset, while the remaining two ( $[eqn]$ and $[eqn]$ ) were used as the test dataset. With a time step size of 0.005 seconds, the temporal snapshots for each waveform are: 201 for $[eqn]$ , 135 for $[eqn]$ , 156 for $[eqn]$ , 151 for $[eqn]$ , and 201 for $[eqn]$ .

ROM construction

POD was applied to the snapshot matrix of the training dataset $[eqn]$ to extract $[eqn]$ . Multiple ROMs were then created by truncating different numbers of modes (Eq 7). The reconstruction errors and RIC associated with these ROMs were calculated (Eq 8b), and shown in Fig 4.

POD-based ROM performance in case study 1: PAD.(Left) Reconstruction error of ROMs retaining the first i modes and (right) their RIC. The selected ROM retains r = 10 modes to achieve a reconstruction error less than 5.00%.

To achieve a reconstruction error below $[eqn]$ , a ROM with r = 10 modes was selected for further ML model development. Note that this $[eqn]$ reconstruction error was calculated from $[eqn]$ using Eq 8b. For the test cases ( $[eqn]$ and $[eqn]$ ), the reconstruction errors were computed using Eq 8a and found to be $[eqn]$ and $[eqn]$ , respectively. These small reconstruction errors in the test cases indicated that the flow fields inside the training dataset effectively captured the important flow characteristics of the test cases.

ML performance

Fig 5 shows a comparison of the performance of the autoregressive model and the flowrate-coefficients mapping model in predicting TAWSS and OSI. Each model was trained 10 times, and the errors displayed are the average values of NMAE and NRMSE for TAWSS and OSI with 95% confidence intervals.

Performance comparison of the flowrate-coefficients mapping and autoregressive models on the test dataset, case study 1: PAD.Each bar shows average NMAE and NRMSE for TAWSS and OSI, with 95% confidence intervals

The results indicate that the flowrate-coefficients mapping model outperforms the autoregressive model in all error metrics. For the autoregressive model, the $[eqn]$ and $[eqn]$ are 15.27 $[eqn]$ 2.93% and 15.68 $[eqn]$ 2.84%, respectively, while for OSI, the $[eqn]$ and $[eqn]$ are significantly higher at 34.47 $[eqn]$ 10.71% and 30.86 $[eqn]$ 8.41%, respectively. In contrast, the flowrate-coefficients mapping model shows much lower errors, with $[eqn]$ and $[eqn]$ of 6.20 $[eqn]$ 0.80% and 6.37 $[eqn]$ 0.82%, respectively, and for OSI, the $[eqn]$ and $[eqn]$ are at 21.48 $[eqn]$ 4.37% and 22.41 $[eqn]$ 3.04%, respectively.

Among the 10 flowrate-coefficients mapping models trained, the best was used for further qualitative analysis in Figs 6–8 with $[eqn]$ case.

Comparison of TAWSS in the PAD under μ5:ML prediction (top), original CFD (middle), and their differences (bottom). The detailed view shows a region with a relatively high magnitude of absolute differences.

Comparison of OSI in the PAD under μ5:ML prediction (top), original CFD (middle), and their differences (bottom). The detailed view shows a region with a relatively high magnitude of absolute differences.

Bland-Altman plots for TAWSS (left) and OSI (right) comparing ML predictions and CFD results in the PAD case under μ5.The mean difference and limits of agreement (±1.96 SD) are indicated. To enhance readability, the graph displays a subset of only 2,000 randomly chosen data points.

Fig 6 presents a 3D comparison of TAWSS derived from the ML model against the original CFD data. Both the ML prediction and original CFD data exhibit similar TAWSS distribution along the artery, with high values appearing in similar regions. This indicates that the ML model is generally effective in capturing essential flow dynamics, without any consistent trend of under- or over-prediction across the artery. The differences are primarily confined to very small areas around the valves (present in the vein that was used to create this bypass). They increase the artery’s cross-sectional area, likely introducing more flow disturbances, and consequently reducing prediction accuracy. Similarly, the OSI derived from ML predictions closely matches the spatial distribution patterns observed in the original CFD results, as shown in Fig 7. However, the plot reveals more noticeable areas of under- and over-prediction by the ML model. These discrepancies can be attributed to the fact that OSI calculations consider the directional changes and magnitude of WSS over a cycle, and are inherently more complex than TAWSS computations. This complexity can challenge the ML model’s predictive accuracy, as OSI is sensitive to subtle flow dynamics and temporal variations that are more nuanced than the average shear stress measurements. Nevertheless, the regions exhibiting high discrepancy are small compared to the overall artery surface area where the prediction is accurate.

Fig 8 shows the Bland-Altman plots for TAWSS and OSI in the $[eqn]$ case, assessing the agreement between the ML predicted and original CFD-derived values on each node on the vessel wall. Similar to Fig 6, the TAWSS plot shows excellent performance, with an extremely small under-prediction and a mean bias of –0.03% and limits of agreement from –19.68% to –16.92%. For OSI, the mean bias is –17.85%, with wider limits of agreement from –74.80% to 39.11%, highlighting greater variability in OSI values. The model tends to under-predict when OSI values are low, whereas the errors are closer to 0% for higher OSI values. This trend suggests that regions with low WSS fluctuation are estimated to exhibit even less fluctuation. This may be attributed to the construction of the ROM from truncated POD modes. While this approach can effectively capture dominant flow features, it neglects smaller variations in the higher (truncated) modes. This causes inaccuracies in the regions where the flow dynamics are complex but have lower magnitudes of WSS, such as low OSI regions.

Case study 2: AD

Problem description

Type-B AD is a serious vascular condition that can lead to disability or death. It occurs when an intimal tear develops in the wall of the aorta distal to the left subclavian artery. This tear separates the aorta into two distinct blood flow channels: TL and FL. This causes severe pain and frequently leads to organ ischemia such as renal, limb or mesenteric ischemia, and can also progress to aneurysmal degeneration [53, 54]. The complexity of AD is heightened by its patient-specific nature, with significant variability in tear size, location, and progression of TL and FL. This variability affects the blood flow dynamics, making AD modelling significantly more complicated and challenging compared to the PAD case study, which involves more uniform arterial narrowing patterns and streamlined flow.

The data used in this case study was based on the work of Stokes et al. [5] studying the impact of different inlet conditions on key haemodynamic indices involving aneurysmal growth in type-B AD (Fig 9). The dataset was from a 56-year-old male patient diagnosed with chronic Type B AD acquired following an approved ethics protocol (ID 2019-00556, Inselspital, Bern, Switzerland).

Patient-specific geometry of the AD with boundary conditions,showing the inlet velocity profile and three-element Windkessel model at each outlet. The graph on the top left corner shows mass flowrate waveforms for training (μ−25%, μ+25%) and testing (μ+0%) datasets over a cardiac cycle. TL and FL are colored blue and grey, respectively. (Figure modified from Stokes et al. [5] with permission).

The geometry and boundary conditions were obtained from Computed tomography angiography (CTA) and 4D-magnetic resonance imaging (4D-MRI). CFD simulations were conducted using Ansys CFX 2020 (Ansys Inc., PA, USA) to solve the 3D incompressible Unsteady Reynolds-averaged Navier-Stokes (URANS) and continuity equations. Blood was represented as a non-Newtonian fluid following the Carreau–Yasuda model, with simulations assuming rigid wall boundaries. (In a chronic Type B AD the dissected intima is overlayed with a neo-intima resulting in a rather thick flap with little or no motion in comparison to a freshly dissected aortic wall.) The $[eqn]$ SST turbulence model was employed and three-element Windkessels were incorporated at the outlets to simulate peripheral resistance and compliance. The simulation offers a dataset of 128 temporal snapshots, with a time step size of 0.005 s. The WSS data was then obtained using Eq 2.

Stokes et al. [5] presented four simulations with different inlet profiles including a 3D inlet velocity profile (3DIVP, or referred to as +0% case in this study), a flat profile, a through-plane profile, and a condition with a 25% increased flowrate. To broaden the training dataset for this current study, we introduced an additional simulation case, reducing the flowrate by 25%. The flat and through-plane cases were excluded from this study as they share the same inlet mass flowrate waveform as the +0% case. The -25% and +25% cases were used as the training dataset, while the +0% case was selected as the test dataset. These flowrate waveforms are depicted in Fig 9. $[eqn]$ , $[eqn]$ , and $[eqn]$ represents the –25%, +0%, and +25% cases, respectively.

ROM construction

The snapshot matrix for the training dataset $[eqn]$ was constructed and then POD was applied to extract $[eqn]$ . ROMs were then generated by truncating different numbers of modes, and the corresponding reconstruction errors are calculated and displayed in Fig 10. Unlike the PAD case study where only 10 modes were necessary to achieve a reconstruction error below 5%, 26 modes were needed for this case. Subsequently, the ROM with r = 26 modes was tested against the +0% case, resulting in a reconstruction error of 14.56%, significantly higher than that observed in the training dataset. This larger error suggests notable differences in WSS patterns between the training and test datasets, highlighting the complexity of this case study.

POD-based ROM performance in case study 2: AD.(Left) Reconstruction error of ROMs retaining the first i modes and (right) their RIC. The selected ROM retains r = 26 modes to achieve a reconstruction error less than 5.00%.

ML performance

The performance of the ML models was evaluated using the same approach as in the PAD case study. Two models were tested: the flowrate-coefficients mapping model and the autoregressive model. Fig 11 shows the performance comparison between the two models in predicting TAWSS and OSI for the $[eqn]$ (test case). Similar to the PAD case, each model was trained 10 times, and the errors displayed are the average values of NMAE and NRMSE for TAWSS and OSI with 95% confidence intervals.

Performance comparison of the flowrate-coefficients mapping and autoregressive models on the test dataset, case study 2: AD.Each bar shows average NMAE and NRMSE for TAWSS and OSI, with 95% confidence intervals.

The flowrate-coefficients mapping model substantially outperforms the autoregressive model in this case study. The $[eqn]$ and $[eqn]$ for the autoregressive model are 26.34 $[eqn]$ 0.55% and 39.13 $[eqn]$ 0.81%, respectively, while for OSI, the $[eqn]$ and $[eqn]$ are higher at 27.23 $[eqn]$ 0.82% and 30.95 $[eqn]$ 0.34%, respectively. In contrast, the flowrate-coefficients mapping model shows lower errors, with $[eqn]$ and $[eqn]$ at 8.88 $[eqn]$ 0.51% and 14.24 $[eqn]$ 0.66%, respectively, and for OSI, the $[eqn]$ and $[eqn]$ are 26.45 $[eqn]$ 1.81% and 29.15 $[eqn]$ 1.33%, respectively.

Interestingly, while the TAWSS errors from the flowrate-coefficients mapping model are significantly lower than those from the autoregressive model, the OSI errors for both models are relatively close. Overall, the results highlight that the flowrate-coefficients mapping model is more effective and reliable than the autoregressive model in predicting the haemodynamic quantities in AD. This advantage is consistent with findings in the PAD case study, further underscoring the importance of choosing an appropriate model complexity, especially when dealing with limited training data.

Figs 12 and 13 show the comparison between ML predictions from the flowrate-coefficients mapping model and original CFD results for the $[eqn]$ case in terms of TAWSS and OSI, respectively. While both TAWSS and OSI derived from the ML predictions match the CFD-derived values in general, there are notable areas of discrepancy. Specifically, the differences in TAWSS primarily localize near regions of high curvature such as the aortic arch and branch entries. For OSI, variations are more pronounced along the wall of the FL after the tear and distally downstream, where the OSI values are relatively high. Similar to the PAD case study, the discrepancy in OSI is more noticeable than that of TAWSS due to the reasons stated in Case study 1: PAD.

Comparison of TAWSS in the AD under μ+0%:ML prediction (left), original CFD (middle), and their differences (right).

Comparison of OSI in the AD under μ+0%:ML prediction (left), original CFD (middle), and their differences (right).

Fig 14 displays the Bland-Altman plots for TAWSS and OSI in the $[eqn]$ case. The TAWSS plot shows a mean bias of –2.80% with limits of agreement from –27.41% to 21.82%. The OSI plot shows a mean bias of -26.70% with wider limits of agreement from –120.10% to 66.69%. Both plots indicate that the ML model tends to slightly under-predict TAWSS and OSI. The variability is higher for OSI, reflecting the model’s greater challenge in accurately capturing the oscillation of WSS. Unlike the PAD case, the OSI plot here shows no noticeable tendency to under-predict when the OSI is low. The absence of this trend may be due to the presence of larger sources of errors, which obscure the effect of POD truncation mentioned earlier in Case study 1: PAD. There is no noticeable difference between errors in TL and FL.

Bland-Altman plots(left) and OSI (right) comparing ML predictions and CFD results in the TL and FL of AD under μ+0%. The mean difference and limits of agreement (±1.96 SD) are indicated. To enhance readability, the graph displays a subset of 1,000 randomly chosen data points each from the TL and FL (a total of 2,000 data points displayed).

Computational cost

The computational time for the two case studies is presented in Table 1. For case study 1 (PAD), the CFD simulations used the cluster provided by the Department of Computer Science at UCL (Intel Xeon Gold 5118 at 2.3 GHz using 10 processors), taking approximately 36 hours for each simulation. In case study 2 (AD), the simulations ran on an Intel(R) Core(TM) i9-10900X at 3.7 GHz using 10 processors, requiring roughly 3 days per simulation. However, each AD simulation case also required additional manual fine-tuning of Windkessel parameters, thus multiple rounds of simulations were needed before achieving the suitable set of Windkessel parameters, averaging about 3 weeks per simulation case. For this reason, it took as much as about 1,000 hours to generate the training dataset for the AD case study.

Table 1: Computational time.

The ROM and ML tasks were performed on an Intel(R) Core(TM) i9-12900K at 3.2 GHz and Nvidia RTX A2000, respectively. While the ML evaluation time differed in the two case studies because of the difference in the number of time steps in a cycle, these values were very close at 0.030 $[eqn]$ 0.001 s and 0.034 $[eqn]$ 0.004 s per time step. This is because the same ML model architecture with an equal number of layers and neurons was used for both case studies. The only difference is the number of modes in each ROM: more modes are included in Case study 2: AD, thus it took slightly more evaluation time. This showcases the ML model’s capability to handle tasks with varying complexity while maintaining consistent computational demands. The speed-up ratios are approximately 22,000 for Case study 1: PAD and 60,000 for Case study 2: AD. It is crucial to note that these reported speed-up ratios are conservative estimates because different computational devices were used for the training/testing of the ML model and the CFD simulations.

Discussion

The application of POD-based ROM combined with neural network-based ML models showcased different levels of success in the two case studies. As anticipated, the accuracy was high in the simpler PAD case study and decreased in the more complex AD case study. The higher difficulty of the AD case can be attributed to the intricate flow dynamics caused by the complex geometry that separates into TL and FL, along with the presence of a turbulent flow regime. This trend was evident not only in the ML predictive accuracy but also earlier during the ROM construction phase. The ROM for PAD required only 10 modes to reach the reconstruction error of 5% and it generalized well to unseen test cases ( $[eqn]$ and $[eqn]$ ). In contrast, the ROM for AD required 26 modes to achieve similar reconstruction accuracy and showed significantly larger errors when applied to the unseen case ( $[eqn]$ ). Notably, the need for more than twice the number of modes in AD signifies a far more complex WSS landscape. This higher-dimensional ROM basis implies that the ML model was tasked with learning a more intricate coefficient space, which also contributed to the larger prediction errors observed in the AD case. In contrast, PAD’s WSS dynamics were effectively captured with fewer modes, simplifying the learning task and yielding better generalization.

While ML has been widely applied to model haemodynamics in many cardiovascular conditions [19–24, 26], very few past studies have tackled AD haemodynamics modelling using ML. A very recent study by Deneker et al. [55] introduced warm-start physics-informed neural networks (WS-PINNs) to analyze the velocity field inside the FL of Type B AD, showing effectiveness in handling MRI noise. Similar to our findings, their study acknowledged the challenge of accurately predicting complex flow patterns in AD. Their meshfree (point cloud-based) approach offered greater flexibility in dealing with geometrical variations. However, their approach did not utilize dimensionality reduction techniques, leading to considerably longer training times, and each new case required separate retraining (although partially expedited by their transfer learning technique). In contrast, our ROM-based ML models operate using low dimensional representations of the haemodynamic quantity, which significantly decreases computation time while still providing reasonable accuracy and valuable qualitative insights (as evidenced by the good overall agreement of TAWSS and OSI spatial distribution patterns in Fig 12 and Fig 13, respectively).

This makes our approach more practical for clinical applications where fast and reliable predictions are crucial. Such WSS predictions can be valuable in a clinical setting. For a PAD patient, the ability to quickly estimate TAWSS across a reconstructed artery or bypass graft could help identify regions of persistently low WSS – areas prone to restenosis or intimal hyperplasia – and thus support treatment planning or graft design adjustments [6–10]. In the context of AD, being able to predict WSS (with TAWSS and OSI) in near-real-time may aid clinicians’ decisions on monitoring versus intervention, as regions with abnormal shear (e.g., high OSI or low TAWSS zones often linked to aneurysmal growth) can be non-invasively identified [3–5, 12].

In both case studies presented in our work, the simple flowrate-coefficients mapping model outperformed the more advanced autoregressive model in all error metrics. This result may seem unexpected, given the success of autoregressive models in similar prediction tasks in the literature [25, 44, 48]. The key reason for the superior performance of the simpler model in our study was likely due to the highly limited training dataset available in both case studies: 3 and 2 conditions for PAD and AD, respectively. It is important to highlight that this is extremely limited compared to works in the literature that used autoregressive-type models. For example, Drakoulas et al. [25] trained their FastSVD-ML-ROM on a dataset comprising 10 inlet conditions. Maulik et al. [44] used a training dataset consisting of 5 conditions to model the 1D viscous Burgers’ equation, and another dataset with 90 conditions for the 2D inviscid shallow water equations. In the work by Ahmed et al. [48], flow fields with 5 conditions were used to train the model to learn the dynamics of the Marsigli flows, in which a fluid is divided into two sections with different temperatures. When the separating barrier is suddenly removed, the fluids flow over each other, driven by convection and buoyancy forces. These prior studies suggest that autoregressive networks generally require more extensive training data to perform effectively. By nature, more complex neural network models have a higher tendency to overfit, especially with small datasets, whereas simpler models are known to be more robust and generalizable to unseen cases [36]. In particular, the smaller training sets in both PAD and AD gave the autoregressive model insufficient examples to learn time-dependent patterns without overfitting. Meanwhile, the simpler mapping approach effectively leveraged the available data to capture principal relationships between flowrate and WSS, resulting in consistently better error metrics in our study. This type of approach is also far more compatible with the reality of clinical applications where limited datasets are often the norm.

This finding underscores the importance of balancing ML model complexity with the size and quality of the training dataset. In clinical applications, where data acquisition can be challenging and datasets are often limited, simpler models may lead to more reliable predictions. Moreover, simpler models offer greater flexibility for future improvements and expansions. For instance, if additional input features, e.g. geometric parameters of the vessels, are to be incorporated, a simpler model will be easier to modify to accommodate these needs.

Limitations and future work

A key limitation of our current study is the extremely limited training dataset, which restricts the ML model’s effectiveness in capturing the full variability of the time-dependent WSS through its POD coefficients. To address this, future work could involve implementing statistical generative techniques to enlarge the training dataset by creating additional synthetic data. This has been demonstrated in the study of Pajaziti et al. [23] who used Statistical Shape Modelling (SSM) to create 3,000 synthetic aortic geometries from 67 real geometries before using them as the training dataset for ML prediction of steady-state velocity and pressure fields. Similar techniques have also been used in Liang et al. [20] and Du et al. [22]. While SSM is used primarily for analyzing geometries, the core concept of capturing variability in a dataset can be applied to other data types. Tools like PCA can generate synthetic data by perturbing principal component coefficients to create new sets of $[eqn]$ , expanding our training dataset.

Lower performance in complex pathologies (e.g. AD case) is another key limitation of this work. The complex flow dynamics, caused by turbulent effects and severe flow splitting between the true and false lumens, appear to be more difficult to capture by the current workflow. This shortcoming may arise from several factors: the inherent complexity of the haemodynamics in such geometries, the very limited training data available for these complex conditions, and the reliance on traditional linear dimensionality reduction techniques (in this case, POD) that may not adequately represent the finer-scale, non-linear features of the WSS data. To address this limitation, future research could focus on enriching the training datasets specifically for complex pathologies using statistical generative techniques as mentioned above. Additionally, incorporating more advanced non-linear dimensionality reduction methods such as autoencoders may better capture the intricate flow features [20, 25, 26]. It may also be beneficial to enhance the ML model architecture by incorporating physics-informed constraints to help the model learn the flow dynamics from the governing equation [29, 30, 46, 55].

Another important limitation is that the ML models presented in this work are restricted to only a single patient-specific vessel geometry. Incorporating geometric parameters into the prediction process to improve model applicability to real-world scenarios is another avenue of future work. It should be noted that most published studies have focused on:

The prediction of steady-state flow fields in different blood vessel geometries [20, 22, 23, 33];The prediction of the time-dependent flow field in a fixed geometry under different flow conditions [25, 32, 45, 48].

To the best of the authors’ knowledge, Siena et al. [24] is the only study to develop an ML model to predict time-dependent flow fields in blood vessels with geometric variations, albeit considering the degree of stenosis as the only geometric parameter. Therefore, developing an ML model for time-dependent flow field prediction in blood vessels with practical geometric variation remains a novel and challenging task.

Conclusions

This study demonstrates the effectiveness of integrating POD-based ROM and neural network-based ML to predict WSS in blood vessels affected by vascular diseases. High-fidelity CFD simulations generated WSS data, which was then processed through POD to construct the ROM. The ML models were trained to predict the ROM coefficients from the inlet flowrate waveform which is a quantity that can be measured in the clinics. Two ML models were explored: the relatively simple flowrate-coefficients mapping model and the more advanced autoregressive model. Both ML models were then tested against two case studies: flow in PAD and flow in AD. The former served as a simpler case study, and the latter represented a more complex one. The key findings and implications of this study are summarized below:

The combined POD-ROM and ML framework successfully predicted WSS in patient-specific diseased vessels, demonstrating robust performance even with extremely limited training data.Between the two ML models tested, the simple flowrate-to-coefficient mapping model consistently outperformed the advanced autoregressive model across all error metrics in both case studies, suggesting that less complex models generalize better under data-scarce conditions.Prediction accuracy was high for the simpler PAD case and somewhat reduced for the complex AD case, indicating that flow complexity and data availability influence model performance.Despite the sparse training dataset, the ML-ROM approach reproduced the overall patterns of TAWSS and OSI in both vascular conditions, underscoring its reliability in capturing clinically relevant haemodynamic metrics.The proposed approach achieved speed-ups in the order of 10^4^ compared to full CFD simulations, dramatically reducing computation time and enabling near-real-time WSS predictions, which is highly relevant for clinical applications.

Future work could focus on expanding the training dataset using statistical generative techniques and incorporating geometric parameters to enhance model generalizability. Overall, this study highlights the promise of using ML models for fast, accurate predictions of haemodynamic quantities, potentially aiding in the diagnosis and treatment planning of cardiovascular diseases.

Supporting information

S1 FigReconstruction of time-dependent WSS for Case 1.(a) 3D reconstructions of WSS in the PAD under $[eqn]$ (test case) at four states of the cardiac cycle: Acceleration, peak systole, deceleration, and diastole shown in (b). (c) Mean absolute error from the prediction of WSS over a cardiac cycle. Gray area shows the range between 75th and 25th percentile of error.).(TIFF)

S2 FigReconstruction of time-dependent WSS for Case 2.(a) 3D reconstructions of WSS in the AD under $[eqn]$ (test case) at four states of the cardiac cycle: Acceleration, peak systole, deceleration, and diastole shown in (b). (c) Mean absolute error from the prediction of WSS over a cardiac cycle. Gray area shows the range between 75th and 25th percentile of error.(TIFF)