# Nonlinear Viscoelastic Properties of 3D-Printed Tissue Mimicking Materials and Metrics to Determine the Best Printed Material Match to Tissue Mechanical Behavior Adam S. Verga, et al.

May 14, 2022

## Introduction

3D printing has become ubiquitous in clinical applications for a wide variety of disciplines. These clinical applications include 3D-printed implanted medical devices and tissue engineering scaffolds (Zopf et al., 2013; Di Prima et al., 2016; Ricles et al., 2018; Zhou et al., 2018) as well as models for surgical and procedure planning (Sommer et al., 2013; Pacione et al., 2016; Gocer et al., 2019; Chen et al., 2020a; Chen et al., 2020b; Bellia-Munzon et al., 2020; Capucha et al., 2020; Cho et al., 2020; Damon et al., 2020; Ghazi and Teplitz, 2020; Illmann et al., 2020; Javan et al., 2020; Kim et al., 2020; Levin et al., 2020; Liu et al., 2020; McMillan et al., 2020; Meglioli et al., 2020; Stramiello et al., 2020; Vannucci et al., 2020; Wang et al., 2020; Witowski et al., 2020; Ghazi et al., 2021). Although these applications may appear disparate at first glance, they are connected on two levels. First, tissue engineering scaffold development will benefit from simulating fixation and mechanical performance in tissue-mimicking materials 3D printed to simulate clinical defects. Second, both fields hypothesize that improving the match of the scaffold or model mechanical behavior to that of the tissue being replaced/simulated increases clinical benefit. For tissue engineering, this hypothesis states that a scaffold more closely matching the mechanics of the tissue to be replaced will produce better tissue regeneration (see for example Lin et al., 2004; Koh et al., 2019; Mardling et al., 2020; Chao et al., 2021; Ghorbani et al., 2021). For surgical planning, the hypothesis states that the closer a 3D-printed model matches tissue mechanics the better the clinical and training outcome. This is especially critical for example in haptic feedback for training (Grillo et al., 2018; Ratinam et al., 2019; Pietrabissa et al., 2020; Tejo-Otero et al., 2020), for practicing complex surgical procedures like transseptal puncture in cardiology (Bezek et al., 2020), and in developing cardiac valves tested in flow loops to achieve realistic hemodynamic and fluid–structure interaction (Vukicevic et al., 2017; Ferrari et al., 2019; Ferrari et al., 2020; Levin et al., 2020; Vukicevic et al., 2020).

Rigorously testing these hypotheses requires two fundamentals. First, the properties of the tissues of interest and the 3D-printed scaffolds and models must be characterized using the most appropriate and comparable constitutive models. Especially for the soft tissues, this would include nonlinear elastic and viscoelastic models. Second, there must be appropriate metrics that can be used to objectively compare constitutive properties for many different materials, even when these properties are very different. The goal of this study is to address both issues by 1) fitting nonlinear viscoelastic models to 3D-printed materials using an inverse finite element optimization approach by FEBio(Maas et al., 2012; Maas et al., 2017) and 2) developing metrics to compare the mechanical behavior tissue and 3D-printed materials using both nonlinear elastic and nonlinear viscoelastic constitutive models. We specifically focus in this study on using 3D-printed Polyjet materials.

Matching complex and varied tissue mechanical behavior has motivated the use of multi-material Polyjet 3D printers like the Stratasys Connex line and Stratasys J750/J850 Polyjet printers to create tissue-mimicking phantoms. Since detailed information on constitutive models for these materials is not readily available, recent studies have sought to characterize Polyjet material mechanical properties in relation to tissue properties. Cloonan et al. (2014) characterized the nonlinear elastic behavior of Connex 3D-printed Tango material using three to five term Ogden strain energy functions. They showed qualitative agreement with the nonlinear elastic Abdominal Aorta Aneurysm (AAA) tissue stress distributions (Raghavan and Vorp, 2000). Ruiz de Galarreta et al. (2017) used an exponential anisotropic strain energy function to compare Stratasys TangoPlus material stress distributions to those generated using the same strain energy function developed by Choi and Vito for canine pericardium (Choi and Vito, 1990). Finally, other studies have used simpler constitutive models like linear elastic modulus or Shore hardness values (Severseike et al., 2019; Bezek et al., 2020; Tejo-Otero et al., 2020), to relate Polyjet material to tissue mechanics.

While an important step, these studies do not provide the extensive constitutive modeling needed to completely characterize Polyjet materials in relation to tissue behavior for four reasons. First, Polyjet materials exhibit nonlinear viscoelastic behavior with stress relaxation, large deformation, and nonlinear stress-strain behavior. Second, tissues exhibit nonlinear anisotropic viscoelastic behavior with an extremely broad range of constitutive parameters, see for example (Puso and Weiss, 1998; Holzapfel et al., 2004; Motallebzadeh et al., 2013; Polzer et al., 2015; Chang et al., 2020). Since both J750 materials and tissues exhibit nonlinear viscoelastic behavior, it is necessary to characterize the Polyjet materials using nonlinear viscoelastic constitutive models to provide the most complete comparison of 3D-printed material to tissues. The nonlinear viscoelastic models can be used to compare with the tissue nonlinear elastic constitutive models as well. Finally, given that the Stratasys J750 with Digital Anatomy Printing (DAP) contains 96 distinct materials, a methodology to automatically determine the best-fit 3D-printed material from a large material database of nonlinear viscoelastic materials to a given tissue is needed. For these purposes, a complete constitutive characterization is needed for the printed materials to generate the needed force-displacement, force-time, stress-strain, and/or stress-time curves that can be used to compute metrics comparing 3D-printed behavior to tissue behavior.

In this study we present compressible and incompressible isotropic nonlinear viscoelastic constitutive models fit for the Stratasys J750 fundamental Shore, the DAP blood vessel materials, and the DAP structural heart materials. These are 29 of the 96 (at present) materials available in this system. We also present three metrics (least square difference, Kolmogorov–Smirnov statistics, and the area under force/stress-time/displacement/strain curves) computed in a MATLAB (MathWorks) program from a database of 3D-printed materials to find the best-matched material to mimic a given tissue mechanical behavior. Examples are shown for the tympanic membrane, the nasal cartilage, and the aortic wall tissue.

## Materials and Methods

Test specimen geometry for both the 3D-printed specimens and finite element models for optimization was based on the ASTM D638 “Standard Test Method for Tensile Properties of Plastics” Type V specimen geometry (length: 63.5 mm, maximum width: 11 mm, thickness: 2 mm). Shore 30/35/40/50/60/70/85/95, DAP Blood Vessel materials, and DAP Structural Heart Materials (29 materials total) were printed and tested with n = 4-6 specimens per group.

### 3D Material Printing

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.

### Specimen Mechanical Testing

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.

### Nonlinear Viscoelastic Constitutive Model Theory

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=l1l‘1; λ2=l2l‘2; λ3=l3l‘3;$

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 Fij is defined in terms of the stretch ratios as:

$Fij=[λ1000λ2000λ3]$

The right Cauchy Deformation tensor Cij is defined in terms of the deformation gradient tensor Fij as:

The scalar volume change J is the determinant of Fij;

For materials that are nearly incompressible (nearly zero volume change), the deviatoric portion of the deformation measures are defined as:

$F˜ij=J−13Fij ⇒ λ˜i=J−13λi$

$C˜ij=F˜ijF˜kj =J−23Cij$

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= E4(1+ν)(λ12+λ22+λ32−3)−E4(1+ν)ln⁡J+νE2(1+ν)(1−2ν)(ln⁡J)2$

Incompressible Neo-Hookean model:

$W= c1(λ˜12+λ˜22+λ˜32−3)+12K(ln⁡J)2=μ2(λ˜12+λ˜22+λ˜32−3)+12K(ln⁡J)2$

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 c1 (Eq. (5b)) and(5μ) denotes a shear modulus (2*c1) and K is a bulk modulus chosen as 500–1,000*c1 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:

$μ=E2(1+ν)⇒c1=E4(1+ν)$

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 c1. The second Piola–Kirchoff (second PK) elastic stress tensor

$Sije$

may be derived from the strain energy functions for compressible and incompressible materials by:

Compressible material:

$Sije= 2∂W∂Cij$

Incompressible material:

$Sije= −pJCij−1+J−23Dev[2∂W∂C˜ij]$

where Dev indicates the deviatoric component defined as:

$Dev(Sij)=Sij−13δijSkk$

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=1Nγies−tτi=1+γ1es−tτ1+γ2es−tτ2+γ3es−tτ3$

where γ1, γ2, and γ3 are model coefficients, τ1, τ2, and τ3 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:

$Sij(t) = ∫−∞tG(t−s)dSijedsds$

where

$Sije$

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:

$Sij(t) = pJCij−1+J−23∫−∞tG(t−s)d(Dev[Sije])dsds$

where p is the hydrostatic pressure, J is the volume change, Cij is the right Cauchy Deformation tensor,

$Sije$

is the elastic stress from Eq. (8b), and

$G(t−s)$

is the stress relaxation function from Eq. (9).

### Constitutive Model Parameter Fitting Using Inverse Finite Element Modeling

For both nonlinear viscoelastic models, we fit not only the nonlinear elastic coefficients (c1 or E and ν) but also the stress relaxation coefficients γ1, γ2, γ3, τ1, τ2, and τ3 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:

$Min(c1 or E&v),γ1,γ2,γ3,τ1,τ2,τ3 ∑t=1#time steps for 360 sec⁡onds[Fe(t)−Fm(t)]2$

where Fe denotes the experimental reaction force and Fm 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.

FIGURE 1. (A) Finite element model of ASTM D638 dumbbell test specimen used to test J750 materials. (B) Actual 3D printed J750 Shore material undergoing tensile stress relaxation test in Instron. The finite element model was created using FEBio Preview 1.5. The model was analyzed using either FEBio 2.9.1 or FEBio 3.3.1 using the parameter optimization tool to fit the nonlinear viscoelastic model parameters. Rigid bodies at the end of the specimen allow for tension only contact mimic specimen gripping. The right platen was displaced to match the time-displacement curve from experimental tests.

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 (R2 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 R2 value greater than 0.95 is considered a good fit for the constitutive model of experimental data (Humphrey, 2002).

### Metrics for Best Matched 3D-Printed Material to Tissue Mechanics

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.

FIGURE 2. Finite element model was used to simulate biaxial testing of the porcine aorta wall tissue specimens reported by Polzer et al. Specimen dimensions matched those given by Polzer and displacements of rigid test grips were set to match the stretch ratios and loading rate. Cauchy stress versus stretch was calculated for all J750 constitutive models using the model. Model of aorta wall specimen contained 5,184 8-node hexahedral elements and 6,845 nodes. Sliding elastic tension contact was assumed between rigid grips and specimen.

## Results

Complete Coeff Set 1 and Coeff Set 2 parameter fits and corresponding R2 values are presented for the following using Neo-Hookean Prony Series Viscoelastic models under a compressible or incompressible assumption as noted:

2) J750 Shore materials; incompressible (Table 2)

3) J750 DAP Blood Vessel Wall materials; compressible (Table 3)

4) J750 DAP Blood Vessel Wall Materials; incompressible (Table 4)

5) J750 DAP Structural Heart materials; compressible (Table 5)

6) J750 DAP Structural Heart materials; incompressible (Table 6)

TABLE 1. Nonlinear viscoelastic compressible Neo-Hookean fit to J750 Shore materials. Table of compressible nonlinear viscoelastic compressible Neo-Hookean parameters for the J750 Shore materials along with the R2 values of the fit between the constitutive model and experimental data (n = 6 specimens per material except for Shore 30 with n = 5 specimens per material). Coefficient Set refers to the initial guess (first) with wide bounds on the parameters followed by a second optimization run with the bounds being the mean parameters from the first optimization run ±10% of the mean.

TABLE 2. Nonlinear viscoelastic incompressible Neo-Hookean fit to J750 Shore materials. Table of compressible nonlinear viscoelastic compressible Neo-Hookean parameters for the J750 Shore materials along with the R2 values of the fit between the constitutive model and experimental data (n = 6 specimens per material except for Shore 30 with n = 5 specimens per material). Coefficient Set refers to the initial guess (first) with wide bounds on the parameters followed by a second optimization run with the bounds being the mean parameters from the first optimization run ±10% of the mean.

TABLE 3. Nonlinear viscoelastic compressible Neo-Hookean fit to J750 Digital Anatomy Printer (DAP) Blood Vessel Materials. Table of compressible nonlinear viscoelastic compressible Neo-Hookean parameters for the DAP, Blood Vessel materials along with the R2 values of the fit between the constitutive model and experimental data (n = 6 specimens per material except for Shore 30 with n = 5 specimens per material). Coefficient Set refers to the initial guess (first) with wide bounds on the parameters followed by a second optimization run with the bounds being the mean parameters from the first optimization run ±10% of the mean.

TABLE 4. Nonlinear viscoelastic incompressible Neo-Hookean fit to J750 Digital Anatomy Printer (DAP) Blood Vessel Materials. Table of compressible nonlinear viscoelastic compressible Neo-Hookean parameters for the DAP, Blood Vessel materials along with the R2 values of the fit between the constitutive model and experimental data (n = 6 specimens per material except for Shore 30 with n = 5 specimens per material). Coefficient Set refers to the initial guess (first) with wide bounds on the parameters followed by a second optimization run with the bounds being the mean parameters from the first optimization run ±10% of the mean.

TABLE 5. Nonlinear viscoelastic compressible Neo-Hookean fit to J750 Digital Anatomy Printer (DAP) Structural Heart Materials. Table of compressible nonlinear viscoelastic compressible Neo-Hookean parameters for the Structural Heart materials along with the R2 values of the fit between the constitutive model and experimental data (n = 4 specimens per). Only the second optimization run with the bounds being the mean parameters from the first optimization run ±10% of the mean is presented.

TABLE 6. Nonlinear viscoelastic incompressible Neo-Hookean fit to J750 Digital Anatomy Printer (DAP) Structural Heart Materials. Table of compressible nonlinear viscoelastic compressible Neo-Hookean parameters for the Structural Heart materials along with the R2 values of the fit between the constitutive model and experimental data (n = 4 specimens per material). Only the second optimization run with the bounds being the mean parameters from the first optimization run ±10% of the mean is presented.

All models fit the J750 Standard and DAP material behavior well with R2 > 0.95 (in most cases R2 > 0.99) except for the incompressible Neo-Hookean Coeff Set 2 for the Shore 30 material, with R2 = 0.922. As expected, Coeff Set 1 starting with a large parameter space fit the experimental data well, but exhibited a wide variation in coefficients with standard deviations up to 100% or greater of the mean values, demonstrating non-uniqueness, or equivalently coefficient identifiability issues. Coeff Set 2 using the average of Coeff Set 1 as the starting guess had much lower resulting variations due to the 10% fitting bound. Despite the narrower bound, R2 fits were equivalent for Coeff Set 2 as Coeff Set 1, except for the Shore 30 material. In this case, the Coeff Set 2 parameters tended towards the lower bound of the parameter space. Examples of fits for both compressible and incompressible nonlinear viscoelastic models are shown in Figure 3.

FIGURE 3. Example of nonlinear viscoelastic constitutive model fit and experimental uniaxial test data for (A) compressible model and Shore 30 material specimen, (B) incompressible model and Shore 30 material specimen, (C) compressible model and Shore 85 material specimen, and (D) incompressible model and Shore 85 material specimen.

Rankings of J750 standard Shore materials and Blood Vessel materials are shown for compressible nonlinear elastic modulus (Eq.(6a); see Figure 4), incompressible nonlinear elastic shear modulus (Eq. (6b); see Figure 5), and percent force relaxation (Figure 6). Note that many J750 materials exhibit significant stress relaxation, between 8 and 50% of peak force (Figure 6). As expected, the compressible nonlinear elastic modulus was nearly six times the incompressible c1 coefficient for all materials: 5.96 ± 0.46.

FIGURE 4. Ranking of J750 3D printed Shore, the DAP blood vessel materials, and the DAP structural heart materials for the elastic modulus of the compressible nonlinear strain energy function of Eq. 6a.

FIGURE 5. Ranking of J750 3D printed Shore, the DAP blood vessel materials, and the DAP structural heart materials for the c1 component of the incompressible Neo-Hookean strain energy function of Eq. 6b.

FIGURE 6. Ranking of J750 3D printed Shore, the DAP bloodvessel materials, and the DAP structural heart materials for the degree of force relaxation over time.

The use of least square differences, KStest2, and difference in total area under the force-time curves revealed similar trends in identifying J750 materials that best matched the nonlinear viscoelastic behavior of native tissues. Figure 5 shows the three best J750 material fit to the nonlinear viscoelastic behavior of the human tympanic membrane (Figure 7A) and the human nasal cartilage (Figure 7B) from published nonlinear viscoelastic constitutive data (Motallebzadeh et al., 2013; Chang et al., 2020) using the dumbbell model (Figure 1). In the case of the tympanic membrane, the metrics provide the best 3D-printed material choice that brackets the tympanic membrane mechanical response.

FIGURE 7. Calculation of least-square difference (LS), Kolmorogorov–Smirnov statistic (KS), and differences in area under the force-time curve between three closest J750 materials and (A) human tympanic membrane nonlinear viscoelastic properties (Motallebzadeh et al., 2013) and (B) human nasal cartilage nonlinear viscoelastic properties (Chang et al., 2020). Note that lower magnitude values of each parameter indicate a J750 parameter that more closely matches native tissue behavior. Results suggest that either the Shore 70 or Rigid Blood Vessel material best match the tympanic membrane (A), while the Shore 30 material best matches human nasal cartilage (B).

The material matching parameters were also calculated for one porcine aortic wall specimen under biaxial tension (Figure 2) (Polzer et al., 2015). Note that this was a ramp test to assess elastic response, not a stress relaxation test to assess viscoelastic response. However, the full viscoelastic 3D-printed material model was used and subject to the loading rate specified by Polzer et al. Parameters were calculated for Cauchy stress versus stretch results in both the axial and circumferential test directions, although plot results are only shown for the axial direction (Figure 8). In this case, the most compliant material tested (Shore 30) was chosen but was stiffer than the aorta wall.

FIGURE 8. Calculation of tissue matching parameters for stress-strain elastic curve of porcine aortic wall specimen (Polzer et al., 2015). Again, the lowest magnitude values for all parameters agree, indicating that the Shore 30 J750 material is closest matched to this aortic wall specimen.

## Discussion

We demonstrated that nonlinear three-term Prony series viscoelastic constitutive models, with either a compressible or incompressible Neo-Hookean strain energy kernel, fit the nonlinear stress relaxation behavior of Stratasys J750 materials well, with R2 values greater than 0.95 for 42 of 43 material fits (2 fits per material for J750 and DAP Blood Vessel and 1 for DAP Structural Heart). The remaining fit had an R2 value of 0.922. We also demonstrated that a two-step optimization approach, initially performing a fit with a broad parameter range, then averaging these parameters to use as the initial starting guess with a ±10% bound for a second optimization run, produced consistent, narrowly bound constitutive parameters that fit experimental data equally well. Although it cannot be proven that the second round fit coefficients are indeed unique, the results are likely to provide more rigorous results for comparing to native tissue properties or performing finite element simulations and subsequent experimental testing for procedure training or medical device development.

We presented new metrics to match 3D-printed material behavior against tissue properties by calculating least-squares differences, Kolmogorov–Smirnov statistics, and differences in curve areas that can be used for force-time curves, force-displacement curves, stress-time curves, or stress-strain data. All these metrics were consistently lower for materials that better-matched tissue mechanical behavior, and thus provide an automated way to choose the best material from a 3D printing database for matching a given tissue behavior. A weighted average of the three parameters may give the best prediction. Such an automated approach becomes more critical when the number of materials becomes exceedingly large, as the DAP catalog for the J750 Polyjet printer alone contains 96 different combinations.

Exactly matching tissue nonlinear viscoelastic properties remains a daunting challenge for a number of reasons. First, tissue properties are extremely varied and complex, exhibiting anisotropy due to embedded fibers with complicated orientation distributions (Holzapfel et al., 2004; Sommer et al., 2013; Limbert, 2017) which 3D-printed materials do not replicate. Tissues exhibit a classic strain stiffening response, with soft materials stiffening at higher strains due to these embedded fibers. In contrast, most synthetic biomaterial polymers for modeling or implantation (Mitsak et al., 2012; Ramaraju et al., 2020), including the J750 3D-printed materials, exhibit Neo-Hookean responses. These differences will clearly create a challenge in matching 3D-printed materials to tissue response. An interesting approach to address this mismatch has been put forward by Chen et al. (2018), who proposed printing stiffer sinusoidal fibers in softer matrices on Polyjet systems to mimic the embedded fibers of tissues.

A second challenge is the extremely varied models in the literature used to characterize tissue nonlinear elastic and viscoelastic properties which makes an “apples to apples” comparison to commonly used nonlinear viscoelastic constitutive models difficult. This widely varied data suggest the need for a standardized characterization of nonlinear elastic and viscoelastic constitutive models for 3D-printed materials, natural tissues, and engineered tissues to create a common database. In other words, common testing and constitutive modeling approaches including common nonlinear elastic strain energy and stress relaxation functions (like the Prony series) need to be applied to all these materials. Finally, although Polyjet and mulitjet printers like the Stratasys J750 have significantly advanced our capability to print and design a wide range of nonlinear material behaviors, printing extremely compliant materials to match tissues like esophagus and blood vessels remains a challenge.

## Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

## Author Contributions

AV: 3D-printed samples, performed mechanical testing, prepared finite element files for the optimization, and wrote portions of the manuscript. ST: 3D-printed samples, performed mechanical testing, and wrote parts of the manuscript. YG: printed samples, performed mechanical testing, and prepared finite element files for the optimization. AP: 3D-printed samples and performed mechanical testing. SH: ran FEBio analyses, fit and interpreted data, wrote portions of the manuscript, and designed the study.

## Funding

This study was supported by the Children’s Healthcare of Atlanta Pediatric Trust Fund.

## Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

## Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

## Acknowledgments

The cluster computing work used the Hive cluster, which is supported by the National Science Foundation under grant number 1828187. Cluster computing was also supported in part through research cyberinfrastructure resources and services provided by the Partnership for an Advanced Computing Environment (PACE) at the Georgia Institute of Technology, Atlanta, Georgia, USA.