Modeling of a Simplified 2D Cardiac Valve by Means of System Identification

Abstract

The function of a heart valve is to maintain the correct circulation of blood from the heart to the entire Cardiovascular System (CVS). The CVS has become a main target in the health science and engineering, due to the importance of the detection anomalies or pathologies. Some techniques that deal with these problems, like Magnetic Resonance Imaging (MRI), Ultrasound and Electrocardiography (ECG) among them, require important investments to acquire measurement instruments. Numerical simulations have arisen as an alternative, by means of Fluid-Structure Interaction (FSI) analysis. The combination of computational methods, such as Finite Element Method (FEM) and Incompressible Computational Fluid Dynamics (ICFD), have allowed the development of FSI software, which leads to design versatility and economic advantages. However, computational time is high even for simplified designs. Lumped models appear as another solution. They consist in simplified zero-dimensional models that allow to mimic the Cardiovascular System (CVS) with low computational cost and good approximations. Nevertheless, they require a parameter identification and the calculation process for each scenario, which becomes increasingly complicated when the number of parameters is high. Therefore, an alternative approach is proposed in the present manuscript. Taking the volumetric flow rate and the Pressure Difference (PD) as the input and output parameters, respectively, of the flow across a heart valve of simple geometry, a System Identification (SI) is performed. The proposed methodology requires low computational cost and provides a good approximation. Results obtained guarantee a 99.76% of goodness of fit.

1 Introduction

Within the CVS, the purpose of a heart valve is to maintain the correct circulation of blood, from the heart to the entire body. Pathologies presented in one or more of the four heart valves can lead to a critical malfunctions of the heart, such as stenotic or back-flow complications. Many of these pathologies arise as consequences of an increment of the magnitude of the shear stresses acting on the valves, due to the blood flow.

For these reasons, many attempts to create artificial valves have been tried, with special interest on Biological Heart Valves (BHV) and Mechanical Heart Valves (MHV). MHV are usually fabricated with inorganic and stiff materials, whereas BHV are obtained directly from animals or fabricated with organic and flexible tissues. According to [25], BHV has outstanding mechanic and hemodynamic properties due to the bio-materials used for their construction. Contrary to BHV, MHV designs minimize back-flow and last longer lifespans, but they produce deficient physiological flow conditions such as high stress levels that lead to the destruction of red blood cells [14, 21].

Considering that the design of MHVs is based on mathematical modeling, they can be tested by using numerical simulations. The methodology employed to simulate the performance of valves is known as computational FSI, and consists of a combination of ICFD to calculate the fluid motion and a force balance to compute the leaflets rotation. This alternative allows us to modify and create new geometries by using a Computer-Aided Design (CAD) software. Studies on the flow across MHVs have been carried out in [2, 9, 10, 12]. Nevertheless, a numerical simulation of flexible leaflet materials is completely different, since their motion depends strongly on the interaction with the fluid. In this case, the dynamics of the mechanical valve can be simulated using time-dependent FEM combined with ICFD. Many works using this strategy are found in literature using different techniques [1, 8], for instance an Arbitraty Lagrange-Eulerian (ALE) [6] method, the Fictitious Domain Method (FDM) and the Boundary Method (IBM) in [4, 5, 21] and [19, 22], respectively.

Another approach to represent the CVS is by using the zero-dimensional lumped models. They are based on simple ordinary differential equations and they consider an average value of the CVS local variables (pressure and velocity) within each of their components (organs and blood vessels). Different variations, such as the Windkessel Model (WM) allows to mimic the systemic arterial tree [3, 7, 17, 18, 20, 23, 24]. However, their main disadvantage in these approaches is the need to identify and calculate all the parameters to adjust to a particular physiological condition.

Considering all previous methods, the main contribution of this work relays on the implementation of an alternative method using SI to estimate the behavior of a cardiac valve under certain blood flow rate inside of a channel. This method allows to replicate the output, pressure difference PD measured across the valve, with a lower computational cost than the one required by FSI, and with no requirement of parameters identification and calculations, as required by the lumped models such as WM.

In Sect. 2, a new approach to study them using AutoRegressive (AR) system identification is proposed. A comparison between the results obtained with FSI and AR is shown in Sect. 3. Finally, future work and conclusions of this paper are listed in Sect. 4.

2 Theoretical Background and Methodology

2.1 Computational FSI

As mentioned in the previous section, computational FSI is one of the main alternatives to understand how the heart valves behave under certain fluid conditions. As a first step, a simplified 2D cardiac valve has been studied by means of a numerical simulation, which corresponds to the geometry and flow conditions presented in a previous experimental work [11]. Parameters used for the experiment can be described as (a) Fluid domain $(\varOmega _{f}) \in$ $\lbrace U_{f}, V, f, \rho _{f}, \mu \rbrace$, where the elements are inlet velocity, stroke volume, frequency, fluid density and dynamic viscosity, respectively, and (b) Solid domain $(\varOmega _{s}) \in$ $\lbrace d, l, E, \rho _{s}, PR \rbrace$, where the elements are leaflet thickness, length, elastic modulus, material density and Poisson ratio, respectively.

The simulation has been performed with the following boundary conditions for the fluid domain (see Fig. 1): Inlet Wall representing the entrance for the fluid, Outlet Wall where the fluid goes away and the pressure reaches zero, Non-Slip Condition Wall representing that the fluid has zero velocity relative to the boundary. The valve is also configured as a Non-Slip Condition Object.

Fig. 1.

Boundary conditions used for a channel with a simplified cardiac valve design inside. The channel is in terms of semi-height h, where $h=15$ mm, width $w=50$ mm.

Inlet velocity curve is formed by the two stages during the cardiac cycle. Systole $t_{s}$ is the positive velocity pushing forward the volume V during 0.35Tc and diastolic part $t_{d}$ is a momentary negative pulse that moves back the fluid and then goes to zero velocity at the end of cardiac cycle. Magnitude of the negative pulse is controlled by a scale factor $\frac{\alpha }{\gamma }$, where $\alpha = \frac{ts}{td}$ and $\gamma$ is a scalar positive value. The velocity used for the experiment is obtained by initially considering the following volume formula:

$$\begin{aligned} V= A \int _{0}^{Tc} U \cdot g(t) dt \end{aligned}$$

(1)

where $g(t)=1$ during $t_{d}$ (represents a unitary pulse) and $g(t)=-\alpha$ during $t_{s}$ (see Fig. 2a), A is the cross-sectional area $A=hw$. By solving the integral divided into the two stages, we obtain the following:

$$\begin{aligned} U_{s}&= \frac{V}{A t_{s}}&U_{d}&= \frac{V}{A \alpha t_{d}} \end{aligned}$$

(2)

2.2 System Identification

The hereditary identification principle uses the same idea by taking the experimental expectation instead of the mathematical one [16]. In this case, the distance among estimated trajectories is minimized. Hence,

$$\begin{aligned} E^t[x_\bullet ]=\frac{1}{t}\sum _{\tau =1}^tx_\tau \end{aligned}$$

(3)

is the trajectory $x_{\tau }$ experimental expectation over [1, t] and $x_\bullet =\{x_\tau ,\tau =1,...,t\}$. Considering the optimal estimator trajectory $\hat{y}_\tau ^t$ at instant t for data over [1, t] and the past estimator values:

$$\begin{aligned} \hat{y}_{\tau |\tau -1}^t=\sum _{i=1}^{n}a_{i,t}\hat{y}_{\tau -i|\tau -i-1}^{t-i}+\sum _{i=1}^{n}k_{i,t}(y_{t-i}-\hat{y}_{t-i|t-i-1}^{t-i}) \end{aligned}$$

(4)

If the regression vector trajectory is defined by:

$$\begin{aligned} \varphi _\tau ^{t-1}=[\hat{y}_{\tau -1|\tau -2}^{t-1},\cdots ,\hat{y}_{\tau -n|\tau -n-1}^{t-n},\tilde{y}_{\tau -1|\tau -2}^{t-1},\cdots ,\tilde{y}_{\tau -n|\tau -n-1}^{t-n}] \end{aligned}$$

(5)

minimizing the estimation experimental error yields:

$$\begin{aligned} \theta _t=E^t[\varphi _\bullet ^{t-1}(\varphi _\bullet ^{t-1})^T]^{-1}E^t[\varphi _\bullet ^{t-1}y_\bullet ] \end{aligned}$$

(6)

The parameters computation is hereditary; therefore, it is required to know the experimental autocorrelations of $\varphi _\tau ^{t-1}$, plus $\varphi _\tau ^{t-1}$ and $y_\tau$ inter-correlations. Hence, to compute the system parameters at instant t, it is necessary to keep in memory all the past trajectories $\hat{y}_{\tau |\tau -1}^{t-1},\cdots ,\hat{y}_{\tau |\tau -1}^{t-n}$.

2.3 ARMAX Model

The hereditary algorithm can be extended when a system exogenous measured input is available [15]. It is then necessary to consider the following estimator model:

$$\begin{aligned} \begin{array}{c} \hat{y}_{\tau |\tau -1}^t=\displaystyle \sum _{i=1}^{n}a_{i,t}\hat{y}_{\tau -i|\tau -i-1}^{t-i}+\sum _{i=1}^{n}k_{i,t}(y_{t-i}-\hat{y}_{t-i|t-i-1}^{t-i})\\ +\displaystyle \sum _{i=0}^nb_{i,t}u_{t-i} \end{array} \end{aligned}$$

(7)

If we define two vectors $x_\tau$, $y_\tau$ external product such as:

$$\begin{aligned} \langle x_\bullet |y_\bullet \rangle ^t=\frac{1}{t}\sum _{\tau =1}^tx_\tau ^Ty_\tau \end{aligned}$$

(8)

The optimal estimator can be seen as the output trajectory projection over the following trajectories:

$$\begin{aligned} \hat{y}_\tau ^t=P_\tau ^t(y_\bullet |\hat{y}_{\bullet -1|\bullet -2}^{t-1},\tilde{y}_{\bullet -1|\bullet -2}^{t-1},u_\bullet ,\cdots ,\hat{y}_{\bullet -n|\bullet -n-1}^{t-n},u_{\bullet -n+1}) \end{aligned}$$

(9)

It is important to note that the criterion minimization task is quadratic; as a result, it is not necessary to derivate it. The nonlinear parameters dependence is solved by the hereditary algorithm linear augmenting memory.

2.4 Hereditary Identification Algorithm

• 1. Initialization. As data before $t=1$ is unknown, the projection over it is zero (division by a covariance à priori infinite). It is then possible to think of it as zero:

  $$\begin{aligned} \hat{y}_{\tau |\tau -1}^{1-i}, \quad \forall \tau <1, \forall i=1,\cdots ,n \end{aligned}$$

  (10)

• 2. At instant $t-1$ we have data from the already optimized trajectories $\hat{y}_{\tau |\tau -1}^{t-i}, \forall \tau =1,\cdots ,t-1, \forall i=1,\cdots ,n$

• 3. At instant t, new data $y_t, u_t$ are available in order to update the experimental correlation matrix and the inter-correlation vector, both defined in Eq. (6), by hereditary computation:

  $$\begin{aligned} E^t[\varphi _\bullet ^{t-1}(\varphi _\bullet ^{t-1})^T]=\frac{1}{t}\displaystyle \sum _{\tau =1}^t\varphi _\bullet ^{t-1}(\varphi _\bullet ^{t-1})^T \end{aligned}$$
  $$\begin{aligned} E^t[y_\bullet \varphi _\bullet ^{t-1}]=\frac{1}{t}\displaystyle \sum _{\tau =1}^ty_\bullet \varphi _\bullet ^{t-1} \end{aligned}$$

  (11)

• 4. We obtain the new system parameters at instant t by inverting:

  $$\begin{aligned} \theta _t\hat{=}E^t[\varphi _\bullet ^{t-1}(\varphi _\bullet ^{t-1})^T]^{-1} \times E^t[y_\bullet \varphi _\bullet ^{t-1}] \end{aligned}$$

  (12)

• 5. Hereditary part. With the obtained parameters, it is then possible to compute the new estimator trajectories $\hat{y}_{\tau |\tau -1}^{t-1}$ based on the trajectory $\hat{y}_\tau ^{t-i}$ computed at the preceding stage by employing the model $\forall \tau =1,\cdots ,t$:lh

  $$\begin{aligned} \hat{y}_{\tau |\tau -1}^t=\displaystyle \sum _{i=1}^{n}a_{i,t}\hat{y}_{\tau -i|\tau -i-1}^{t-i}+ \displaystyle \sum _{i=1}^{n}k_{i,t}(y_{t-i}-\hat{y}_{t-i|t-i-1}^{t-i})+\sum _{i=0}^nb_{i,t}u_{t-i} \end{aligned}$$

  (13)

• 6. Return to step 2.

3 Results

Parameters for this simulation are the following. $(\varOmega _{f})$ = $\lbrace$0.09641 m/s, 3.8 $x10^-5$ m$^3$, 0.666 Hz, 1000 kg/m$^3$, 0.001 Pa$\cdot$ s$\rbrace$ and $(\varOmega _{s})$ = $\lbrace$0.0004 m, 0.26 m, 2.15 $x 10^6$ Pa, 1070 kg/m$^3$, 0.47$\rbrace$. The cardiac cycle for this test is Tc=1.5015 seconds and by substituting Eq. (2) we can find that $\vert U_{s} \vert \approx \vert U_{d} \vert = 0.09641$ m/s. The curve is designed using hyperbolic tangents trying to fit the respective systole and diastole duration (see Fig. 2a). The magnitude of this curve is unitary and multiplied by $U_{s}$ and $\frac{U_{d} \alpha }{\gamma }$, depending on the stage.

FSI results are compared with experiments in [11] by measuring two parameters. Opening Area (OA), indicates how much the valves opened during systole stage (see Fig. 2b). Pressure difference is also measured into two locations of the channel, $x=[-5h,5h]$ and $y=0$ (before and after the valves).

Fig. 2.

(a) The curve describes two stages of the cardiac cycle. ts and td represents systole and diastole time respectively. $\alpha {/}\gamma$ represents the scaling factor for the negative pulse. (b) Opening and closure steps of a heart valve during cardiac cycle.

Considering that our simplified valve only moves along the $x-y$ plane, the fluid only interacts into two directions and not along the width as in the experiment. This could lead to some variations as the ones seen in magnitude for the OA observed in Fig. 3a or the falling part during diastole in the pressure difference (blue curve) of Fig. 3b. Another consideration to consider is that the experiment results are taken from a Particle Image Velocimetry (PIV) system during certain time-steps and the number of samples of the output signal is very limited compared to the ones obtained in simulations. Despite these considerations, results seen are qualitative correlated to the experimental 3D set-up.

Fig. 3.

(a) As it can be seen in Figure 2b the maximum opening area occurs between $0.15Tc - 0.3 Tc$ which approximately represents half of the systole. (b) PD observed in both curves during the first stage is almost identical. The falling part seen during diastole in the simulation curve is cleaner due to lack of external uncontrollable disturbances that were present in experimental test. Experimental data is retrieved from [11]. (Color figure online)

In order to implement the SI method, it is necessary to conduct more simulations to have a test/validation process. Hence, a set of simulations are conducted by changing the frequency of the fluid $f=[0.1,0.2,0.3,...,0.666,...,0.8,0.9,1] Hz$ (these f change the velocity $U_{f}$ as well) and leaving all the remaining parameters as they were.

3.1 Identification Results

The figure of merit considered in this work, called the multiple correlation coefficient, assesses the algorithm quality by indicating the percentage of $y_\tau$ that has been well explained by $\hat{y}_\tau ^t$ [13]:

$$\begin{aligned} R=\left( 1-\frac{\displaystyle \sum _{\tau =1}^t(y_\tau -\hat{y}_\tau ^t)^2}{\displaystyle \sum _{\tau =1}^t y_\tau ^2}\right) \times 100 \end{aligned}$$

(14)

The 11 set of samples obtained from the set of simulations were concatenated in order to use 45% of them for training the model and 55% to validate it afterwards.

3.2 Validation

The identified ARMAX transfer function (see Eq. (15)) used for validation deals with two components, one for the input plus another for tracking the error or innovations within the process and its respective polynomials are given in Table 1.

$$\begin{aligned} ARMAX(s)&=\frac{B(s)}{1-[A(s)-K(s)]} + \frac{K(s)}{1-[A(s)-K(s)]} \end{aligned}$$

(15)

Table 1. Polynomials for validation

The identification process employed a dimension n $=$ 2, given the CVS behavior did not require further exploration (Fig. 4).

Fig. 4.

(a) Difference of pressure observed for the cardiac valve during the first two sets of simulation. (b) Estimation of the output signal for the second set of simulations.

4 Conclusions

SI has proved to be a suitable alternative to estimate how a valve is going to behave under certain conditions. By using an inlet flow rate curve as the only input, the method can be trained to obtain a good approximation on the pressure difference, the only output, without any complications. It is important to remark that this method only works if any of the parameters inside $(\varOmega _{f})$ changes; hence, by changing only the frequency f ($\vert U_{f} \vert$ is directly affected), the proposed method can obtain results with a 99.76% of goodness of fit with the FSI simulations, between a set of frequencies established for the same elasticity modulus. In order to solve this limitation, future work includes the implementation of a method that considers changes not only in the fluid but in the solid as well. This will guarantee to achieve estimations; for example, for different elasticity modulus or densities and even changes in physical properties such as density or length.