Data-Driven Pulsatile Blood Flow Physics with Dynamic Mode Decomposition

: Dynamic mode decomposition (DMD) is a purely data-driven and equation-free technique for reduced-order modeling of dynamical systems and ﬂuid ﬂow. DMD ﬁnds a best ﬁt linear reduced-order model that represents any given spatiotemporal data. In DMD, each mode evolves with a ﬁxed frequency and therefore DMD modes represent physically meaningful structures that are ranked based on their dynamics. The application of DMD to patient-speciﬁc cardiovascular ﬂow data is challenging. First, the input ﬂow rate is unsteady and pulsatile. Second, the ﬂow topology can change signiﬁcantly in different phases of the cardiac cycle. Finally, blood ﬂow in patient-speciﬁc diseased arteries is complex and often chaotic. The objective of this study was to overcome these challenges using our proposed multistage dynamic mode decomposition with control (mDMDc) method and use this technique to study patient-speciﬁc blood ﬂow physics. The inlet ﬂow rate was considered as the controller input to the systems. Blood ﬂow data were divided into different stages based on the inlet ﬂow waveform and DMD with control was applied to each stage. The system was augmented to consider both velocity and wall shear stress (WSS) vector data, and therefore study the interaction between the coherent structures in velocity and near-wall coherent structures in WSS. First, it was shown that DMD modes can exactly represent the analytical Womersley solution for incompressible pulsatile ﬂow in tubes. Next, our method was applied to image-based coronary artery stenosis and cerebral aneurysm models where complex blood ﬂow patterns are anticipated. The ﬂow patterns were studied using the mDMDc modes and the reconstruction errors were reported. Our augmented mDMDc framework could capture coherent structures in velocity and WSS with a fewer number of modes compared to the traditional DMD approach and demonstrated a close connection between the velocity and WSS modes. Author Contributions: Conceptualization, A.A. and S.T.M.D.; methodology, M.H. and S.T.M.D. and A.A.; formal analysis, M.H.; data curation, M.H.; writing—original draft preparation, M.H.; writing—review and editing, A.A. and S.T.M.D. All authors read agreed the


Introduction
Patient-specific computational fluid dynamics (CFD) has evolved over the past two decades to become a leading modality for obtaining 3D time-resolved blood flow data [1,2]. Advances in these simulations have recently led to the first ever food and drug administration (FDA) approval of computer simulations in evaluating the severity of coronary artery disease using CFD derived fractional flow reserve [3]. Thanks to today's computational power and high-performance computing, we are capable of generating hemodynamic data with extraordinarily high spatiotemporal resolution consisting of meshes with tens of millions of degrees of freedom and spanning tens of thousands of intra-cardiac time steps [4,5]. While these high-resolution simulations significantly improve the spatial in evaluating the performance of reduced-order models, and the authors highlighted the importance of preserving DMD modes with Ritz values close to the unit circle. These studies motivate the utility of DMD in cardiovascular flows, however, 3D patient-specific models of diseased arteries produce complicated flow structures and represent a challenge for the direct application of DMD. First, we deal with unsteady flow at the inlet due to the pulsatile nature of blood flow. Second, different parts of the cardiac cycle can have topologically very different flow patterns [41], which could challenge the linearization behind DMD. Finally, in cardiovascular flow problems, we are not only interested in coherent structures in velocity but also in the spatiotemporal dynamics of wall shear stress (WSS) vectors to study cardiovascular disease, which is typically localized at the vessel wall. To address the above issues, we propose a multistage DMD with control (mDMDc) approach and augment velocity data with WSS data. The unsteady flow rate produced by the heart is considered as a controller and DMD with control (DMDc) [27] is used to account for the incoming unsteady flow rate at the inlet. Subsequently, the cardiac cycle is divided into different stages based on the flow rate waveform and DMDc is applied to each stage separately. The goal of this study is to demonstrate that mDMDc could reveal the hidden low dimensionality in patient-specific blood flow data and therefore facilitate quantitative flow physics analysis.
In this study, we use 3D image-based models of a cerebral aneurysm and a coronary artery stenosis to evaluate our mDMDc approach. Interestingly, cavity flow [34] and jet flow [33] were among the early applications of DMD, which represent a simplified version of an aneurysm and stenosis, respectively. We first show that DMD modes match with the analytical Womersley solution [42] for pulsatile blood flow in tubes. Next, we demonstrate the performance of our mDMDc model in flow reconstruction, and subsequently, use our augmented mDMDc model to study coherent structures in velocity and WSS.

Computational Fluid Dynamics (CFD)
A cerebral aneurysm and a coronary artery stenosis model are considered in this study. Both CFD simulations are done using the Oasis solver, an open-source, minimally dissipative, finite-element CFD solver developed in FEniCS [43]. Oasis has been used in previous studies to obtain high-resolution blood flow data [5,44]. Herein, quadratic tetrahedral elements for velocity and first-order linear elements for pressure are used for high accuracy. Blood is assumed as Newtonian homogenous fluid and the walls are assumed rigid. Pulsatile parabolic profiles are used at the inlet. Both simulations are done for three cardiac cycles, and the last cardiac cycle is used in the data processing. The models and the flow waveforms at the inlet are shown in Figure 1.

Cerebral Aneurysm
An image-based cerebral aneurysm model from the online Aneurisk data set is selected. A 3.5 M tetrahedral element mesh with three layers of boundary meshing is created in SimVascular [45]. Thanks to FEniCS, quadratic shape functions are used, and therefore the mesh is approximately equivalent to using 28 M linear elements. A pulsatile flow rate waveform reported in a previous study [46] has been scaled in the range of the mean flow rate for our selected model according to the scaling law in Reference [47]. Zero traction is prescribed at the outlets. Each cardiac cycle is divided into 12,000 time steps.

Coronary Artery Stenosis
The coronary artery model used in a previous study, which includes the left anterior descending artery is considered [48]. The model is publicly available on the SimVascular website. A 50% artificial stenosis is created in the geometry. The mesh is created in SimVascular and consists of 3.31 M elements (equivalent to around 26 M linear elements once quadratic shape functions are applied in FEniCS).
A scaled version of the pulsatile flow rate reported in a previous study [49] is used for inlet boundary condition. The cardiac cycle is divided into 10,000 time steps and zero traction is prescribed at the outlets. Thanks to quadratic elements and high-performance computing, our spatial and temporal resolution for both models are an order of magnitude more resolved than most studies in the literature (see Reference [4] for an excellent discussion on this topic).

Dynamic Mode Decomposition with Control (DMDc)
DMDc was developed by Proctor et al. [27] as an extension of the original DMD method [50] to account for actuators and inputs to the system. This method can be described as a linear discrete state space system for a feedforward control problem as follows: where A ∈ R M×M and B ∈ R M×Q and A models the inherent dynamics of the system and B describes the effect of the input controller on the system. The state vector is defined as x ∈ R M and the input vector is defined as ψ ∈ R Q . The state vector x could be the variable of interest such as velocity or WSS and ψ represents the unsteady input flow. The snapshots of data in time (the selected time steps) are arranged in a matrix X and a shifted matrix X where X = x 2 x 3 · · · x N and X = x 1 x 2 · · · x N−1 . The x i columns could represent velocity vectors or velocity and WSS vectors combined as discussed below. A matrix Ψ is created based on the inputs as Ψ = ψ 1 ψ 2 · · · ψ N−1 . Herein, ψ i is a single entry column and represents the average velocity (or flow rate) at the inlet corresponding to the time-point of x i . The data-driven DMDc algorithm finds A and B matrices purely based on snapshots of collected data and input actuation. Based on the above definitions, Equation (1) could be rewritten as where C = [A, B] and Ω = X Ψ .
The best fit linear model C is obtained by minimizing the Frobenius norm of data reconstruction ||X − CΩ|| F . This problem could be solved using the truncated singular value decomposition (SVD) of the augmented matrix Ω = UΣV * ≈ŨΣṼ * , where Σ contains the singular values, the columns of U and V are the left and right singular vectors, respectively, * denotes the complex conjugate transpose, and ∼ specifies the truncated matrices from truncated SVD. This leads to the following calculation: Matrices A and B can be approximated by separating the complex conjugate transpose of the linear operatorŨ into two distinct components: whereŨ 1 ∈ R M×q ,Ũ 2 ∈ R Q×q . Finally, To find the reduced-order subspace of the snapshot matrix the following algorithm is used [27,31]: 1. Given data, construct snapshot matrices X, X , and Ψ. Then construct the augmented matrix Ω based on X and Ψ. 2. Compute the SVD of the augmented matrix Ω to get Ω ≈ŨΣṼ * with truncation value q.
Subsequently, based on the SVD results and the truncation value q, divide the complex conjugate transpose ofŨ into two distinct componentsŨ * = [Ũ * 1 ,Ũ * 2 ]. 3. Compute the SVD of X and obtain the truncated SVD decomposition X ≈ÛΣV * with truncation value r (Σ ∈ R r×r ). 4. Compute a reduced-order approximation of the linear operators A and B using: 5. Compute the spectral decomposition ofÃ asÃW = WΛ to find its eigenvectors (W) and eigenvalues (Λ). 6. Extract the dynamic modes of the operator A A more detailed discussion is provided elsewhere [10,27].

Dynamic Mode Decomposition in Cardiovascular Flows
As explained above, applying DMD directly to cardiovascular flows is challenging. An appropriate DMD model should account for the unsteadiness of blood flow due to the pulsatile nature of the heart. Additionally, different flow topologies during different phases of the cardiac cycle challenge finding optimal linear modes. Finally, we would like to study the coherent structures in velocity and WSS. We discuss how these issues could be handled.

Multistage dynamic mode decomposition with control (mDMDc)
To model the effect of pulsatile flow ejected by the heart, the input matrix could be constructed as Ψ = V mean 1 V mean 2 · · · V mean N−1 where V mean = Q/A inlet is the mean velocity at the inlet and is computed from Q, the volumetric flow rate at the inlet boundary. Therefore, the mean velocity at the inlet is assumed to be the input controller. To account for different flow topologies in different parts of the cardiac cycle, we use multiple time windows and apply DMDc to each window separately. The different stages of the cardiac cycle are shown in Figure 2 where each stage is colored differently. The coronary artery stenosis flow rate is divided into three parts (low flow rate, acceleration, and deceleration). Since the blood flow in the aneurysm model is more complex and chaotic, we used six stages to capture distinct flow topologies. •

Augmented mDMDc
To model coherent structures (from velocity) and near-wall coherent structures (from WSS) separately, we consider an augmented matrix of velocity and WSS to construct the state space matrix as: Subsequently, the augmented multistage dynamic mode decomposition with control (mDMDc) leads to the following equation:

Performance Evaluation
To evaluate the performance of mDMDc, the reconstruction error for velocity and WSS in the augmented system is calculated during different stages. After training the mDMDc model with velocity, WSS, and input data, the mDMDc algorithm can approximate A and B matrices in each stage based on Equation (5). Consequently, the velocity and WSS data can be reconstructed and predicted by Equation (1). A spatiotemporal error measure is defined (ST error ) by calculating the spatiotemporal average of the reconstruction error. Subsequently, the error is normalized by the spatiotemporal average of the data (ST data ) to obtain a normalized error (N error ) and compute a measure of accuracy (Acc). That is, where n is the number of predicted snapshots (time steps), m is the number of spatial points, and X CFD (i, j) and X mDMDc (i, j) are the ith snapshot (time) for the jth point predicted from CFD simulation and the mDMDc approach, respectively. The above analysis was done for the velocity and WSS data separately.

Results
In this Section, we first demonstrate the relationship between DMD and the classical Womersley analytical solution for incompressible pulsatile flow in straight tubes. Subsequently, the number of modes required for a specific level of accuracy is reported for velocity and WSS reconstruction. Finally, velocity and WSS streamlines for different modes are shown during different stages.

DMD and Womersley's Analytical Solution
Womersley derived an analytical solution for homogeneous, incompressible, Newtonian, and pulsatile flow in rigid and straight tubes [42]. The Womersley solution serves as a simple idealized model for pulsatile blood flow. The solution provides an analytical representation for spatial variation of velocity along the radial direction using Bessel functions and temporal variation using periodic Fourier modes. The solution is based on a given input pressure gradient waveform at the inlet, which needs to be represented based on Fourier series. Herein, we demonstrate that Womersley's analytical solution could be represented precisely using DMD modes.
In Figure 3, we sketch the Womersley solution and show its relationship with DMD modes. The dimensions for the tube and the input pressure gradient waveform are taken from a previous study [51]. The Fourier series frequencies and coefficients used to represent the pressure gradient waveform are shown. Additionally, the Womersley analytical solution and the resulting velocity profile at t = 0.3 s are shown. The final part of the figure reviews the DMD algorithm and shows the DMD results. The data matrix in DMD consists of velocity data at 101 temporal snapshots where each snapshot is sampled at 100 radial positions. The location of the resulting Ritz values (discrete eigenvalues) in the unit circle is shown. The velocity profile reconstructed based on each mode and eigenvalue is also shown. Finally, it is shown that by assembling all of the modes, the resulting DMD velocity profile exactly matches the analytical Womersley solution. The frequency parts of the dominant eigenvalues are listed, which match the frequencies in the pressure gradient waveform and Womersley solution with negligible error.

Data Reconstruction Accuracy
To evaluate the performance of our model, the minimum number of modes required to reach a certain level of accuracy is reported for each stage. The minimum number of modes for velocity and WSS in the coronary artery stenosis model is shown in Tables 1 and 2, respectively. Similar tables are also reported for velocity and WSS reconstruction in the cerebral aneurysm model (Tables 3 and 4). In general, velocity requires fewer number of modes compared to WSS for the same level of accuracy. In all cases, fewer than 10 modes are required for 90% accuracy. To compare the results with the traditional DMD method, Table 5 shows the minimum number of modes required for data reconstruction in each model for the DMD method without dividing the data into multiple stages, and two modes are shown in Figure 4. As expected, a large number of modes are required to reconstruct the data.

Flow Physics Using mDMDc
In the following figures, representative mDMDc velocity and WSS modes are shown in the coronary artery stenosis and aneurysm models. Each mode represents a different pattern within the data. Namely, each mode consists of the spatial pattern shown (eigenvector) and an associated temporal dynamics (eigenvalue). That is, based on its eigenvalue, each mode will oscillate and/or grow/decay in time with a frequency and rate specified by the eigenvalue. Therefore, using mDMDc we could characterize the spatiotemporal flow physics by dividing the flow into different modes.

Flow Physics Using mDMDc Modes in Coronary Artery Stenosis
The mDMDc modes for the coronary artery stenosis model during different stages are shown in Figures 5-7. Three modes are selected for each stage such that they exhibit mostly different flow patterns and can describe the dynamics of periodic flow (Ritz values lying close to the unit circle) [40]. The red dashed part of the waveform shows the stage of the cardiac cycle for each figure. The figures show the selected CFD data, reconstructed data with 95% accuracy, three mDMDc modes, and the location of the Ritz values (discrete eigenvalues) corresponding to the selected modes. The velocity and WSS streamlines are shown for front and back views. The first stage represents the systole phase of the cardiac cycle where low flow exists in the coronary arteries. Different flow and WSS patterns are seen in each mode. In each mode, regions of high WSS and velocity magnitude at the stenosis throat coincide. Additionally, the WSS vector directions match flow patterns near the wall. In the back view, a close correspondence between helical velocity streamlines and WSS streamlines could be seen in the second mode.
The second stage (flow acceleration) is shown in Figure 6. The close connection between flow and WSS topology could be seen in each mode. The comparison of each velocity and WSS modes could be used to observe how different patterns in velocity induce different WSS patterns. For example, in mode three, a stable fixed point in WSS could be seen, which corresponds to upwelling motion in the velocity streamlines [52]. Finally, similar results are presented for the third stage (flow deceleration) in Figure 7. Different patterns could be seen in different modes and a clear connection between near-wall helical flow and helical WSS patterns with spiral WSS fixed points could be seen in mode one.

Flow Physics Using mDMDc Modes in Cerebral Aneurysm
The mDMDc modes for the cerebral aneurysm model during three representative stages are shown in

Discussion
DMD is a purely data-driven and equation-free technique to reduce the snapshots of data into dynamically coherent modes. We have shown how DMD could be properly modified and applied to pulsatile cardiovascular flows. The incoming transient pulsatile flow waveform in arteries makes the direct application of DMD challenging. DMD seeks to find an optimal linear model based on given data. It may be perceived that linearizing nonlinear flow where the flow behavior jumps between different nonlinear regimes is not an easy task. Dividing the cardiac cycle into different stages based on the flow waveform mitigates this inherent difficulty in applying DMD to these types of flows.
Additionally, using DMDc allows one to separate the effect of unsteady inflow and the inherent dynamics of the system.
The application of DMD to cardiovascular flow physics raises a major question: why is DMD an appropriate framework to represent velocity and WSS data from patient-specific cardiovascular flows?
To answer this question we demonstrated that the classical Womersley solution could be represented using DMD modes and DMD analysis precisely identifies the frequencies in the Womersley solution.
In the Womersley solution, the temporal variation of velocity is represented using Fourier series where the frequencies are based on the frequencies of the input pressure or flow waveform. DMD modes are known to approximate Fourier modes [21,53], therefore we may think of DMD as a data-driven counterpart of the Womersley solution in patient-specific geometries. In other words, in each mode, the eigenvectors represent the spatial approximation to the solution and the eigenvalues correspond to the frequencies. In DMD, each temporal frequency in the solution has its unique spatial velocity pattern, whereas, in the Womersley solution, all frequencies are associated with Bessel functions representing the spatial velocity profile. Having an analytical solution is desirable because it guides a quantitative interpretation of the flow physics. Obviously, it is not possible to find an analytical solution for patient-specific geometries, however, given numerical or experimental solution, we may expect to approximate it with DMD modes to provide a quantitative model for the evolution of coherent structures.
Our results showed that applying mDMDc significantly improves flow reconstruction compared to the traditional DMD approach. Namely, as shown in Table 5, a considerably large number of modes are required to reconstruct the velocity data in both models when DMD is applied directly to the entire data. For example, the cerebral aneurysm model required 155 modes for 98% velocity reconstruction accuracy. Interestingly, a prior DMD study in 2D left ventricles required less than 39 modes for 98% accuracy in capturing the kinetic energy [40], which demonstrates the complexity in our 3D patient-specific models. Nevertheless, using the mDMDc approach reduces the number of required modes for our aneurysm model to 3-13 modes depending on the stage. Additionally, the use of an augmented system enables more accurate mode reconstruction for WSS. Specifically, applying mDMDc directly to WSS data in the coronary artery stenosis model resulted in highly unstable modes with high-frequency structures that compromised flow reconstruction. For the cerebral aneurysm model, we observed that the augmented system required on average 6 and 8 fewer number of modes to obtain 99% and 99.5% accuracy in WSS, respectively (results not shown). We also observed that uniformly scaling the input controller (u) could improve reconstruction error in some cases. Such scaling is equivalent to changing the weight in the optimization problem for finding matrices A and B in DMDc.
The relationship between DMD and mDMDc modes in complex data is not trivial. We observed that the mDMDc modes/eigenvalues, in general, are different from any subset of DMD modes/eigenvalues. It is not guaranteed to obtain the same results as a multistage variant of DMD just from knowing which DMD modes (identified from the full dataset) are active in a certain stage. In the ideal cases such as the Womersley solution (sinusoidal data), we saw that mDMDc and DMD modes are aligned together. However, this could not be seen in our patient-specific models.
The mDMDc results shown in the figures show how the spatiotemporal evolution of disturbed blood flow in diseased arteries is made up of different coherent structures (modes) with different flow features. For instance, different types of fixed points in the WSS vector field could be seen in different modes. These WSS fixed points are closely related to flow topology away from the wall [52] as could be seen in the results. For instance, source and sink type fixed points in the WSS vector field correspond to flow towards the wall (impingement) and flow away from the wall (upwelling motion). In general, the WSS modes without fixed points correspond to unidirectional attached flow in the velocity modes. On the other hand, the appearance of WSS fixed points is related to more complex flow patterns as seen in the velocity streamlines.
The application of DMD to diseased arterial flows could provide different types of information. From a mathematical viewpoint, our mDMDc results show that despite the apparent spatiotemporal complexity of blood flow in aneurysms and stenoses, the spatiotemporal evolution of hemodynamic data could be represented with a few modes where each mode evolves in time according to the linear model discovered with DMD. From a flow physics viewpoint, dividing the flow into distinct modes where each mode evolves in time with a fixed oscillation frequency and growth/decay rate significantly simplifies the flow physics. One could identify the modes that decay in time and therefore identify flow structures that dynamically do not persist in the flow. Additionally, the combination of different modes with eigenvalues close to the unit circle represents the flow patterns that persist and oscillate in time. These modes were selected and shown in the results. It can be observed that most modes posses a unique flow pattern. What distinguishes DMD from other tools in studying flow physics is that in DMD each mode is equipped with an eigenvalue that defines the quantitative temporal evolution of the mode. In other words, DMD modes and the flow patterns extracted from each mode are quantitative. For example, if a vortex is identified from one mode, its evolution in time could be quantitatively explained. This is something that other flow physics identification tools, for example, Lagrangian coherent structures (LCS) [54,55], do not offer. One could, of course, compute the LCS from the velocity [41] and WSS [56] modes to rigorously identify flow structures and use the eigenvalues associated with each DMD mode to equip the qualitative flow structures identified from LCS with a quantitative model.
Another interesting potential application of DMD is future forecasting. Once the DMD model is established, one could use the model to predict blood flow beyond the available data. For example, the waveform could be extended (e.g., increase in the peak flow rate) and the DMD model could be used to extend the data based on the new waveform. Dividing the cardiac cycle into multiple stages enables one to extend the waveform in different parts of the flow rate, and therefore study variations in the flow waveform. Of course, it is highly likely that only a short time window could be predicted before the DMD model deviates from the physical system. Nevertheless, such a model could be used in uncertainty quantification where small perturbations to the flow rate need to be examined. Additionally, using DMD on the fly in a CFD code could be used to design a new time integration strategy where the solver could switch between the high-fidelity CFD solver and the DMD model during different time intervals.
There are different possibilities to improve and extend our mDMDc model in future work. Our preliminary investigation revealed that mDMDc accuracy is reduced when the flow rate is sufficiently increased. Transitional and unstable flow patterns have been observed in prior high-resolution CFD studies [57]. These unstable flow patterns challenge the task of finding an optimal linear model in DMD. One solution is to develop hybrid reduced-order models leveraging the optimal reconstruction property in POD [21,58,59]. Additionally, using robust principal component analysis (RPCA) in conjunction with DMD could potentially remove high-frequency components in velocity if they are spatially sparsely distributed [60]. Other developments in DMD such as multiresolution DMD [28] and extended DMD [61] could also be incorporated within the mDMDc framework to improve the method's accuracy. Finally, it is possible to use other parameters as the input controller in the DMDc framework. For instance, viscosity or vessel stiffness (in a fluid-structure interaction problem) could be varied to build a reduced-order model with these parameters as the input in DMDc.

Conclusions
Dynamic mode decomposition (DMD) was used to reconstruct hemodynamic data in image-based coronary artery stenosis and cerebral aneurysm models. A close relationship between DMD and Womersley solution was shown. A new multistage DMD with control (mDMDc) technique was proposed to overcome challenges in the direct application of DMD to pulsatile blood flow data. The results demonstrated that in most stages of the cardiac cycle velocity and WSS data could be reconstructed using less than 10 modes with an error of less than 10%. The flow patterns were studied using DMD modes to extract different spatiotemporal flow structures in velocity and WSS. A close correspondence between coherent structures identified in velocity modes and near-wall coherent structures identified in WSS modes was observed. In summary, this work demonstrated that DMD could be customized to study flow physics in 3D patient-specific cardiovascular flows.