Modelling the anisotropic inelastic response of polymeric scaffolds for in situ tissue engineering applications

1. Introduction

Cardiovascular diseases, which affect the heart and blood vessels, remain the leading cause of death worldwide. In particular, the aortic valve, which connects the aorta with the left ventricle and prevents the backflow of blood into the heart during diastole, may be affected by valve stenosis and regurgitation [1]. Today, these pathological conditions are often treated with a minimally invasive surgical procedure called transcatheter aortic valve implantation, in which a replacement heart valve is positioned through a catheter [2]. The most common solutions include mechanical and bioprosthetic heart valves manufactured from processed animal tissue. However, both categories have several drawbacks, including poor biocompatibility and age-related structural deterioration [3]. Innovative solutions are made available through in situ tissue engineering, also known as endogenous tissue restoration (ETR), in which a temporary bioabsorbable polymeric scaffold provides the biomechanical structural characteristics until the native tissue is produced [4–8]. In addition to heart valves, ETR is a promising technique for replacement vascular grafts [9,10]. The key benefits of such an approach are that it avoids the long-term presence of foreign materials in the body, provides a replacement heart valve with a native-like response to haemodynamic loads and has the ability to remodel and regenerate [4].

Cardiovascular tissues appear as sophisticated assemblies of fibres and cells that respond to a range of biomechanical and biochemical stimuli and can adapt to physiological and pathological conditions [11]. Therefore, the design of functional heart valves is a complex task that requires knowledge of the structural and mechanical characteristic of the native tissue as well as the biomechanical loads in vivo [12]. Furthermore, the success of ETR depends on the cellular response to the scaffold, e.g. with regard to infiltration, adhesion and immunogenic response [8]. In terms of mechanical behaviour, heart valves require excellent fatigue resistance to continuous exposure from cyclic pressures, shear stresses and strains resulting from haemodynamic loads [4]. In this context, the amount and spatial arrangement of collagen fibres are of fundamental importance, and it has been observed that the formation and organization of restored native tissue can be controlled by incorporating anisotropy into the scaffold [13]. In addition, fibre diameter and alignment correlate with pore size and thus have a direct influence on cell infiltration [14]. Synthetic polymers produced by electrospinning are among the most promising materials for biodegradable scaffolds with non-woven meshes featuring micrometre-sized fibres [15,16].

Given the large number of variables involved in tissue engineering design, simulation-based methods are highly attractive [17–19]. Computational tools allow the prediction of the entire ETR process. However, the accuracy of their results is highly dependent on constitutive models used to describe the short-term mechanical response of the biodegradable scaffold, its interaction with the haemodynamic environment, the process of cell infiltration and growth of native tissue, and finally the scaffold degradation [20–23]. The mechanical behaviour of the scaffold is similar to other fibre-reinforced soft composites, including biological tissues [24], textiles [25], fibre-reinforced polymers [26] and biopolymer gels [27]. These materials are endowed with peculiar features that are a direct expression of the arrangement of their microstructure and the interactions within it. More specifically, the mechanical response of electrospun polymeric meshes is characterized by large deformations, nonlinearity and anisotropy [14,28,29]. In addition, numerous inelastic effects are also observed in the experiments, including stress softening, permanent deformation and rate dependency [30]. In fact, these features are also found in vascular tissues. For example, rate dependence is particularly relevant for the aortic valve, which in vivo is subject to high deformation rates [31].

Here, we present a computational model to describe the short-term inelastic response of polymeric scaffolds for tissue engineering applications. The constitutive framework is based on a recently formulated continuum mechanics approach that combines anisotropy, stress softening, permanent deformations, and rate dependency [32]. This base model served as the foundation for modelling the ETR process in a biodegradable scaffold embedded in a thorough verification, validation and uncertainty quantification plan [33,34], as developed in the Horizon 2020 Project of the European Union SimInSitu. Aiming to simulate the entire ETR process, fluid–structure interaction, growth and remodelling, and patient-specific geometry, an efficient numerical implementation in industry-standard finite element (FE) codes is a key goal of this research. Hence, we derive simple algorithmic expressions for the computation of the spatial stress tensor and the consistent tangent modulus needed for standard implementations in implicit solvers developed in the updated-Lagrangian formulation. The FE code verification is performed by comparing the implementation in user-material subroutines for the commercial Abaqus software [35] with analytical solutions for elementary deformations. We also consider a large-scale representative example consisting of the simulation of a parallel-plate compression test with cyclic loading. In particular, we describe the experimental procedure for the mechanical characterization of the electrospun polymeric scaffold and use FE simulation to validate the proposed computational model.

2. Computational model

3. Representative examples

To show the performance of the numerical model, some selected representative examples are presented in this section. For this purpose, the numerical model was implemented as a user-defined material for the static implicit solver of the commercial FE software Abaqus Standard [35]. All examples assume that the material is perfectly incompressible. Mixed displacement–stress element formulations, with linear interpolation for displacements and constant pressure, are employed.

3.1. Homogeneous deformations

In order to verify the numerical model and to illustrate the physical behaviour described by the constitutive framework, we have selected three boundary-value problems. In the case of perfect incompressibility, these can be solved explicitly, since the deformation gradient is completely prescribed, and the unknown Lagrange multiplier p in (2.11) results from the solution of a system of equations. In the following, the stress–strain curves obtained from a single FE are compared with the analytical solution. We have considered representative material parameters of an electrospun polymeric material (table 1). Due to confidentiality agreements with the manufacturing company of the polymeric scaffold, Xeltis BV (Eindhoven, The Netherlands), the material parameters and the results shown in this section have been normalized.

Table 1. Summary of material parameters for an electrospun polymer for heart valve scaffolds. Stiffness-like quantities are normalized with μ_‖ and relaxation times with τ₁.

Collapse

material parameter	value
structure tensor eigenvalue H11 (–)	0.55
structure tensor eigenvalue H22 (–)	0.27
coaxial stiffness μ‖ (normalized)	1.00
coaxial first shape parameter γ‖ (–)	1.45
coaxial second shape parameter δ‖ (–)	0.01
perpendicular stiffness μ⊥ (normalized)	0.21
perpendicular first shape parameter γ⊥ (–)	2.49
perpendicular second shape parameter δ⊥ (–)	0.01
Mullins maximum damage parameter ηm0 (–)	0.78
Mullins damage evolution parameter αm (–)	27.26
maximum residual stress parameter ηr0 (–)	0.19
residual stress evolution parameter αr (–)	2.38
viscoelastic first relaxation time τ1 (normalized)	1.00
viscoelastic first inverse stiffness ${\mu }_1^{{\hspace{0.17em}}^{\prime }-1}$ (normalized)	1.54
viscoelastic second relaxation time τ2 (normalized)	15.98
viscoelastic second inverse stiffness ${\mu }_2^{{\hspace{0.17em}}^{\prime }-1}$ (normalized)	1.31
viscoelastic third relaxation time τ3 (normalized)	145.01
viscoelastic third inverse stiffness ${\mu }_3^{{\hspace{0.17em}}^{\prime }-1}$ (normalized)	0.44

The results shown below apply to deformations in the axial–circumferential plane (E₁, E₂). The full set of deformations is reported in the electronic supplementary material, figures S2–S4, in addition to another benchmark example to verify the convergence behaviour (electronic supplementary material, figures S5 and S6).

3.1.1. Biaxial extension

A biaxial deformation for an incompressible material is prescribed by the following deformation gradient:

$$\mathbf{F}={\lambda }_1{\mathit{e}}_1\otimes {\mathit{E}}_1+{\lambda }_2{\mathit{e}}_2\otimes {\mathit{E}}_2+{({\lambda }_1{\lambda }_2)}^{-1}{\mathit{e}}_3\otimes {\mathit{E}}_3,$$3.1

where λ_i are the principal stretches, and we have assumed that ${\mathit{E}}_1$ and ${\mathit{E}}_2$ are the principal material directions aligned with the axes of biaxial extension. Using the plane stress condition (2.12), the Cauchy stress tensor is given by

$$\mathit{\sigma }({\lambda }_1,{\lambda }_2)=[\mathbb{I}-\mathbf{I}\otimes ({\mathit{e}}_3\otimes {\mathit{e}}_3)]\hspace{0.17em}:\hspace{0.17em}{\mathit{\sigma }}_{\mathrm{iso}}({\lambda }_1,{\lambda }_2),$$3.2

where ${\mathit{\sigma }}_{\mathrm{iso}}=\mathbf{F}{\mathbf{S}}_{\mathrm{iso}}{\mathbf{F}}^{\mathrm{T}}$.

The results are illustrated in figure 1a for a cyclic loading protocol with two different strain rates and equibiaxial tension λ₁ = λ₂ = λ.

• Download figure
• Open in new tab
• Download PowerPoint

3.1.2. Uniaxial extension

A uniaxial deformation for an incompressible material is defined by the deformation gradient in (3.1). In addition, we prescribe λ₁ = λ for a uniaxial extension along the principal direction ${\mathit{E}}_1$. The Cauchy stress tensor for an anisotropic material is obtained from the solution of

$$\begin{array}{rl}\mathit{\sigma }(\lambda ,{\lambda }_2)& =[\mathbb{I}-\mathbf{I}\otimes ({\mathit{e}}_3\otimes {\mathit{e}}_3)]\hspace{0.17em}:\hspace{0.17em}{\mathit{\sigma }}_{\mathrm{iso}}[\lambda ,{\lambda }_2(\lambda )],\\ & \hspace{0.17em}⟹\hspace{0.17em}f\hspace{0.17em}:\hspace{0.17em}\lambda \to {\lambda }_2.\end{array}$$3.3

A uniaxial extension along the principal direction ${\mathit{E}}_2$ can be computed analogously, with λ₂ = λ and g : λ → λ₁. The results are illustrated in figure 1b.

3.1.3. Simple shear

A simple shear deformation can be prescribed by the following deformation gradient:

$$\mathbf{F}={\mathit{e}}_1\otimes {\mathit{E}}_1+{\mathit{e}}_2\otimes {\mathit{E}}_2+{\mathit{e}}_3\otimes {\mathit{E}}_3+\gamma {\mathit{e}}_1\otimes {\mathit{E}}_2,$$3.4

where γ is the amount of shear along ${\mathit{E}}_1$. By using the plane stress condition (2.12), the Cauchy stress tensor is given by

$$\mathit{\sigma }(\gamma )=[\mathbb{I}-\mathbf{I}\otimes ({\mathit{e}}_3\otimes {\mathit{e}}_3)]\hspace{0.17em}:\hspace{0.17em}{\mathit{\sigma }}_{\mathrm{iso}}(\gamma ).$$3.5

The results are illustrated in figure 1c.

3.2. Parallel-plate test validation

The parallel-plate test is an experimental procedure commonly used to characterize the elastic and inelastic material response of polymeric materials. In this section, a FE model is developed to simulate this experimental procedure applied to an electrospun polymeric material. This example was selected as a validation of the proposed theoretical and numerical framework. For this purpose, we determined the material parameters of the polymeric material from cyclic uniaxial and biaxial extension tests and stress relaxation experiments. We then employ these parameters in the numerical simulation of the validation experiment. The uncertainty in the experimental measurements was taken into account by developing a probabilistic framework based on Bayesian inference [47], which will be the subject of a follow-up study. Here, we focus on the description of the validation experiment and the numerical simulation in which the uncertainty in the experimental data is transferred to the numerical model.

3.2.1. Experimental procedure

We measured the force–displacement curves of two different sample geometries under cyclic compressive loading. The ring-shaped samples were tested in the parallel-plate configuration with displacement control at constant velocity. Displacements between 0.25 and 0.75, defined as a fraction of the initial diameter of each sample d₀, were applied. The material was provided by the company Xeltis BV (Eindhoven, The Netherlands) and used as received. To prevent unintended movements during loading and unloading, the samples were fixed at the top and bottom with a thin strip (d ≈ 1 mm) of double-sided adhesive tape. The nominal dimensions of the tested samples and the detailed experimental loading protocols used during work are summarized in table 2. All mechanical tests were performed with a displacement-controlled triaxial testing device with a vertical stroke resolution of 0.04 μm [48]. The nominal force of the installed sensor was ±2 N with an accuracy class of 0.5%. To account for the dispersion in sample lengths (l₀ = 10–16 mm), the recorded forces were normalized to the actual length of each sample for further processing. The original device was modified and updated with a 3D-printed sample stage with a side length of 70 mm. The dimensions of the stage were chosen so large that there is no overhang during testing and to ensure a homogeneous deformation across the whole sample is guaranteed. All tests were performed under dry conditions at room temperature. A total of n = 6 samples were tested for each of the two diameters. A scheme of the experimental configuration is shown in figure 2a.

Table 2. Nominal dimensions and loading protocols of the parallel-plate tests.

Collapse

geometry
diameter (internal) d0	6 mm	23 mm
thickness t0	0.5 mm	0.5 mm
length l0	10 mm	10 mm
loading
load steps Δ/d0	0.25–0.50	0.25–0.50–0.75
loading rate $\dot{\mathrm{\Delta }}$	10 mm min−1	20 mm min−1
no. cycles	5	5

• Download figure
• Open in new tab
• Download PowerPoint

3.2.2. Finite element model

The FE model is defined by taking advantage of the double symmetry of the problem such that only a quarter of the geometry is modelled as a three-dimensional deformable body. A local cylindrical reference system is defined by the orthogonal basis ${({\mathit{E}}_i)}_{i=1}^3$, where E₁ and E₂ indicate the axial and circumferential directions, respectively. The metal plate is defined as an analytical surface and is associated with a rigid body reference node. Then, a contact interaction between the analytical surface, which acts as the target, and the outer surface of the cylinder, which acts as the contact, is defined using the contact pair algorithm of Abaqus Standard FE solver with surface-to-surface contact discretization. Contact properties are defined by a frictionless canonical pressure–overclosure relationship. The effect of glue between plate and specimen is simulated by imposing a non-separation behaviour on a strip with a width approximately equal to the width of the tape used in the experiments. This specification prevents nodes from loosening after contact and thus creating tensile forces during unloading. Boundary conditions are defined to prevent rigid body motion and enforce symmetry. The loading is assigned by imposing a vertical displacement on the rigid body reference node while constraining all other degrees of freedom and using an amplitude curve to define the loading history as in the experimental protocol. A sketch of the reference configuration of the FE model is shown in figure 2b. A preliminary mesh convergence study is performed by tracking the resultant force and maximum true principal strain at the largest diameter. The optimal mesh has a relative characteristic element size of h_el/d₀ = 0.03 and h_el/d₀ = 0.01 for d₀ = 6 mm and d₀ = 23 mm, respectively.

In order to validate the numerical model with the experiments, uncertainties related to the experimental geometries should also be considered. For the diameter of the cylindrical sample, we have assumed a uniform distribution within the limits d₀ ± 1 mm. It was assumed that the thickness follows a normal distribution with a mean of t₀ and a standard deviation of 0.01 mm (table 2). A total of 800 simulations were run for each diameter, randomly selecting geometric and material properties from the specified distributions. In order to reduce the computational time, a plane strain assumption in the axial direction was introduced, so that a two-dimensional model is obtained in the circumferential–radial plane (E₂, E₃) (figure 2b). A comparison between representative three-dimensional and plane strain solutions is reported in the electronic supplementary material, figure S7. Representative results of the simulations with the diameter d₀ = 23 mm are shown in figure 3, where the material parameters given in table 1 are used with mean geometric values. The normalized reaction force $\overline{F}=F/({\mu }_{\Vert }{d}_0)$ and the prescribed displacement were extracted from the reference node of the rigid plate and doubled to account for symmetry. Figure 4 shows the comparison with the experimental results, using only the data from the last unloading curve of the last loading step. The FE result is represented by the median within a 90% confidence interval. The experimental results are visualized as a mean along with the maximum and minimum measurements. The corresponding plots for the diameter d₀ = 6 mm can be found in the electronic supplementary material, figure S9.

• Download figure
• Open in new tab
• Download PowerPoint

• Download figure
• Open in new tab
• Download PowerPoint

4. Discussion

Polymeric scaffolds have received increased popularity in tissue engineering because of the innovative design paradigm they have introduced. However, the full potential of the in situ tissue engineering approach is still hampered by the large number of design variables that can potentially affect the final outcome in terms of functionality and safety of the replacement implant. In our research, we focused on the mechanical behaviour of the scaffold material before the degradation process and native tissue formation.

The new constitutive model [32] recently introduced by the authors encompasses the full set of nonlinear and inelastic effects observed in typical electrospun polymers used for the scaffolds in a continuum mechanics framework. As a further and fundamental step, we have proposed an efficient numerical implementation of the model for the FE method in this work. The code was developed for the commercial software Abaqus and can be adapted to other commercial FE packages, allowing a straightforward use in the biomedical industry. In addition, the numerical implementation was verified by comparing the FE results with semi-analytical solutions for elementary deformations and validated through a large-scale bending problem with contact.

The verification tasks were successfully completed, demonstrating the correctness of the numerical implementation. Overall, the constitutive behaviour observed in response to elementary deformations indicates a strongly anisotropic response of the scaffold material in the local axial–circumferential plane (figure 1). Furthermore, the expected rate-stiffening effect is also correctly captured by the model.

The validation problem was selected to allow large-scale application of the model. Furthermore, additional complexity is introduced by the kinematic nonlinearities arising from contact interactions. Qualitatively, the simulation can reproduce the large hysteresis in the force–displacement curves and the different features observed between the two experimental geometries. The quantitative comparison is also satisfactory (figure 4), although there are some differences worth discussing. Firstly, the selected example is very sensitive to geometrical effects. As a rough measure of the influence, the reaction force when bending a circular ring scales with the third power of the thickness in the linear elastic theory [49]. For this reason, we included a quantification of the uncertainties on the experimental geometry in the simulation. In addition, the different boundary-value problem of the validation with regard to material parameter identification could also explain part of the observed discrepancy. Material parameters were identified using elementary uniaxial and biaxial deformations with stretches up to 50%, but no experiment was specifically designed to test the compressive response of the material. On the contrary, the validation problem involves extensive compression, very large displacements, but only moderate deformations (electronic supplementary material, figure S8). These observations underscore the importance of an appropriate selection of experimental protocols, a constitutive modelling framework and a validation problem in the development of biomedical solutions. Despite the limitations, the proposed constitutive model and its numerical implementation show great potential for more complex simulations in the future.

Ethics

This work did not require ethical approval from a human subject or animal welfare committee.

Data accessibility

The source code of the Abaqus Standard user-material subroutine is available from: https://doi.org/10.3217/bhxxe-ng262. Additional derivations and results can be found in the electronic supplementary material available online.

The data are provided in electronic supplementary material [51].

Declaration of AI use

We have not used AI-assisted technologies in creating this article.

Authors' contributions

M.T.: conceptualization, formal analysis, methodology, software, visualization and writing—original draft; M.P.W.: conceptualization, formal analysis, methodology and writing—review and editing; M.P.K.: investigation, methodology and writing—original draft; M.R.-P.: methodology, project administration and writing—review and editing; N.G.: project administration and writing—review and editing; G.A.H.: funding acquisition, methodology, supervision and writing—review and editing.

All authors gave final approval for publication and agreed to be held accountable for the work performed therein.

Conflict of interest declaration

We declare we have no competing interests.

Funding

We gratefully acknowledge the financial support provided by the European Union’s Horizon 2020 programme for research and innovation (grant no. 101017523).