The ASTM D638 dogbone STL file was loaded and duplicated into n = 6 per shore group using GrabCAD (Stratasys). The shore levels used were 30, 35, 40, 50, 60, 70, 85, and 95. Each comprises different mixed ratios of the Stratasys rigid Vero material and flexible Agilus material. For blood vessel specimens, the same dogbone STL file was loaded and duplicated into n = 6 specimens per group into GrabCAD and assigned the six different blood vessel wall anatomy function groups: compliant, moderately compliant, slightly compliant, low compliant, semi-rigid, and rigid. Structural heart specimens consisted of myocardium, leaflet, chordae, and vessel wall with n = 4 specimens per group.

Printed dogbones were post-processed to remove all support structures. Each group (n = 4–6) was mechanically tested using an Instron 5,944 Single Column 2 kN universal Tensile Tester. Specimens were ramped in tension up to 10 mm displacement (1 mm for stiffer materials) over 60 s and then held at 10 mm displacement (1 mm for stiffer materials) for an additional 300 s with the resulting force, displacement, and time data exported to an Excel file. Once all specimens were printed and mechanically tested, they were fit using the constitutive models described in Nonlinear Viscoelastic Constitutive Model Theory using the methods described in Constitutive Model Parameter Fitting Using Inverse Finite Element Modeling.

The large deformation and nonlinear stress-strain relationship exhibited by the Stratasys J750 materials requires nonlinear constitutive models to account for this behavior. Adding rate-dependent and stress relaxation behavior additionally requires that we model J750 materials as viscoelastic. We sequentially describe the elastic and stress relaxation components of the resulting constitutive models.

A basic measure of deformation is the stretch ratio λ, a ratio of the specimen length after deformation l to the specimen before deformation l’. Furthermore, we have a stretch ratio in each of three directions, where we denote 1 as the x direction, 2 as the y direction, and 3 as the z direction: λ 1 = l 1 l ' 1 ;   λ 2 = l 2 l ' 2 ;   λ 3 = l 3 l ' 3 ;   (eq. 1) where λ_i is the stretch ratio in the ith testing direction. For deformations in the principal directions (i.e., no shear deformation) the deformation gradient tensor F_ij is defined in terms of the stretch ratios as: F i j = [ λ 1 0 0 0 λ 2 0 0 0 λ 3 ] (eq. 2)

The right Cauchy Deformation tensor C_ij is defined in terms of the deformation gradient tensor F_ij as: C i j = F k i F k j (eq. 3)

The scalar volume change J is the determinant of F_ij; J = det ( F i j ) (eq. 4)

For materials that are nearly incompressible (nearly zero volume change), the deviatoric portion of the deformation measures are defined as: F ˜ i j = J − 1 3 F i j       ⇒     λ ˜ i = J − 1 3 λ i (eq. 5a) C ˜ i j = F ˜ i j F ˜ k j       = J − 2 3 C i j (eq. 5b)

Strain energy functions W are used to characterize nonlinear elastic behavior. We utilized the simplest isotropic nonlinear elastic constitutive models that would characterize the 3D-printed Stratasys behavior, compressible and incompressible Neo-Hookean strain energy functions. These are defined in FEBio (FEBio User Manual 3.3) as:

Compressible Neo-Hookean model: W =   E 4 ( 1 + ν ) ( λ 1 2 + λ 2 2 + λ 3 2 − 3 ) − E 4 ( 1 + ν ) ln ⁡ J + ν E 2 ( 1 + ν ) ( 1 − 2 ν ) ( ln ⁡ J ) 2 (eq. 6a)

Incompressible Neo-Hookean model: W =   c 1 ( λ ˜ 1 2 + λ ˜ 2 2 + λ ˜ 3 2 − 3 ) + 1 2 K ( ln ⁡ J ) 2 = μ 2 ( λ ˜ 1 2 + λ ˜ 2 2 + λ ˜ 3 2 − 3 ) + 1 2 K ( ln ⁡ J ) 2 (eq. 6b) where E denotes Young’s modulus and ν the Poisson’s ratio for the compressible Neo-Hookean material (Eq. (5a)) with stretch ratios λ_i (Eq. 1) and Jacobian (volume change) J (Eq. (4)). For the incompressible Neo-Hookean material often written with the coefficient c₁ (Eq. (5b)) and(5μ) denotes a shear modulus (2*c₁) and K is a bulk modulus chosen as 500–1,000*c₁ to enforce the near incompressibility condition. Note that for the compressible material, the shear modulus, Young’s modulus, and Poisson’s ratio are related through the standard linear elastic relationship: μ = E 2 ( 1 + ν ) ⇒ c 1 = E 4 ( 1 + ν ) (eq. 7)

Thus, as the compressible material approaches incompressibility, Poisson’s ratio in the limit will reach 0.5 and Young’s modulus will be three times the shear modulus, six times c₁. The second Piola–Kirchoff (second PK) elastic stress tensor S i j e may be derived from the strain energy functions for compressible and incompressible materials by:

Compressible material: S i j e =   2 ∂ W ∂ C i j (eq. 8a)

Incompressible material: S i j e =   − p J C i j − 1 + J − 2 3 D e v [ 2 ∂ W ∂ C ˜ i j ] (eq. 8b) where Dev indicates the deviatoric component defined as: D e v ( S i j ) = S i j − 1 3 δ i j S k k with δ_ij being the Kronecker delta, and p is the hydrostatic pressure.

Accounting for the stress relaxation exhibited by the J750 materials requires a viscoelastic model. In this case, viscoelasticity was incorporated using a nonlinear stress-strain function to represent instantaneous elastic response coupled with stress relaxation through a reduced relaxation function that is independent of deformation. Thus, the nonlinear elastic kernels of Eq.(6a) and (6b) are modified by a reduced relaxation function that characterizes the relaxation behavior. Stress at a given time of loading is history-dependent, requiring integration of the strain history over the loading period together with a stress relaxation function that characterizes how stress relaxes overtime for a fixed displacement. We chose a 3-term Prony Series to model stress relaxation: G ( t − s ) = 1 + ∑ i = 1 N γ i e s − t τ i = 1 + γ 1 e s − t τ 1 + γ 2 e s − t τ 2 + γ 3 e s − t τ 3 (eq. 9) where γ₁, γ₂, and γ₃ are model coefficients, τ₁, τ₂, and τ₃ are model coefficients representing characteristic relaxation times, s is the current time in the loading history and t is the final loading time. The final loading time t for all tests was 360 s. The complete second Piola–Kirchoff (second PK) stress tensor as a function of time for a nonlinear compressible viscoelastic material can then be represented as: S i j ( t )   =   ∫ − ∞ t G ( t − s ) d S i j e d s d s (eq. 10) where S i j e is the elastic stress from Eq. (8a) and G ( t − s ) is the stress relaxation function from Eq. (9). For a nonlinear incompressible viscoelastic material, the time-dependent second PK stress tensor becomes: S i j ( t )   =   p J C i j − 1 + J − 2 3 ∫ − ∞ t G ( t − s ) d ( D e v [ S i j e ] ) d s d s (eq. 11) where p is the hydrostatic pressure, J is the volume change, C_ij is the right Cauchy Deformation tensor, S i j e is the elastic stress from Eq. (8b), and G ( t − s ) is the stress relaxation function from Eq. (9).

For both nonlinear viscoelastic models, we fit not only the nonlinear elastic coefficients (c₁ or E and ν) but also the stress relaxation coefficients γ₁, γ₂, γ₃, τ₁, τ₂, and τ₃ using inverse finite element modeling. This approach minimizes the difference in a least-squares sense over the entire 360 s loading period between the finite element model (Figure 1) reaction force at the rigid grips and the experimental force: M i n ( c 1   o r   E & v ) , γ 1 , γ 2 , γ 3 , τ 1 , τ 2 , τ 3     ∑ t = 1 # t i m e   s t e p s   f o r   360     sec ⁡ o n d s [ F e ( t ) − F m ( t ) ] 2 (eq. 12) where F_e denotes the experimental reaction force and F_m denotes the model reaction force. We fit both the compressible and incompressible versions of the Neo-Hookean elastic kernel, as many investigators model soft tissues as incompressible.

All finite element modeling and parameter optimization fitting of constitutive model coefficients were performed using the FEBio software suite, versions 2.91 and 3.1.0 (Maas et al., 2012; Maas et al., 2017). Due to the long computing time (∼10–100 h) for nonlinear viscoelastic material coefficient fitting (resulting from nonlinear iterations for the optimization fit and iteratively solving a nonlinear finite element problem within each optimization iteration), all models were run on the Hive Cluster high performance computing cluster supported by the Partnership for Advanced Computing Environment (PACE) at Georgia Tech. FEBio PreView version 1.5 was used to generate a finite element model of the ASTM D638 dumbbell test specimen using the same STL file used for the J750 printing. The dumbbell finite element model consisted of 14,515 10-node tetrahedral elements with 28,303 nodes (Figure 1A). The parameter optimization module in FEBio uses the Levenberg–Marquardt algorithm to update material parameters that minimize the square difference between model and experimental reaction forces. The testing grips in this model were considered rigid materials under the FEBio sliding elastic tension contact algorithm with the J750 material dumbbell specimen. Default convergence parameters were chosen in FEBio for both the nonlinear optimization iteration in the Levenberg-Marquardt algorithm and the sub-optimization nonlinear stress equilibrium problem, with the exception that the displacement convergence tolerance for the stress equilibrium problem was set to 0.01.

The nonlinear viscoelastic constitutive models applied here contain either eight (compressible Neo-Hookean case) or seven (incompressible Neo-Hookean case) parameters to fit a single uniaxial test. Fitting complex constitutive models using least-squares optimization approaches often leads to non-unique material parameters, i.e. multiple parameter sets will give equivalent fits to experimental data, also known as coefficient identifiability issues (Safa et al., 2021). Safa et al. (2021) recommend using a randomized multi-start least-squares fitting approach. While this process is feasible if the objective for fitting can be analytically derived, it is not feasible to use when each optimization solution uses an inverse nonlinear finite element modeling method requiring significant computing time. Therefore, to address coefficient variability, we averaged initial optimization results for all materials, labeled as Coefficient Set 1 (Coeff Set 1), and used these average values as the initial guess for a second optimization run with a lower and upper bound on the parameters set to ±10% of the initial guess. This second optimization run was labeled Coefficient Set 2 (Coeff Set 2).

A coefficient of determination (R ² value) was calculated in MATLAB (Mathworks) for each fit between reaction forces from the constitutive model fit used for the finite element model (Figure 1) and experimental reaction force data. An R ² value greater than 0.95 is considered a good fit for the constitutive model of experimental data (Humphrey, 2002).

A prime motivation for characterizing 3D-printed materials is to determine which 3D-printed material best mimics a given tissue behavior. This motivation plays into the most basic hypotheses concerning clinical applications of 3D-printed models. In training, we hypothesize that the 3D-printed material which best mimics the target tissue mechanics will provide the most realistic haptic experience, where the training outcome is assessed by survey methods. In medical device development, we hypothesize that device performance is best assessed using 3D-printed tissue models where the 3D-printed material most closely mimics tissue mechanics. For nonlinear viscoelastic materials, it is not a straightforward task to pick the materials that best match or at least bound a given tissue mechanical response from a wide range of materials (29 fit in the current study, but up to 96 in the current DAP database).

We, therefore, propose three metrics for assessing the 3D-printed material that best mimics either the nonlinear elastic or viscoelastic behavior of a given tissue. The first method is the total least square difference between forces from the 3D-printed material compared to the forces (stresses can also be used) from a tissue under the same time-dependent deformation. Note that this parameter is used to drive the initial optimization fitting of Eq. (12). The second method is using a two-sample Kolmogorov–Smirnov test (Press and Teukolsky, 1988). We used the KStest2 Kolmogorov–Smirnov statistics test version in MATLAB which tests the null hypothesis that two data vectors come from the same populations. In our case, we test the null hypothesis that force vectors matched over time for the stress relaxation test come from the same material, i.e., the 3D-printed material has the same force vs. time distribution as the tissue of interest. Note that this can readily be applied to force-displacement curves, stress-strain curves, or stress-time curves. Our hypothesis is that lower KStest2 p values will delineate 3D-printed materials that better match target tissue mechanical behavior. The third parameter is the difference in the total area under the force-time curve for the experimental data versus model fit. This area was calculated using the trapz trapezoidal numerical integration function in MATLAB.

The tissue-mimicking 3D-printed fit metrics may be implemented in two ways. First, if constitutive parameters for a given tissue are reported, which may include nonlinear elastic or nonlinear viscoelastic, isotropic, or anisotropic, these properties may directly be input to the FEBio model of Figure 1. The resulting force-time curve under the ramp displacement curve described in the mechanical testing section can be computed in FEBio and saved in an Excel file. The force-time curve for all 3D-printed Stratasys materials from the mechanical test is also stored in an Excel file. Both datasets are read into a MATLAB code that calculates the summed least-square difference, the KStest2 p value, and the difference in the area under the force-time curve. The program then ranks how well the 3D-printed materials represent tissue behavior based on the three fit parameters. If the material is anisotropic, this matching procedure may be performed over all desired testing directions, and the ranking of 3D-printed materials performed for the aggregate of testing direction responses. We illustrate this approach for the human tympanic membrane (Motallebzadeh et al., 2013) and the human nasal cartilage (Chang et al., 2020) for published nonlinear viscoelastic constitutive parameters, searching the J750 standard and the DAP blood vessel data. Nonlinear elastic tissue data can be readily assessed using this approach as well.

Second, if either force-displacement or stress-strain data is reported or can be derived from published data, the experimental setup from the published tissue setup may be modeled with the Stratasys material parameters to generate a database to compare with tissue results. Polzer et al. (2015) performed biaxial testing of the porcine aortic wall, providing extensive raw first Piola-Kirchoff stress vs. stretch data for axial and circumferential directions. The first PK stresses were converted to Cauchy stress vs. stretch results. We modeled the biaxial test results for one specimen in FEBio using rigid tension contact (Figure 2) and computed axial and circumferential Cauchy stress vs. stretch results for all compressible Neo-Hookean models of the Stratasys materials. We then computed the least-square difference, Kolmogorov-Smirnov statistics, and the area under stress-strain curve to assess which 3D-printed material best matched the experimental porcine aorta wall mechanical behavior.