Erratum: “A Numerical Scheme for Anisotropic Reactive Nonlinear Viscoelasticity” [ASME J. Biomech Eng., 2023, 145(1), p. 011004; DOI: 10.1115/1.4054983]

In our recently published study [1], we reported on a numerical scheme for the framework of reactive viscoelasticity initially presented in our prior theoretical report [2]. We also extended that original theoretical framework to incorporate multiple families of weak bonds, which break and reform into a stress-free state under loading, to be optionally recruited with increasing amount of loading (Sec. 2.3 of the paper [1]). We explained that this weak-bond recruitment could be easily incorporated into the existing framework via Eq. (2.31) of the paper. We provided a pseudo-code that incorporated Eq. (2.31) in Algorithm 2 of the paper. Finally, we illustrated the application of weak-bond recruitment in the curve-fitting results of Fig. 9 of the paper, which analyzed the response of the tissue to consecutive ramp-and-hold displacement profiles.

Upon release of the code in the open-source finite element software FEBio [3], we got feedback from a user (see Acknowledgments below) that this algorithm produced nonphysical responses when a tissue is subjected to uni-axial prescribed cyclical displacement with constant amplitude, such as load amplification along the loading direction, instead of the expected load attenuation. Here, we present a correction to Eq. (2.31), and an update to Algorithm 2 and Fig. 9, which overcome this original limitation.

where Ξ is the strain measure described in Eq. (2.29) of the paper, and $F(Ξ)$ represents the enhancement in weak bond mass fraction relative to the zero strain state, as illustrated in the example constitutive model of Eq. (2.30). The main difference between this corrected model and the original model is that the new model accounts for the enhancement $F(Ξ(t))$ in the weak bond mass fraction based on the strain measure Ξ evaluated at the current time t. The original formulation employed the maximum strain measure achieved in the strain history of the material, up until the first time-step t^v when weak bonds from generation u started breaking.

This formulation adjusts the availability of weak bonds as the strain varies; in particular, it accommodates cyclical strains.

The results corresponding to Fig. 9 of the original paper become updated when employing this new formulation as follows.

We refitted the experimental data to this model to extract the material parameters μ^e and μ^a (shear moduli of the strong and weak bonds); τ₁₀, τ₂₀, τ₂₁, and s₂ for the reduced relaxation function; and μ and $Ξ0$ for the recruitment function. This eight-parameter fit produced $μe=0.40 kPa$⁠, $μa=1.79 kPa, τ10=0.244 s,τ20=205 s$⁠, $τ21=2575 s, s2=0.0096, μ=30.0$⁠, and $Ξ0=0.971$⁠. The nonlinear regression coefficient between the fit and data was $R2=0.991$⁠, and the fit is shown in Fig. 9(a). To better understand the strain-dependence of $τ2(K2)$ and $F(K2)$ for this fit, these functions are shown in Fig. 9(b), calculated using τ₂₀, τ₂₁, and s₂ for $τ2(K2)$ and μ, $Ξ0$ and α for $F(K2)$⁠.

The most significant difference we found in the curve-fit of the experimental data of Fig. 9, when compared to the original publication, was the observation of a smoother decrease with increasing strain, in the relaxation time τ₂ of the Malkin-distortion reduced relaxation function (Fig. 9(b)).

The authors would like to thank Dr. Lance Frazer of the Southwest Research Institute, San Antonio, Texas, for identifying the nonphysical response when modeling cyclical loading of a reactive viscoelastic material with weak-bond recruitment.

• National Institute of General Medical Sciences (Award ID: R01 GM083925; Funder ID: 10.13039/100000057).

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health.