Design of Epicardial Restraints for Optimized Passive Filling of the Right Ventricle

1. Introduction

End-stage heart failure (HF) can be effectively treated with cardiac transplantation, but this approach has limitations and is not available for every patient in need. Alternatively, ventricular restraint therapy provides a non-transplant option that uses prosthetic materials to mechanically stabilize the heart, curb further deterioration, and enhance patient outcomes [1]. The benefits of passive epicardial constraints on the failing left ventricle (LV) have also been exhibited using finite element analysis (FEA) [2]. However, each ventricle responds differently to epicardial constraints and the LV can withstand higher levels of such therapy [3]. It was already demonstrated that the application of adjustable ventricular restraints exclusively to the dilated and failing right ventricle (RV) in a sheep model of right ventricular (RV) failure resulted in an acute improvement of ventricular function with the potential to induce RV reverse remodeling, similar to what was observed in the LV during left ventricular (LV) failure [4].

RV computational modeling has the potential to demonstrate the complex mechanisms of RV adaptation in response to pressure overload triggered by pulmonary arterial hypertension (PAH) [5]. A structure-based modeling framework was used to analyze fiber, tissue, and organ-level adaptations of the RV in PAH, revealing that pressure overload-induced changes in geometry and wall stiffness significantly affect contractility [6,7]. A data assimilation method was used to quantify changes in regional myocardial properties and stresses in PAH, highlighting the relationship between right ventricular free wall (RVFW) contractility and the degree of remodeling [8]. Biventricular finite element models of rats with PAH showed that RV wall thickening, increased heart muscle contractility, and greater stiffness work together to mitigate the increased stress and strain of the wall caused by the pressure excess [9]. The application of more detailed FEA-based computational models to assess HF in individuals with pulmonary hypertension (PH) revealed significant improvements in RV function and reduced risk of failure after mitral valve surgery [10].

Computer-aided simulations revealed that some HF therapies for the LV can have adverse effects on the RV. A computational framework, incorporating both lumped-parameter and FEA, was used to evaluate the performance of the RV after the addition of the continuous-flow left ventricular assist device (CF-LVAD) in patients with end-stage heart failure [11]. Another biventricular FEA model of an end-stage failing heart assisted with a left ventricular assist device (LVAD) revealed that operating the LVAD at high speeds can lead to RV dysfunction due to LV unloading [12]. At speeds higher than 4k rpm, RV contractility decreases, especially when the septal wall is thin [13]. The performance of cardiac devices was analyzed using FEA to evaluate the balance between the hemodynamic and biomechanical impact of RV support in the context of LVAD-associated failure [14].

FEA has already been used to simulate repaired cardiac ventricles. Using an FEA model of the LV, it was demonstrated that increasing the size and stiffness of pericardial patches used in surgical repair of ventricular septal rupture following myocardial infarction (MI) significantly reduces stress in the suture line, thereby minimizing the risk of residual leakage [15]. Computational models of the infarcted LV in canines highlight the potential of mechanical modifications to improve heart function by reducing fiber stress and reshaping ventricular geometry, emphasizing the importance of tailoring therapeutic strategies [16,17]. FEA was also used to evaluate the restoration of LV function after implanting a contractile fibrous cardiac patch for chronic MI, focusing on different fiber orientations [18]. Numerical simulations could help investigate the effects of epicardial force patches on ventricular mechanics [19]. The design and evaluation of the RV external device (RVEX), a biomimetic elastic sleeve, has been documented to reinforce and control RV function in RV failure [20]. This was achieved by demonstrating hemodynamic improvements in both the left and right ventricles. In addition to epicardial restraints and patches [21], some computational studies focused on modeling injectable hydrogels [22] for infarcted ventricles [23,24]. However, since the RV has a relatively thin wall, injection is challenging, and patching is preferable for this ventricle.

Despite the development of various types of cardiac restraints and reinforcements to support failing hearts mechanically [25], there is still a need for assessing the clinical effectiveness of ventricular restraint devices and other mechanical supports in preventing adverse LV remodeling post-MI [26,27]. Further advancement of these therapies requires in-depth research into the changes in the heart following MI through collaboration across multiple disciplines [28]. Although previous computational studies have provided invaluable insights into the effect of patch fiber orientations on the acutely infarcted LV [16,18] or the impact of myofiber structure on the passive filling of the unloaded LV [29,30,31], such evaluations are unexplored for the RV. Moreover, different epicardial restraints have different mechanical properties and the optimal conditions for each should be determined.

Here, we propose a framework for designing optimal passive epicardial restraints for RV failure. Our method will be used to find the most effective restraints applied to unloaded biventricular FEA models using a variety of available materials. To examine the effectiveness of the therapy, the fiber stress and strain of the constrained ventricles will be analyzed. To the best of our knowledge, this is the first study to utilize computational modeling to identify the biomechanical properties and the optimal epicardial restraint for the RV, taking into account the unloaded state of the cardiac ventricles.

2. Materials and Methods

2.1. FEA Models from Imaging

We selected three patients diagnosed with RV failure for our study. Biventricular heart geometries were created based on the segmentation of cardiac magnetic resonance (CMR) images using the freely available software Segment Research (Medviso AB, Lund, Sweden) [32]. After smoothing and simplifying the 3D extracted ventricles, they were meshed with linear tetrahedral elements employing hybrid formulation (C3D4H). The total number of elements and nodes of all FEA models are shown in Table 1. Subsequently, material orientations were defined for each element. It was assumed that fiber angles changed smoothly from -60 to 60 degrees when moving from the endocardium to the epicardium, while sheets were perpendicular to each of these surfaces. A fixed reference point (RP) was created at the basal center of each ventricle. The RP of the LV was transitionally coupled with the basal nodes [33]. Ventricular blood pressure was exerted through the fluid cavity interaction [34] to the RPs. Valvular stiffness was added to the biventricular models using springs [35,36]. The LVEDV and RVEDV were calculated from the segmentation-based models at end-diastole (ED). Ventricular pressures were estimated noninvasively based on the reports provided by echocardiographic measures. For the LV, e/e’ helped us find the LVEDP [37,38,39], while the right ventricular end-diastolic pressure (RVEDP) was approximated through the right atrial pressure (RAP) according to the previous estimation guidelines [40]. RAP was determined from the pressure gradient across the tricuspid valve (ΔP) and the right ventricular systolic pressure (RVSP) [41,42,43,44,45]. End-diastolic (ED) pressures and volumes are shown in Table 1. To examine the effect of RVEDP on the optimization results, a second model of the third patient was developed with an RVEDP of 10 mmHg (referred to as patient 3-H), despite the initial estimate of 5 mmHg for this patient (patient 3). Due to the relatively low RVEDP of patient 3, we decided to examine the performance of his heart under both actual and elevated pressures.

2.2. Passive Mechanical Response

The myocardium was modeled as an anisotropic hyperelastic material using the invariant-based strain energy density function as defined in Equation (1) [46].

$$\psi ={\displaystyle \frac{{\widehat{a}}_{\mathrm{m}}}{2{\widehat{b}}_{\mathrm{m}}}}\left\{\mathrm{e}\mathrm{x}\mathrm{p}\left[{\widehat{b}}_{\mathrm{m}}\left(I_1-3\right)\right]-1\right\}+\sum _{i=\mathrm{f},\mathrm{s}}{\displaystyle \frac{{\widehat{a}}_i}{2{\widehat{b}}_i}}\left\{\mathrm{e}\mathrm{x}\mathrm{p}\left[{\widehat{b}}_i{\left(\mathrm{m}\mathrm{a}\mathrm{x}(I_{4i},1)-1\right)}^2\right]-1\right\}+{\displaystyle \frac{{\widehat{a}}_{\mathrm{fs}}}{2{\widehat{b}}_{\mathrm{fs}}}}\left\{\mathrm{e}\mathrm{x}\mathrm{p}\left[{\widehat{b}}_{\mathrm{f}\mathrm{s}}{I_{8\mathrm{f}\mathrm{s}}}^2\right]-1\right\}$$

(1)

Where the ${\widehat{a}}_{\mathrm{m}}$, ${\widehat{b}}_{\mathrm{m}}$, ${\widehat{a}}_f$, ${\widehat{b}}_f$, ${\widehat{a}}_s$, ${\widehat{b}}_s$ , ${\widehat{a}}_{\mathrm{f}\mathrm{s}}$, and ${\widehat{b}}_{\mathrm{f}\mathrm{s}}$ are eight positive material constants. The $\widehat{a}$ parameters have stress dimensions, whereas the $\widehat{b}$ parameters are dimensionless. An orthonormal basis consisting of fibers ($f_0$) and sheets ($s_0$) directions was used in the formation of the first ($I_1$), fourth ($I_{4f},I_{4s}$), and eighth ($I_{8\mathrm{f}\mathrm{s}}$) invariants of the right Cauchy-Green tensor. For $j$ representing $\mathrm{m}$, $\mathrm{f}$, $\mathrm{s}$, and $\mathrm{f}\mathrm{s}$, scale factors $A$ and $B$ were defined so that ${\widehat{a}}_j=A{a}_j$ and ${\widehat{b}}_j=B{b}_j$. Basic material parameters $a_j$ and $b_j$ came from curve fitting to the experimental data from six shear modes of cubic specimens obtained from the myocardial wall, being $a_m=1.05\mathrm{k}\mathrm{P}\mathrm{a},b_m=7.542,a_f=3.465\mathrm{k}\mathrm{P}\mathrm{a},b_f=14.472,a_s=0.481\mathrm{k}\mathrm{P}\mathrm{a},b_s=12.548,a_{fs}=0.283\mathrm{k}\mathrm{P}\mathrm{a},b_{fs}=3.088$ [47,48]. The restraint was assumed to be isotropic and follow the neo-Hookean behavior, with the constant $C_{10}$ representing the material stiffness (Equation (2)).

$$\psi =C_{10}(I_1-3)$$

(2)

$C_{10}$ can be determined from the Young’s modulus ($E$) and the Poisson ratio ($\upsilon$) using $C_{10}=E/\left(4\right(1+\upsilon \left)\right)$. The Cauchy stress tensor ($\mathit{\sigma }$) was calculated for the myocardium as shown in Equation (3). Additionally, incompressibility was maintained by using $p$ as the Lagrange multiplier.

$$\mathit{\sigma }=\sum _{k=\mathrm{1,4}\mathrm{f},4\mathrm{s},8\mathrm{f}\mathrm{s}}{\displaystyle \frac{\partial \psi }{\partial I_k}}{\displaystyle \frac{\partial I_k}{\partial \mathit{F}}}\cdot {\mathit{F}}^{\mathrm{t}}-p\mathit{I}$$

(3)

The deformation gradient and identity tensors are represented by $\mathit{F}$ and $\mathit{I}$, respectively.

The constitutive behavior of the myocardium was defined utilizing the UANISOHYPER_INV user subroutine [49,50] in the ABAQUS FEA software (Simulia, RI, USA).

2.3. Coupled Shape-Material Estimations

Shape and material optimization were combined to determine biomechanical scale factors and the unloaded ventricles based on the RV end-diastolic pressure-volume relationships (EDPVR) and ED configuration. This method was similar to the technique we previously introduced for the LV EDPVR in ischemic cardiomyopathy [29]. Initially, the target RV EDPVR was established for each patient, assuming that the Klotz curve [51] could also serve as the target RV EDPVR [52,53]. Bayesian optimization was employed for the simulated RV EDPVRs to align them as closely as possible with the target curves. This alignment was achieved by using a pair of biomechanical scale factors (as inputs of the objective function) to adjust the base experimental material parameters. By minimizing the output of the target function ($O_{\mathrm{v}}$), which is an amplification of the Normalized Root Mean Squared Deviation of RV volumes (NRMSD_V) by $\beta$, both scale factors and unloaded shapes will be established. Amplifying the volume deviations helped us achieve accurate results from the optimizations.

$$O_{\mathrm{v}}=\beta {\mathrm{N}\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{D}}_{\mathrm{v}}={\displaystyle \frac{\beta }{V_{D-1}^{\mathrm{p}\mathrm{r}\mathrm{e}}-V_0^{\mathrm{p}\mathrm{r}\mathrm{e}}}}\times \sqrt{{\displaystyle \frac{\sum _{i=0}^{D-1}{(V_i^{\mathrm{s}\mathrm{i}\mathrm{m}}-V_i^{\mathrm{p}\mathrm{r}\mathrm{e}})}^2}{D}}}$$

(4)

In Equation (4), the total number of $D$ simulated ($V_i^{\mathrm{s}\mathrm{i}\mathrm{m}}$) and predicted ($V_i^{\mathrm{p}\mathrm{r}\mathrm{e}}$) RV volumes were used to calculate the output. The index $i=D$ represents RVEDV, where the simulated and predicted cavity volumes of the RV were calibrated during the estimation of the unloaded shapes. Therefore, RVEDV ($=V_D$) was not included in Equation (4). The simulated RV EDPVR was generated when the unloaded FEA models of the ventricles were loaded according to the estimated ED pressures. The best pair of scale factors was selected that resulted in the lowest possible NRMSD_V, indicating that the simulated and target RV EDPVRs were closely matched. This was achieved by minimizing the output of Equation (4).

The reference (unloaded) geometry was calculated using the augmented iterative method [54], where the difference between the loaded node configurations at the current iteration $\left({\mathit{x}}^{iter}\right)$ and the target node coordinates at ED (${\mathit{x}}^{\mathrm{E}\mathrm{D}}$) was defined as ${\mathit{R}}^{iter}$. The undeformed configuration at each iteration (${\mathit{X}}^{iter}$) was updated to ${\mathit{X}}^{iter+1}$ according to Equation (5).

$${\mathit{X}}^{iter+1}={\mathit{X}}^{iter}-\gamma {\mathit{R}}^{iter}$$

(5)

Where $\gamma$ is the augmentation parameter used to facilitate the convergence rate. It was initially set to be one and then updated according to Equation (6).

$$\gamma \leftarrow -\gamma {\displaystyle \frac{{\mathit{R}}^{iter-1}:({\mathit{R}}^{iter}-{\mathit{R}}^{iter-1})}{({\mathit{R}}^{iter}-{\mathit{R}}^{iter-1}):({\mathit{R}}^{iter}-{\mathit{R}}^{iter-1})}}$$

(6)

2.4. Addition of the Epicardial Restraint

The effect of the epicardial restraint was incorporated into the unloaded FEA models by adding a layer of linear hybrid (C3D6H) elements. This was achieved by extending the epicardial element faces using the Bottom-Up mesh tool with the offset meshing method. The offset value was equal to the thickness of the epicardial constraint. As a result, all the added elements shared common nodes with the ventricular elements at the epicardial surface.

The restraint causes the RV EDPVR to shift to the left by reducing the RV volume. It also decreases the unconstrained (original) ${\mathrm{R}\mathrm{V}\mathrm{E}\mathrm{D}\mathrm{V}}_{\mathrm{o}\mathrm{r}\mathrm{i}\mathrm{g}}$ to ${\mathrm{R}\mathrm{V}\mathrm{E}\mathrm{D}\mathrm{V}}_{\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}\mathrm{e}\mathrm{d}}$ through a reduction ratio ($r$) described in Equation (7).

$${\mathrm{R}\mathrm{V}\mathrm{E}\mathrm{D}\mathrm{V}}_{\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}\mathrm{e}\mathrm{d}}=V_{\mathrm{u}\mathrm{n}\mathrm{l}\mathrm{o}\mathrm{a}\mathrm{d}\mathrm{e}\mathrm{d}}+\left({\mathrm{R}\mathrm{V}\mathrm{E}\mathrm{D}\mathrm{V}}_{\mathrm{o}\mathrm{r}\mathrm{i}\mathrm{g}}-V_{\mathrm{u}\mathrm{n}\mathrm{l}\mathrm{o}\mathrm{a}\mathrm{d}\mathrm{e}\mathrm{d}}\right)\times \left(1-{\displaystyle \frac{r}{100}}\right)$$

(7)

It is noteworthy that “$r$” is defined in a way that reduces the volume change of the RV EDPVR $\left({\mathrm{R}\mathrm{V}\mathrm{E}\mathrm{D}\mathrm{V}}_{\mathrm{o}\mathrm{r}\mathrm{i}\mathrm{g}}-V_{\mathrm{u}\mathrm{n}\mathrm{l}\mathrm{o}\mathrm{a}\mathrm{d}\mathrm{e}\mathrm{d}}\right)$ and not just ${\mathrm{R}\mathrm{V}\mathrm{E}\mathrm{D}\mathrm{V}}_{\mathrm{o}\mathrm{r}\mathrm{i}\mathrm{g}}$. Therefore, when $r=100$%, Equation (7) results in zero volume change in RV EDPVR (${\mathrm{R}\mathrm{V}\mathrm{E}\mathrm{D}\mathrm{V}}_{\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}\mathrm{e}\mathrm{d}}=V_{\mathrm{u}\mathrm{n}\mathrm{l}\mathrm{o}\mathrm{a}\mathrm{d}\mathrm{e}\mathrm{d}}$).

We defined the target restrained RV EDPVR assuming two different reduction methods, both of which maintain the nonlinearity of the shifted curve. Firstly, we considered an exponential volume reduction of the RV. The RV EDPVR curve of the restrained RV was determined using a scaling scheme that redistributes the volumes between the unloaded ($V_{\mathrm{u}\mathrm{n}\mathrm{l}\mathrm{o}\mathrm{a}\mathrm{d}\mathrm{e}\mathrm{d}}$) and the previously moved (${\mathrm{R}\mathrm{V}\mathrm{E}\mathrm{D}\mathrm{V}}_{\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}\mathrm{e}\mathrm{d}}$) volumes as described in Equation (8).

$$V_{\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}\mathrm{e}\mathrm{d}}(n)=V_{\mathrm{u}\mathrm{n}\mathrm{l}\mathrm{o}\mathrm{a}\mathrm{d}\mathrm{e}\mathrm{d}}\times \mathrm{e}\mathrm{x}\mathrm{p}\left(\mathrm{l}\mathrm{n}\left({\displaystyle \frac{{\mathrm{R}\mathrm{V}\mathrm{E}\mathrm{D}\mathrm{V}}_{\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}\mathrm{e}\mathrm{d}}}{V_{\mathrm{u}\mathrm{n}\mathrm{l}\mathrm{o}\mathrm{a}\mathrm{d}\mathrm{e}\mathrm{d}}}}\right)\times {\displaystyle \frac{n}{\mathrm{N}-1}}\right)$$

(8)

$n$ changes from zero to $\mathrm{N}-1$, where $\mathrm{N}$ is the total number of available PV data points that make up the RV EDPVR ($\mathrm{N}>1$). Figure 1a illustrates the effect of exponential reduction on the target RV EDPVR with varying values of $r$. It is important to recognize that these reduced RV EDPVRs will serve as the target curves for optimizing the RV constraints.

Equation (9) describes the second reduction strategy for the target RV EDPVR, which involves proportionally reducing all volumes of the unconstrained RV EDPVR ($V_{\mathrm{o}\mathrm{r}\mathrm{i}\mathrm{g}}\left(n\right)$). As shown in Figure 1b, this reduction results in a smoother and more systematic progression, with consistent spacing between the original and reduced curves.

$$V_{\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}\mathrm{e}\mathrm{d}}\left(n\right)=V_{\mathrm{o}\mathrm{r}\mathrm{i}\mathrm{g}}\left(n\right)-({\mathrm{R}\mathrm{V}\mathrm{E}\mathrm{D}\mathrm{V}}_{\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}\mathrm{e}\mathrm{d}}-V_{\mathrm{u}\mathrm{n}\mathrm{l}\mathrm{o}\mathrm{a}\mathrm{d}\mathrm{e}\mathrm{d}})\times \left({\displaystyle \frac{n}{\mathrm{N}-1}}\right)\times \left({\displaystyle \frac{r}{100-r}}\right)$$

(9)

Bayesian optimization was used not only to estimate material scale factors but also to determine the optimal thickness of the RV epicardial restraints. During the optimization of the restraints, the output of the target function was defined similarly to Equation (4) with $V_i^{\mathrm{p}\mathrm{r}\mathrm{e}}$ being determined according to the target volumes of either Equation (8) (exponential reduction) or Equation (9) (smooth or proportional reduction).

20 initial points were generated through the Latin hypercube sampling (LHS) method. The acquisition function was probabilistically selected from the lower confidence bound (LCB), the expected improvement (EI), or the probability of improvement (PI) at each iteration. The optimization stopped successfully when the change in $O_{\mathrm{v}}$ was either below 0.1 in 15 iterations or below 0.01 in 5 iterations, with a maximum allowable number of 200 function evaluations.

When the convergence criteria remain constant, increasing the amplification coefficient ($\beta$) provides more accurate results at the expense of more iterations during Bayesian optimization. To achieve a reliable fit and relatively fast convergence, $\beta$ was selected to be 10000. The choice of the amplification depends on the convergence criteria and it is possible to opt for a lower $\beta$ together with stricter convergence criteria.

The Scikit-Optimize Python package [55] was used to perform Bayesian optimizations in the present study. We determined the optimal restraint thickness for seven different reduction values ($r=$ 10%, 15%, 20%, 25%, 30%, 35%, and 40%) and for eight different restraint material stiffness values ($C_{10}=$ 0.005, 0.01, 0.02, 0.03, 0.05, 0.1, 0.3, and 0.5 MPa). The mentioned stiffness range came from the reported passive mechanical properties of the epicardial restraints in the literature [21,28,56]. This resulted in a total of 56 restraint optimizations per patient.

2.5. Sensitivity Analysis

The sensitivity of the optimized restraint thickness to the choice of initial evaluations was analyzed by comparing the results obtained from four random sets of initial points. Additionally, we examined the changes in the RV EDPVRs with a wide range of restraint thicknesses and materials. The selected thickness range for this purpose was 0.1 to 12 mm, while $C_{10}$ (representative of restraint stiffness) varied between 0.01 and 10 MPa.

2.6. Computational Framework in Brief

Our computational design method has the capability of being implemented in hospitals to provide optimal RV restraint therapy. The patient should initially undergo medical imaging (CMR and echocardiography). The extracted geometry and conditions will then be used to create an FEA model of the ventricles at ED. A combined material and shape optimization process will be used to estimate the mechanical properties and reference configurations simultaneously. Personalized virtual epicardial restraint will then be determined in order to achieve a target passive filling that guarantees shifting the RVEDV to a desired curve. The optimal restraint can be prepared and applied to the patient during cardiac surgery at hospitals (see Figure 2).

3. Results

3.1. Computational Performance

Simulations were performed on an HP Z440 Workstation equipped with an Intel^® Xeon^® E5-2680 v4 CPU, which has 14 cores with a maximum turbo frequency of 3.30 GHz, 32 GB DDR4 RAM, and an NVMe M.2 SSD running the Linux operating system. To improve resource usage, two simulations were run simultaneously, with half of the CPU cores dedicated to each one. Table 2 compares the average timing, number of iterations, estimated material scale factors, and the goodness of fit for each restraint optimization.

3.2. Sensitivity Results

To assess the robustness of Bayesian optimizations with respect to the selection of initial function evaluations, we compared the optimization results obtained from four separate random sets of initial points for both the first and second patients in all 56 cases (as detailed in section 0). In 98.8% of simulations for the first patient and 99.3% of simulations for the second patient, the coefficient of variation was below 1%.

Figure 3 demonstrates the sensitivity of the second patient’s RV EDPVR to the restraint thickness with a $C_{10}$ of 0.03 MPa (Figure 3a) and to the restraint stiffness with a thickness of 0.5 mm (Figure 3b). The black and yellow dashed lines show the target RV EDPVR and its exponentially reduced curve with $r=30\%$, respectively. The target value for $V_{\mathrm{u}\mathrm{n}\mathrm{l}\mathrm{o}\mathrm{a}\mathrm{d}\mathrm{e}\mathrm{d}}$ was 48.8 ml, which was calculated from the predictions of [51]. Using Equation (7) with a target $r=30\%$, the goal RVEDV for the second patient experienced a decrease from 90.3 ml to 77.85 ml. By thickening the RV restraint from zero to 12 mm at a constant stiffness ($C_{10}=0.03MPa)$, the simulated RVEDV decreased from 91.33 ml to 71.47 ml with an unchanged $V_{\mathrm{u}\mathrm{n}\mathrm{l}\mathrm{o}\mathrm{a}\mathrm{d}\mathrm{e}\mathrm{d}}$ of 48.36 ml (Table 2 and Figure 3a). From Equation (7), the $r$ became 46.22% considering volume changes of FEA-based RV under these conditions. Similarly, when the constraint was stiffened by increasing $C_{10}$ to 10 MPa while its thickness remained constant at 0.5 mm, the simulated RVEDV declined to 68 ml resulting in $r=54.29\%$ (Figure 3b). As shown in Figure 3, a similar decrease was also evident for the whole RV EDPVR except for the unloaded volume. This volume reduction continues until the reduced figures converge to specific curves in both figures.

3.3. RV EDPVRs Before Restraint

Compared to the target RV EDPVRs, the simulated curves, produced from shape-material optimizations, showed good agreement before the addition of restraints, i.e., $r=0\%$ or $t=0$ (Figure 4).

3.4. Optimal Restraints

The 3D surface plots in Figure 5 illustrate the relationship between restraint stiffness, the expected reduction in RVEDV, and the estimated optimal thickness of the restraint. The plots also display the NRMSD_V associated with each optimal value. It was observed that, at a constant reduction ratio, the optimized restraint thickness decreased as $C_{10}$ increased. Conversely, for a specific type of restraint, a higher optimized thickness was achieved for a higher desired $r$. While the surface plots exhibited similar trends across all cases, the optimized thicknesses and NRMSD_V values varied for each patient. Notably, the third patient (with an RVEDP of 5mmHg) had a higher NRMSD_V. A low NRMSD_V (approaching zero) indicates a perfect fit, whereas high values suggest a significant deviation between the goal and achieved RV EDPVRs. Interestingly, all patients in the study demonstrated lower NRMSD_V in the second row, which is related to the exponential reduction in target curves.

Fiber Green strain ($E_{\mathrm{f}\mathrm{f}}$) and fiber Cauchy stress (${\sigma }_{\mathrm{f}\mathrm{f}}$) are compared in Figure 6 for different desired $r$ values in a representative biventricular model of the second patient assuming $C_{10}=0.03$ MPa. In general, a higher $r$ value requires a greater restraint thickness with a lower $E_{\mathrm{f}\mathrm{f}}$ and ${\sigma }_{\mathrm{f}\mathrm{f}}$ in the RVFW.

4. Discussion

We have developed a computer-based design method for determining the optimal value of epicardial passive restraint in patients with RV failure. Our Bayesian optimization method demonstrates high robustness, regardless of the initial points selected. Low NRMSD_V values were achieved for all material estimation simulations, indicating an excellent fit for the ventricular material scale factor estimations (see Table 2 and Figure 4). A more accurate fit has already been obtained for ischemic ventricles assuming LV EDPVR as the target [29]. In addition, the scale factors in healthy regions of the myocardium in that study were considerably lower than those found in the present study based on the RV EDPVR. However, due to the necessity of optimizing additional ischemic parameters in that study, the average total solution time for shape-material estimations here was approximately half of those reported previously for a biventricular model running on a relatively similar PC [29] (see Table 2). It is worth mentioning that the computational time depends on the complexity of the models and the available computational resources. Additionally, the shape-material estimation times were higher than the fitting run-times of the LV models reported for supercomputer resources with 96 CPU cores [31] and significantly lower than those reported for an LV being simulated on a regular workstation [57]. Among the three RV failure patients included in our study, the third patient had a longer running time. This patient’s FEA model contained the largest number of elements and nodes compared to the others.

The constraint sensitivity analysis revealed that while increased restraint thickness leads to reduced RV passive volumes, the rate of this decline is not consistent. Finally, the RV EDPVR converges to a final curve, after which increased thickness does not result in reduced RV volumes. Similar leftward shifts have already been reported in the literature using FEA to model the Acorn CorCap Cardiac Support Device (Acorn Cardiovascular Inc, St. Paul, MN) on a biventricular model [58], a constrained LV with pressure on the epicardium [2], and a longitudinal patch added to the epicardium of the acutely ischemic LV [16]. However, small changes in LV EDPVR have been reported after the addition of a longitudinal patch [16]. Another experimental study on dogs showed that longitudinal support of acute anteroapical infarcts does not affect LV EDPVR while improving cardiac output (CO) [59]. In addition, computational modeling showed that isotropic stiffening of the aforementioned infarcts moves LV EDPVR to the left and improves contractility, although it does not substantially influence stroke volume [17]. Additionally, it was reported that the RVEDV did not significantly alter when higher levels of quantitative ventricular restraint were applied solely to the RV epicardium [4]. The present study showed that constraint thickness and stiffness determine the degree of such ventricular volume reductions if passive constraints are applied to the epicardium of the unloaded RV (Figure 3).

Surface plots from Figure 5 indicate that the optimal restraint thickness is dependent on RVEDP, with patients having lower RVEDP requiring slimmer restraints. These surface plots also demonstrate that optimized restraint values provide a better fit for the exponentially reduced target RV EDPVRs compared to the proportional ones, particularly at higher RVEDP levels. Also, it is not advisable to use stiff restraints for low desired $r$ values due to the high NRMSD_V, as indicated by the lower reddish corners of the surface plots in Figure 5. It is noteworthy that these conclusions are based on Klotz’s EDPVR as the goal curve for the RV passive filling, and there is potential for further exploration in future studies to assess the impact of measured target RV EDPVRs as opposed to the empirical Klotz curves utilized in this study.

Figure 6 illustrates a reduction in fiber stress and strain as restraint thickness increases for patient 2, with similar findings observed in the other two cases. Previous studies have shown decreased LV fiber stress at ED following the implementation of epicardial constraints [2,58]. In addition, restrained LV has been associated with reduced ED fiber strain [16], while isotropic reinforcement substantially decreased both circumferential and longitudinal ED strains [17].

According to the law of Laplace, reducing RV diameter and increasing its wall thickness after restraint both lead to decreased wall stress. This is similar to the changes in fiber stress presented in Figure 6. While this law can offer valuable insights into stress variations, it may not always provide precise estimations of LV stress changes following the treatment of dilated cardiomyopathy using passive devices [60]. The superiority of FEA over Laplace's law in stress calculations for the LV has been well-documented [61]. FEA not only predicts stress and strain alterations, but it can also determine their spatial distributions across the ventricles.

Limitations and Suggestions

Despite all the mentioned benefits, the methodology proposed in this research may have a few downsides. Our FEA and optimization were conducted to assess mechanical support provided by isotropic epicardial patches for the enlarged RV. The study focused on passive restraining, excluding the evaluation of therapeutic ingredients and healing effects [62]. While incorporating systolic behavior is crucial for understanding the impact of constraint thickness on CO, our modeling primarily centered on the diastolic filling of the RV. Also, the method utilized in this study is only able to predict the acute effects of restraints on passive filling (RV EDPVR), assuming consistent RVEDP before and after stiffening, and disregarding RVEDP changes [4,63]. Our assumption of a smooth variation of myofibers in the ventricular wall may have influenced the results negatively, given the structural changes and remodeling in HF. A sensitivity analysis has been conducted previously to demonstrate the effect of fiber orientations on estimated material parameters in ischemic ventricles [29]. For more accurate results, detailed whole-organ models should be considered; however, this may significantly increase solution times. The uniqueness of the estimated unloaded shapes is not yet guaranteed [31]. An additional limiting assumption could be the lack of consideration for prestretch when simulating epicardial restraints. While prestretch may be observed when passive biventricular restraints cover both ventricles, we do not anticipate high levels of prestretch in RV-specific applications. In addition to the aforementioned declarations, there are a few more limitations to our material estimations and unloading method that have already been discussed [29].

5. Conclusions

A novel and easy-to-implement computational framework has been proposed in the present study to predict the passive effect of constraining the RV epicardial surface noninvasively. Our robust method not only successfully estimated both the reference stress-free configuration and the passive biomechanical properties based on experiments and medical imaging, but it also found the optimal thickness for restraining the RV to desired volume reduction ratios using a wide range of available restraint materials. Low NRMSD_V values were achieved for both material estimations and restraint optimizations, demonstrating the accuracy of the results. It was observed that the exponentially reduced target RV EDPVR can better mimic the behavior of the restrained ventricle. A sensitivity study revealed that increasing constraint thickness and stiffness results in a leftward shift of the RV EDPVR. As thickness and stiffness rise, the altered curves would eventually converge to a final state, beyond which further stiffening and thickening could not shift the curve to the left.

Patient-specific surface plots were generated for each patient to show the best thickness of a wide range of available constraint materials for different desired volume reduction values. With a rise in the target reduction, the optimal thickness increased, while fiber stress and strain decreased accordingly. One advantage of the current study is the consideration of the unloaded ventricular shape, personalized material estimations based on the RV EDPVR, and valvular effects. Although the primary focus of this research was modeling the passive constraints of the RV, similar techniques can also be applied to explore the effect of epicardial restraints on the passive behavior of the LV.