Computational Design of a Soft Robotic Myocardium for Biomimetic Motion and Function

Soft robotic devices containing multiple actuating elements have successfully recapitulated complex biological motion, leading to their utility in biomedical applications. However, there are inherent nonlinear mechanics associated with soft composite materials where soft actuators are embedded in elastomeric matrices. Predicting their overall behavior prior to fabrication and subsequent experimental characterization can therefore present a hurdle in the design process and in efficiently satisfying functional requirements and specifications. In this work, a computational design framework for optimizing the motion and function of biomimetic soft robotic composites is demonstrated by conducting a design case study of soft robotic cardiac muscle (myocardium) with a particular focus on applications including replicating and assisting cardiac motion and function. A finite element model of a soft robotic myocardium is built, in which actuators are prescribed with anisotropic strain to simulate local deformation, and various design parameters are investigated by evaluating the performance of each configuration in terms of ventricular twist, volumetric output, and pressure generation. Then, an optimized design is proposed that recapitulates the physiological motion and hemodynamics of the heart, and its thrombogenicity is further explored using a fluid‐structure interaction model. This framework has broader utility in predicting the performance of other soft robotic embedded composites.


Introduction
With the emergence of soft robotics over the past decade, there has been significant development of biomimetic robots using soft materials, offering life-like motion and safe interaction with humans. Nature provides excellent inspiration for many soft robots with unique functions, ranging from spiders [1] and fish [2] to worms. [3] In the biomedical field, soft devices comprising multiple actuating elements have been developed for prosthetics [4,5] or assist devices (including gloves, [6,7] exoskeletons, [8,9] cardiac assist devices [10,11] ), body-part simulators [12][13][14][15] and artificial organs. [16,17] For cardiovascular applications, traditional devices for mechanical circulatory support systems (i.e., ventricular assist devices (VADs)) use pumps with rigid components [18] to produce centrifugal (e.g., HeartWare HVAD (HeartWare), HeartMate III (Thoratec)) or axial flow (e.g., HeartAssist 5 (ReliantHeart), Heart-Mate II (Thoratec), Impella (Abiomed)) to replace cardiac function. Benchtop cardiac models use positive displacement pumps to recapitulate the function of the heart. [19][20][21] While these approaches may be effective in improving or recreating blood flow, the imparted biomechanics often differ from those of the native heart. High-speed rotary pumps or continuous flow VADs are associated with high shear stress and non-physiologic flows, leading to complications including gastrointestinal bleeding, arteriovenous malformations, hemolysis, and pump thrombosis. [22][23][24] In vitro simulators that replicate dynamic flow often suffer from inverse pressure-volume relationships, and unnatural hemodynamics and flow patterns. [25] Devices that reinstate physiologic pressure-volume relationships by external compression of the heart have been developed. Extracardiac assist devices, such as Car-dioSupport System (Cardio Technologies Inc), Heart Booster (AbioMed Inc), HeartPatch (Heart Assist Technologies), CorInnova (CorInnova Inc), and BEAT (AdjuCor), are pneumatically or hydraulically powered to inflate pockets surrounding the heart and cyclically compress the heart to provide partial or biventricular assistance. Some total artificial hearts (TAHs), such as SynCardia TAH (SynCardia Inc) and soft TAH [16,17] have been developed using similar mechanisms.
Recently, biomimicry and bioinspired design approaches using soft robotic techniques have been proposed by our group and other groups [10,14,15] to better replicate cardiac motion. Inspired by the fiber architecture of the native cardiac muscles, these macroscale soft robotic cardiac devices have multiple soft actuators analogous to biological muscles that undergo anisotropic motion (axial contraction and radial expansion) embedded in an elastomeric matrix. By spatially arranging the actuators along muscle fiber orientations, soft robotic devices achieve complex cardiac motion, such as compression and twisting, upon actuation, which is challenging to achieve with rigid robots. Yet, these do not replicate the hemodynamic function of the heart in its entirety. Several challenges limit their broader use for medical applications in which specific functional requirements must be met. First, optimization of the overall behavior of the device prior to fabrication and testing due is challenging owing to the nonlinear mechanics inherent to soft materials and complex interactions between multiple soft components. Second, these devices often have long manufacturing processes involving multi-step casting, making it difficult to troubleshoot or rapidly iterate through multiple designs. Third, the design space is often very broad due to the diverse choices of actuators and matrices and multiple combinations of possible spatial architecture. Therefore, for applications with specific, predetermined functional requirements, there is a need for a systematic approach and framework to design and optimization. [26] In this work, we demonstrate a computational design framework for optimization of motion and hemodynamics of a biomimetic soft robotic myocardium towards benchtop or total artificial heart applications ( Figure 1). First, we build a finite element (FE) model composed of linearly contracting soft actuators arranged in a soft matrix in an idealized left ventricle (LV) geometry. The actuators are prescribed an orthotropic strain to simulate the local deformation, and the resulting global behavior of the soft robotic myocardium is evaluated computationally. A variety of design parameters including architecture, material stiffness and type of soft actuators are varied by measuring the performance of each configuration with a specific focus on left ventricular twist, volumetric output, and pressure generation. Based on the understanding of the design variables, we sought to merge biomimicry with mechanical design principles to maximize biological relevance while preserving mechanical efficiency and simplicity. We propose an optimized design of soft robotic myocardium for replication of physiological motion and hemodynamics. We simulate left ventricular twist of 8.5°, ejection fraction of 68.5% and systolic pressure generation of up to 121 mmHg, which are within the physiological ranges. Further, we perform a two-way fluid-structure interaction (FSI) method by coupling FE and computational fluid dynamics (CFD) analysis. We show that a simplified design with reduced biomimicry to previous work [10,14,15] is sufficient to recapitulate the motion and function of the heart. This framework has broader utility in predicting the performance of other soft robotic embedded composites.

Design Requirements
The scope of this case study is to design a soft robotic myocardium that closely captures the biomechanics of the healthy LV during the systolic phase. The LV of a healthy individual generates an intraventricular pressure of 100-140 mmHg during the systolic phase ejecting ≈50-70% of ejection fraction by squeezing and twisting the cardiac muscles. A small amount of twisting motion (physiological range: 11.2 ± 4.1° of total rotation [27] ) is important for pumping efficiency and also ideal to avoid blood stagnation.
In this work, we aim to meet the following design requirements for soft robotic myocardium: 1) Left ventricular ejection fraction of 50% or above 2) Left ventricular pressure generation of at least 100 mmHg 3) Left ventricular twisting motion (apical rotation) of 7 to 15°d uring the systolic phase, given an idealized LV geometry in an ellipsoidal shape with an average myocardial wall thickness of 13.4 mm.

Design Variables
Users have numerous options in selecting the type of actuators, and soft matrix materials. Some common soft contracting actuators with fluidic actuation include McKibben pneumatic artificial muscles (PAMs), [28] flat pleated PAMs, [29,30] fluid-driven origami-inspired artificial muscles (FOAMs), [31] and vacuum-actuated muscle-inspired pneumatic structures Figure 1. Overview of design process of soft robotic composites for cardiac applications. Image of the fiber tractography of the ex vivo porcine heart is readapted with permission from Froeling et al., 2022, Springer Nature. [48] Image of the mechanical pulsatile pump is readapted from Köhne et al., and image of the soft total artificial heart is readapted with permission from Cohrs et al., 2022, John Wiley and Sons. [16,49] (VAMPs). [32] Others include dielectric elastomer actuators, [33,34] thermal actuators, [35] chemically responsive actuators, and [36] magneto-responsive actuators. [37] All may vary in size (i.e., initial dimensions), deformation characteristics (i.e., degree of expansion and contraction), and the amount of force generation relative to the amount of contraction (i.e., active stiffness) at the individual actuator unit level. The spatial distribution and placement of these actuators add to the complexity to design choices. There are also a wide variety of soft elastomeric materials to choose from. Some common soft matrix materials include a low modulus polyvinyl alcohol hydrogel, low modulus Ecoflex, and higher moduli SmoothSil silicones, which affect the motion of actuators and corresponding force transmission through the medium. In some cases, addition of reinforcement or restraining layers can induce motion in a preferential direction. These properties, in aggregate, determine the overall motion and deformation of the soft robotic myocardium translating to fluid left ventricular volume reduction (ΔLVV) and left ventricular pressure generation (ΔLVP) in the chamber.
As a result, the design space is wide with a vast number of possible solutions and therefore it is important to understand the contribution of individual design variables to meet the design requirements for a soft robotic myocardium. The major design variables investigated in this study are outlined here (

Effect of Actuator Spatial Density
Unlike the native myocardium, which is constituted by a continuum of myocardial fibers, the soft robotic myocardium consists of a discrete number of active elements embedded in a soft matrix. To examine the contribution of actuator spatial density in a given matrix to the overall behavior of the soft robotic myocardium, we varied the number of actuators in the circumferential and helical arrangements (Figure 3). Varying the number of actuators in the soft robotic myocardium led to changes not only in the amount of ΔLVV and motion but also in the degree of ΔLVP in the soft robotic myocardium. When actuators were circumferentially arranged, adding more actuators resulted in a linear increase in both the ΔLVV and ΔLVP of the soft robotic myocardium ( Figure 3A).
In contrast, opposite trends were observed in helical arrangements, where actuators were oriented at 60° from the basal plane (as described in Roche et al.) [15] with lower ΔLVV and ΔLVP as more actuators were being added ( Figure 3B). The difference may be attributed to the differences in the motion constraints on the actuator radial expansion. Circumferentially arranged actuators produced limited outward motion from the longitudinal axis despite actuator expansion since the actuators are bound in ring-like configurations undergoing circumferential shortening while helically arranged actuators had the degree of freedom to expand in directions that minimize the strain energy. With additional units, helically arranged actuators favored outward motion to avoid crowding and produced less ΔLVV and ΔLVP. However, adding actuators contributed to more subapical rotation in the helical arrangement ( Figure 3B).

Effect of Actuator Motion Characteristics
Soft contracting actuators have unique deformation properties with a wide spectrum of combinations for axial contraction and radial expansion (Figure 4A,B). To understand how the changes in the actuator deformation geometry during actuation affect the global configuration of the soft robotic myocardium, we performed a parametric study varying the amount of axial and radial strains of actuators in a representative design ( Figure 4C) in our FE model. Initial geometry and active material stiffness of a standard McKibben actuator were held constant across the tests for a comparative analysis focusing on the geometric contribution from individual actuators on the overall motion of the soft robotic myocardium. Figure 4D shows that when only contraction or radial expansion were applied to actuators, the soft robotic myocardium was not able to achieve a physiological level of ΔLVV, regardless of the degree of contraction or expansion. This limitation was overcome by combining contraction and radial expansion, resulting in a physiological level of ΔLVV. In this simple circumferential configuration, actuator contraction motion was more effective than expansion because the actuator expansion not only added to inward motion but also longitudinal elongation of the LV, countering the effect of ΔLVV in the radial direction.

Effect of Material Stiffness
When designing a soft robotic myocardium, there are ranges of actuator types and soft matrices with varied material stiffness to choose from. For actuators, there may be a wide range of active material stiffness (defined here as stiffness in active contraction measured in the axial direction) even amongst McKibben actuators depending on the initial diameter and inner bladder materials. Flat pleated pneumatic artificial muscles (fPPAMs), a different type of soft contracting actuator, have about 30 times lower active material stiffness than the 3-mm initial diameter McKibben actuator that was used as our baseline actuator. Similarly for the matrix material, choices are ranging from hydrogel, low-modulus to high-modulus elastomers with up to a 70-fold difference in stiffness.
To better understand how these material properties determine the overall behavior of soft robotic myocardium, we studied the effect of actuator and matrix material stiffnesses in a circumferential design with three actuators. The parametric sweep was conducted across a wide range of stiffness values while holding the actuator deformation characteristics constant as the baseline actuator. In the lower range of active material stiffness values for actuators, increasing the stiffness resulted in the increase of both ΔLVV and ΔLVP, but the effect lessened beyond the baseline stiffness ( Figure 5A,B).
Interestingly for the matrix, changing the material stiffness significantly increased the ΔLVP with negligible effect on the ΔLVV ( Figure 5C,D). The force generated from actuator contraction is transmitted more efficiently in a stiffer medium than in a soft medium, which absorbs the energy. This indicates that ΔLVP may be modulated independent of ΔLVV by increasing the matrix stiffness when optimizing for hemodynamic performance.

Effect of Passive Restraint Layer
In cardiac anatomy, there is a stiff fibrous membrane called the pericardium that encloses the ventricular myocardium  ( Figure 6A). It surrounds the epicardial surface of the ventricle and limits volume expansion, affecting the hemodynamics in both healthy and diseased states. Inspired by nature, a passive restraint layer with high material stiffness was applied on the surface of the soft robotic myocardium with the circumferential design to assess whether function could be augmented for the simplest design. In our computational model, a 0.3 mm-thickness thermoplastic urethane (Young's modulus = 18.8 MPa) sheet was selected to model a synthetic substitute for human pericardial tissue (Young's modulus = 20.4 ± 1.9 MPa [38] ), and a rigid outer shell of 2 mm thickness (Young's modulus = 2.5 GPa) was chosen to model a more extreme case of inextensible restraint. As expected, we found that by passively restraining the epicardial surface, the soft robotic myocardium preferentially produces more inwards motion towards the central long axis increasing the ΔLVV and improving ΔLVP ( Figure 6B) by mitigating the LV elongation and outward radial expansion observed without epicardial restraint ( Figure 6C).

Bioinspired Design
A bioinspired design replicating the fiber architecture of the native LV myocardium was first considered (Figure 7A). Given the geometric constraint, three layers of actuators were placed transmurally, oriented helically in the endocardial ("endo") layer (α endo = 60°), circumferentially in the middle ("mid") layer (α circ = 0°) and helically in the epicardial ("epi") layer (α epi = −60°). In each layer, three actuators were used to generate desired motion based on the results discussed in Section 3.1. To identify the dominant contribution from actuator placement, FE models with different combinations of layers ("endo"/"mid"/"epi") were built and examined for their motion and function as summarized in Figure 7B,C and Table S1 (Supporting Information).
Among the single-layer arrangements, the circumferential arrangement ("mid only") was the most efficient in terms of    ΔLVV and ΔLVP by driving circumferential shortening. Helical arrangements in both endocardial ("endo only") and epicardial layers ("epi only") contributed to large degrees of subapical rotations but failed to induce adequate ΔLVV or ΔLVP. While "endo only" and "epi only" were symmetrically about the axis, placing the actuators in the endocardial space ("endo only") was notably more efficient compared to placing them farther away in the epicardial region ("epi only") since in the latter, the radial expansion of actuators contributed to an outward motion away from the axis rather than an inward motion increasing ΔLVV.
For double actuator layer arrangements, the combination of helical actuators in the endocardial and epicardial layers ("endo + epi") demonstrated enhanced effects compared to the summed performance of the individual endocardial and epicardial single-layers ("endo only" + "epi only"). Here, the double helical actuators ("endo + epi") contributed to the shortening of the LV in the longitudinal direction rather than just circumferential contraction as demonstrated in the circumferential arrangement in the single-layer ("mid only"). As expected, including the circumferential actuators was beneficial and efficiently augmented the amount of ΔLVV and ΔLVP, achieving values within the physiological range ('"endo + mid", " mid + epi"). With the addition of the circumferential layer to either the endocardial or epicardial layer ("endo + mid", " mid + epi"), the resulting volume or pressure changes were more pronounced than the summation of the single helical and circumferential arrangements ("endo only" + "mid only" or "epi only" + "mid only"). It is likely that the lengthening of the LV in the circumferential arrangements, which results from the radial expansion of McKibben actuators, was mitigated by the longitudinal shortening provided by the helical actuators.
Consistent with the observations of the helical arrangements in single-layers, helical actuator placement in the endocardial layer ("endo" + "mid") was far more effective than placement in the epicardial layer ("mid" + "epi") when combined with circumferential actuators, resulting in ΔLVV of 85.2% and ΔLVP of 256.5 mmHg. The simulated ΔLVV and ΔLVP values were not very different when all three layers were combined ("endo + mid + epi"), indicating that there was very little contribution from the helical actuators in the epicardial layer in this arrangement. In fact, "endo + mid" produced a subapical rotation of 15.9°, which is much closer to the physiological range while "mid + epi" and "endo + mid + epi" produced a supraphysiological subapical rotation of 43.8° and 34.1°, respectively.

Simplified 2-Layer Design
Based on our understanding of the contribution of each layer (as described in the bioinspired design), a simplified design  containing only two layers with circumferential and endocardial actuators was modeled ( Figure 8A) and the resulting motion and function are summarized in Figure 8B,C and Table S2 (Supporting Information). To generate an overall rotation in the counterclockwise direction in the absence of epicardial actuators, we placed helical actuators in the endocardial layer at a helix angle about the long axis that would be more typical for the epicardial layer. (α endo < 0). α endo was varied to adjust the amount of rotational motion. Varying the helical angles in the endocardial layer maintained a similar level of ΔLVV and ΔLVP compared to "endo" + "mid" in the bioinspired threelayer design. However, there was a wide range of subapical rotation generated, ranging from 22.4 to 6.59° with less rotation induced for endocardial actuators with a higher helix angle (closer to a longitudinal orientation). In the simplified two-layer design, the overall amount of twisting may be easily modulated by changing the α endo , suggesting flexibility in the design.

Design Candidates using Combinations of Different Actuator Types
Based on our understanding of the design parameters of the soft robotic myocardium, we generated multiple design candidates with combinations of different actuator types. Figure 9 shows a summary of the candidates organized by the number of layers. Combinatorial designs were explored using different actuator types with orthotropic expansion properties, such as flattened McKibbens and fPPAMs, which were incorporated in the orientation intended to produce the desired motion and function. Other actuator types were modeled with their approximate initial and final dimensions, deformation characteristics and material properties, as 3D McKibbens were modeled previously (See Figure S1, Supporting Information). Figure 9A-C shows the simplest configurations, a singlelayer of circumferential arrangement with three actuators, consisting of cylindrical McKibbens (3DMcKs), flattened McKibbens (2DMcKs), and fPPAMs, respectively. A reasonable amount of ΔLVV and ΔLVP, with 2DMcKs being most efficient by inducing radial displacement and minimizing longitudinal extension compared to 3DMcKs or fPPAMs that are highly expanding in both directions. However, there is no twisting motion due to the axisymmetric nature of the design. Addition of helical actuators in the endocardial layer at reversed angle (placement similar to Figure 8A) to the circumferential actuators ( Figure 9D-F) generated rotational motions near the physiological level while also significantly increasing ΔLVV and ΔLVP compared to the single-layer designs ( Figure 9A-C). All dual-layer designs were capable of displacing a volume greater than physiological ranges and reaching near physiological values for subapical rotation as well; Compared to using only 3DMcKs for both layers  ( Figure 9D), using 2DMcKs for the circumferential layer augmented radial displacement while mitigating the longitudinal elongation ( Figure 9E). This design was also shown to be superior over the most crowded design with triple-layers of 3DMcK actuators (similar to Figure 7A), which induces a comparable amount of ΔLVV and pressure but generates supraphysiological apical rotation of 34.1°.

Hemodynamic Evaluation and Validation
To evaluate the hemodynamic performance of the design candidate (referring to Figure 9E), we simulated the motion of the soft robotic myocardium under physiological aortic pressure boundary conditions with real-time dynamics. Figure 10A and Video S1 (Supporting Information) show the structural deformation of the optimized design candidate obtained from the FE model, demonstrating physiological ΔLVV and subapical rotation of 8.5°. A two-way FSI model integrating the structural simulation and particle-based computational fluid dynamics simulation was performed to predict the hemodynamic performance of the design during the systolic phase. As a ramp of actuator actuation was applied, the soft robotic myocardium underwent isovolumetric contraction, building up the pressure (i.e., LVP). As the LVP surpassed the boundary pressure (i.e., aortic pressure), it opened the virtual unidirectional valve  placed at the outlet, allowing fluid ejection during that LVV was significantly reduced. A physiological level of hemodynamic performance with an ejection fraction of 68.5% and peak LVP of 121 mmHg was achieved with the optimized design candidate containing a circumferential layer with 2DMcKs and an inner twisting layer with 3DMcKs ( Figure 10B).
Considering the relationship between velocity, shear, and thrombogenicity, [39] we conducted a hemodynamic assessment on thrombogenic potential of the optimized design based on fluid particle velocity and time-averaged wall shear stress (TAWSS) in comparison to our control design ( Figure 9A)-the simplest version of the axisymmetric design with circumferential actuators without helical actuators. The TAWSS index is used because it integrates the velocity field with the low shear flow pattern and identifies places at high risk of thrombus development. [39,40] Figure 11A and Video S2 (Supporting Information) capture a sequence of fluid particle velocity vectors in the two designs, showing a larger amount of fluid ejection for the optimized design due to higher volume displacement and more rapid pressure development compared to the control design.A velocity threshold value of 0.2 m s −1 was prescribed to represent the low-flow regions associated with high thrombogenic risk [41] and the peak-ejection low-velocity particle fraction was determined by the number of fluid particles under the threshold divided by the total number of pre-filled fluid particles. While 65.7% of the total fluid particles were below the low-velocity threshold (0.2 m s −1 ) in the control design, there were only 19.8% of low velocity fluid particles for the optimized design, indicating reduced risk for intraventricular thrombogenicity. Additionally, Figure 11B shows the relative velocity distribution of the fluid particles remaining in the cavities at the peak-ejection phase for the control and optimized designs, showing that the optimized design has a larger portion of the remaining fluid particles in the high velocity ranges (above 0.2 m s −1 ). In Figure 11C, we visualize the spatial distributions of the TAWSS for the control and optimized designs at the peakejection phase. Low TAWSS (red color) implies an increased risk of thrombus development in the region. At the end of the systolic phase, there was a larger area fraction of 42.4% of the endocardial surface exposed to TAWSS values below 0.2 Pa in the control design while only 17.4% of the area was exposed to low TAWSS in the optimized design. High thrombogenic risk regions indicated by the lowest values of TAWSS were located near the apex region of the endocardial surface of the control design, where blood is "trapped" by the actuator motion during the systolic phase. Increased blood flow velocities resulted in higher TAWSS values in the optimized design, which is consistent with the velocity distribution findings in Figure 11B.

Discussion and Conclusion
Here, we have described a design framework for optimizing the biomechanics and hemodynamics of a synthetic soft robotic myocardium. A finite element simulation was used to simulate the soft robotic composite consisting of multiple soft linearly contracting elements and a passive elastomeric matrix. We investigated the contribution of major design variables to the systolic capability of soft robotic myocardium, which is challenging experimentally, using the FE model as a virtual experimental platform. Such an approach may be applied generally for large-scale soft robotic devices containing multiple actuating units integrated to a passive soft matrix. Since we use a thermal expansion analogy to prescribe the deformation of actuating units, this method is not limited to modeling actuators driven by soft fluidic actuation but may also be applied to other types of soft contracting actuators powered by other energy sources such as thermal, electrical, and chemical.
Leveraging our understanding of the relevant design parameters for a soft robotic myocardium, we optimized the design using a combination of multiple actuator types. We demonstrated that the motion and hemodynamics produced by the optimized soft robotic myocardium are within the physiological range of a healthy LV during the systolic phase in a two-way FSI model. A relative assessment for thrombogenic risk potential  for our optimized design compared to the control design was developed based on blood velocity and TAWSS, which are associated with thrombogenicity. The optimized design showed superior fluid ejection due to high ΔLVV and had less lowvelocity particles remaining in the cavity at the end of the systolic phase. Additionally, it was demonstrated that the high thrombogenic risk regions near the apical surface indicated by low TAWSS in the control design containing only circumferential actuators was eliminated with the addition of helical actuators that break the axisymmetry, shorten the longitudinal aspect and prevent pooling/trapping of blood particles near the apex in the optimized design.
There are some limitations to this work. We do not fully capture the fluid dynamics of soft robotic myocardium and do not claim that our optimized design is free from and risk of thrombogenicity. Prior to implementation of our design for in vivo applications, a similar workflow of FSI may be used with specific solid structures representing the anatomy (such as heart valves and other internal cardiac structures) over a few cardiac cycles including the diastolic (filling) phase to fully capture the complex intraventricular flow patterns and to perform a comprehensive assessment of thrombogenic risk.
Overall, we showed that our computational workflow could be used to predict the behavior of a multi-unit soft robotic composite and optimize the desired behavior prior to prototyping and experimentation. Guided by a bioinspired design, we then optimized the design of soft robotic myocardium to achieve the desired functional requirements and cardiac biomechanics in a simplified variant of the design, which is more practical for implementation.

Experimental Section
Finite Element Model of Soft Robotic Myocardium: In our FE model, the soft robotic myocardium composed of passive soft matrix and active pneumatic artificial muscles was modeled. The matrix, which is composed of Ecoflex 00-30 (Smooth-On), was modeled as a Neo-Hookean hyperelastic material (C10 = 0.0113333, D1 = 1.96) using a 10-node modified quadratic tetrahedron (C3D10M) element. Soft contracting actuators that act as pneumatic artificial muscle units were modeled as linearly elastic materials (Figure 12A). A Young's modulus of 15.34 MPa for McKibben actuators and 0.5 MPa for flat pleated PAMs were used as inputs to describe their mechanical properties in the computational model. As previously described in Roche et al., [15] actuator geometric deformations during active contraction were replicated by implementing orthotropic expansion using thermal analogy ( Figure S1, Supporting Information).
To predict the post-actuation state, quasi-static analysis was performed using the Abaqus/Explicit (Dassault Systemes) solver while ensuring the kinetic energy to be less than 1% of the internal energy. Mass scaling was used to speed up the computation time while ensuring the percent mass change to be less than 10%. Mesh size was optimized by performing a mesh sensitivity test on the most complex model containing ≈100,000 elements ( Figure S2, Supporting Information). For each full-scale soft robotic myocardium, we assess its systolic performance under "no load" conditions to identify important design variables ( Figure 12B). The soft robotic myocardium was allowed to deform freely to observe the contribution from the matrix and actuator components to the change in the ventricular chamber volume and the degree of rotation. A linear ramp of actuation was applied over time, starting with an initial condition of no temperature and external pressure applied and the nodes on the basal plane remained fixed throughout the deformation. A closed artificial fluid cavity was defined between the inner ventricular wall and the artificial cap created at the basal plane from which the volume was calculated during contraction. To measure the degree of rotation of the soft robotic myocardium, a total of four markers were placed at the subapical plane on the epicardial surface and tracked over time. In separate sets of simulations, the outlet was blocked, and the fluid cavity model was used to sense the ΔLVP in the ventricle generated by the soft robotic myocardium. For the fluid cavity model, a bulk modulus of 2.4 MPa was used to achieve numerical stability.
Actuator Fabrication and Characterization: McKibben actuators were fabricated using a 6.35 mm-diameter PET braided mesh (TechFlex) and a 12 mm-diameter inner bladder that was made from thermally bonding two sheets of 38 µm-thick thermoplastic elastomer (Stretchlon 200; Airtech International). fPPAMs were fabricated as reported previously. [14] To characterize their active material stiffness, the ends of actuators were mounted onto the grippers of a mechanical tester (Instron 5944) and supplied with constant pressure for contraction. One end of grippers was lowered to allow full contraction at the zero-force state and then slowly pulled away at a rate of 60 mm min −1 to obtain the force-displacement curve from which the effective Young's moduli for each actuator were derived as previously described. [15] Fluid-Structure Interaction Modeling: A fully coupled FSI modeling approach was used to predict the hemodynamic response of the optimized design prototypes. The flow dynamics within the soft robotic myocardium were modeled in XFlow software (Dassault Systemes). While conventional numerical approaches use Navier-Stokes solver-based FE methods, they are still tied to major drawbacks such as imprecise meshing and highly empirical modeling of turbulence (RANS). [42,43] As the use of alternative formulations for the fluid flow equations becomes prevalent in the FSI applications, the Lattice Boltzmann method is becoming increasingly employed due to the benefits of the mesh-free method and the underlying fixed Cartesian grids. [44] Hence, this study relied on the Lattice Boltzmann formulation using the large-eddy simulation (LES) technique for turbulence modeling.
To perform a fully coupled FSI simulation, the fluid and structural domains were linked through the cardiac wall. In the structural domain, the soft robotic myocardium was initially at rest and a ramp of actuation was applied for 300 ms to simulate the cardiac contraction during the systolic phase. For the fluid domain, the cardiac cavity was pre-filled with particles in the initial state of the simulation (Figure 13). The cardiac wall was defined as the moving boundary condition to assess transient flow dynamics and cardiac wall deformation over the period. Time-dependent physiological aortic pressure waveform was assigned to the artificial cap created at the basal plane to represent an outflow boundary condition between LV and aorta. Unidirectional flow only was allowed to prevent backflow into the soft robotic myocardium during the systole. Blood was modeled as an incompressible Newtonian fluid with a density of 1050 kg m −3 and dynamic viscosity of 0.0035 Pa s). Due to its consistency and near-wall behavior, the wall-adapting local eddy model (WALE) is used. [45] Blood was modeled as an incompressible Newtonian fluid with a density of 1050 kg m −3 and dynamic viscosity of 0.0035 Pa s). For numerical stability, the time step Δt = 1 × 10 −5 s was applied to Adv. Funct. Mater. 2022, 32, 2206734 Figure 12. A) Modeling of the matrix-actuator composite using a thermal expansion analogy. B) FE prediction of the behavior of full-scale soft robotic myocardium at no load conditions. Fluid volume and pressure are measured at a reference point (RP) in an artificial fluid cavity and twisting motion is measured by tracking markers on the epicardial surface at the sub-apical plane. each simulation scenario, with a total simulation duration of 0.3 s representing the ventricular ejection. Grid independence was reached at 0.8 mm resolution with approximately 153,000 elements. To achieve an adequate wall boundary layer, the adaptive refinement method with 0.2 mm resolution was performed near the walls in the no-slip condition. The cardiac wall deformation and pressure/velocity characteristics were monitored for each full-scale soft robotic myocardium. Pressures were measured using virtual pressure sensors and the ejection volume was obtained by tracking the outflow at the outlet boundary and by integrating the flow rate over time.
Assessment of Thrombogenic Risk Potential: For relative assessment for thrombogenic risk, TAWSS was computed as a local hemodynamic factor as it is commonly associated with thrombus formation. [39,40,46,47] TAWSS hemodynamic index was computed using fluid simulations to determine regions with low blood flow velocities and endocardial surfaces exposed to low oscillatory shear flows. TAWSS magnitude can be determined based on the given equation where T is the integration period. TAWSS values would be lower in areas with low blood flow velocities (i.e., flow stagnation) precisely where the risk of thrombus formation is higher. [40] Experimental Prototype and Validation: A soft robotic myocardium prototype was built using a multistep casting process, previously described in Roche et al. [15] Briefly, multi-component molds were 3D printed for casting the matrix layers (inner, mid, and outermost) and actuators. Prior to casting, the braided meshes ("1/4" PET expandable sleeving, Techflex) were pre-coated with Ecoflex 00-30 (Smooth-On) over actuator cavity molds that were 3D printed in cylindrical and planar paths for the 3D and flattened McKibben actuators, respectively, to keep the airline access clear during silicone molding. Pre-fabricated inner bladders that were thermosealed using 38 µm thick thermoplastic elastomer sheets at 20 mm width were assembled within the siliconecoated mesh to form the McKibben actuators. The actuators were then aligned and secured in position with respect to the inner, middle, and outer layers of an idealized left ventricular geometry that were cast in the 3D-printed myocardial molds and were integrated embedded within an elastomeric matrix using uncured Ecoflex 00-30 as the glue to make a soft robotic myocardium prototype. Once the prototype was built, then the actuators were supplied with pressurized air from a custommade electropneumatic control system described in Roche et al. [10] and Park et al., [14] as the actuators used for the soft robotic myocardium were pneumatically actuated to generate desired anisotropic strain and contraction forces. The LV cavity was initially filled with water and then volume displacement was characterized as the soft robotic myocardium underwent actuation at incremental pressures (Figure 14).

Supporting Information
Supporting Information is available from the Wiley Online Library or from the author.