A computational model to predict cell traction-mediated prestretch in the mitral valve

Abstract

Prestretch is observed in many soft biological tissues, directly influencing the mechanical behavior of the tissue in question. The development of this prestretch occurs through complex growth and remodeling phenomena, which yet remain to be elucidated. In the present study it was investigated whether local cell-mediated traction forces can explain the development of global anisotropic tissue prestretch in the mitral valve. Towards this end, a model predicting actin stress fiber-generated traction forces was implemented in a finite element framework of the mitral valve. The overall predicted magnitude of prestretch induced valvular contraction after release of in vivo boundary constraints was in good agreement with data reported on valvular retraction after excision from the heart. Next, by using a systematic variation of model parameters and structural properties, a more anisotropic prestretch development in the valve could be obtained, which was also similar to physiological values. In conclusion, this study shows that cell-generated traction forces could explain prestretch magnitude and anisotropy in the mitral valve.

Keywords:

• Prestretch
• mitral valve
• cell-traction forces
• finite element method

1. Introduction

Prestretch is the stretch present in many biological tissues which, upon release of in vivo boundary constraints, leads to retraction. Prestretch distinguishes itself from growth in the sense that it describes geometrical changes while preserving tissue volume. This prestretch directly and indirectly influences the biomechanical behavior of the tissue in question. First, the presence of prestretch directly affects the apparent in vivo stiffness of soft biological tissues. In fact, Rausch and Kuhl (2013) showed that when studying the mechanical properties of the mitral valve, taking tissue prestretch into account could explain great discrepancies (up to three orders of magnitude) in reported stiffness values in the literature. Moreover, Cardamone et al. (2009) showed that a linear distribution of elastin prestretch along the arterial wall can explain uniform homeostatic stresses in vivo. Second, prestretch can influence growth and remodeling processes, leading to different structural adaptation in the long run in these tissues (Grenier et al. 2005; Mol et al. 2005).

For these reasons, it is evident that quantification of this prestretch is imperative when studying the load bearing capacity or adaptation of soft biological tissues. Unfortunately, accurately determining the degree of prestretch in vivo is often cumbersome, if possible at all, due to the complex geometries and loading conditions, and always demands sacrificing the tissue. Many have attempted to quantify prestretch in cardiovascular tissues using such invasive measurements. Already in 1975, Dobrin et al. observed prestretch-induced longitudinal retraction of arteries after excision (Dobrin et al. 1975). In addition, Chuong and Fung (1986), with their famous opening angle experiment, were the first to quantify residual stresses in arteries, even after prestretch was already alleviated. Also in in vitro tissue engineered heart valves, immediate retraction is observed upon release of the leaflets (Mol et al. 2005), which was quantified by van Vlimmeren et al. (2011). In a later study, Van Vlimmeren et al. (2012) showed that this tissue retraction process can be split in an immediate passive response, and an additional cell-mediated active retraction of the leaflets. Finally, Amini et al. (2012) quantified leaflet stretches in the mitral valve after excision from the heart, referenced to the in vivo loaded configuration at minimum left ventricular pressure, and showed different prestretch magnitudes in the circumferential and radial direction.

Although these studies are very insightful, they still require invasive measurements and give little information on the underlying mechanisms of how tissue prestretch arises in vivo. Computational models can potentially be used to gain a mechanistic understanding of the development of tissue prestretch on the one hand, and aid in accurately predicting tissue prestretch magnitude and (an)isotropy in vivo on the other hand. Different studies have proposed models in which prestretch was prescribed explicitly (Nagel and Kelly 2012; Soares et al. 2014) or assumed a priori (Cardamone et al. 2009). A similar study was also conducted by Rausch and Kuhl (2013), in which they prescribed different degrees of prestretch in the mitral valve and used an inverse finite element method with in vivo recorded membrane displacements and boundary conditions to determine the valve’s material properties. The predicted membrane stiffness turned out to vary four orders of magnitude depending on the applied level of prestretch. This demonstrates that an accurate estimation of prestretch is crucial to quantify the material properties of biological tissues. However, in this model it was assumed that prestretch is mapped homogeneously and isotropically over the valve, while Amini et al. (2012) suggested that mitral valve prestretch is rather anisotropic. Furthermore, the biological mechanism through which prestretch is induced was not addressed.

Tissue prestretch originates under the influence of various, complex growth and remodeling processes, which are only partly understood (Ambrosi et al. 2011). In this study we hypothesize that cells embedded in the tissue are primarily responsible for the development of tissue prestretch, by directly compacting fibers in a ‘hand-over-hand’ fashion on a local level (Meshel et al. 2005). On a global tissue level, these traction forces may lead to the development of in vivo prestretch. These so-called cellular traction forces are already known to influence tissue organization directly and indirectly. First, cell-generated forces are known to compact a tissue in the unconstrained direction, leading to an anisotropic alignment of existing fibers in the constrained direction (Brown et al. (1998); Meshel et al. (2005); Weidenhamer and Tranquillo (2013). Second, besides affecting the existing extracellular matrix, newly synthesized matrix components are deposited primarily in the direction in which the cells apply the largest amount of traction forces (Wang et al. 2003; Gealy et al. 2009).

In the present study we aim to understand how cell-mediated traction forces may lead to the development of anisotropic tissue prestretch in the mitral valve. Towards this end, a model is required which predicts the development of traction forces by cellular actin stress fibers. Different models have been developed which use physically-motivated remodeling laws to predict cellular actin stress fiber remodeling. These models rely on stress and strain homeostasis (Deshpande et al. 2006, 2007; Vernerey and Farsad 2011; Obbink-Huizer et al. 2014) or on a thermodynamic equilibrium (Foucard and Vernerey 2012; Vigliotti et al. 2016) to predict stress fiber assembly and dissociation in response to topological and mechanical cues. Loerakker et al. (2014) coupled the cell-mediated remodeling laws of Obbink-Huizer et al. (2014) to an algorithm for collagen remodeling, and showed that cellular contractility is a very important affector of remodeling in tissue engineered heart valves in the pulmonary position (Loerakker et al. 2016).

In this study the framework by Rausch et al. (Rausch and Kuhl 2013; Rausch et al. 2013) was adapted by incorporating the model of Loerakker et al. (2014, 2016) in order to understand if local cell-generated traction forces can explain the global prestretch magnitude and direction in the mitral valve. First, using the original model parameters and structural organization by Rausch et al. (2013), a similar prestretch magnitude as reported by Amini et al. (2012) was obtained. However, contrary to their findings, the obtained prestretch appeared to be rather isotropic. Second, the effect of a systematic investigation in both model parameters and structural organization on the developed prestretch was investigated using insightful virtual stress-stretch experiments and circular sheets. Finally, using these variations, a prestretch magnitude and anisotropy and resulting mitral valve leaflet contraction (18% circumferential, 22% radial) is predicted, in which structural organization and the relative stiffness of each direction were big determinants. The predicted retractions were comparable to the study of Amini et al. (2012).

2. Methods

2.1. Kinematics

In order to model prestretch development in the mitral valve, we used a multiplicative split of the total elastic deformation tensor ${\mathbf{F}}_e$, into a prestretch induced tensor ${\mathbf{F}}_p$, and a deformation induced gradient (e.g. applied pressure and due to deformation in other elements) F (Johnson and Hoger 1995; Rausch et al. 2013; Genet et al. 2015): (1) $${\mathbf{F}}_e=\mathbf{F}·{\mathbf{F}}_p$$(1) where the elastic tensor was split into a volumetric (${\mathbf{F}}_e^{\mathit{vol}}$) and isochoric (${\overline{\mathbf{F}}}_e$) part: (2) $${\mathbf{F}}_e={\mathbf{F}}_e^{\mathit{vol}}·{\overline{\mathbf{F}}}_e\text{with}{\mathbf{F}}_e^{\mathit{vol}}={\left(J_e\right)}^{\mathrm{1/3}}\mathbf{I}\text{and}{\overline{\mathbf{F}}}_e={\left(J_e\right)}^{-\mathrm{1/3}}{\mathbf{F}}_e$$(2) with unity tensor $\mathbf{I}$ and J_e= det(${\mathbf{F}}_e$). The elastic right Cauchy-Green tensor ${\mathbf{C}}_e$ and its relation to the isochoric part ${\overline{\mathbf{C}}}_e$, are given by: (3) $${\mathbf{C}}_e={\mathbf{F}}_e^T·{\mathbf{F}}_e={\left(J^e\right)}^{\mathrm{2/3}}{\overline{\mathbf{C}}}_e$$(3)

2.2. Material behavior of the mitral valve

The constitutive response of the tissue was modeled using the same approach as in Rausch and Kuhl (2013). Briefly, a hyperelastic strain energy density function characterized by one family of fibers was used: (4) $$\psi =U\left(J_e\right)+\overline{\psi }({\overline{I}}_{1e},{\overline{I}}_{4e})$$(4)

Here, the total strain energy density function ψ is split into a volumetric part $U(J^e)$, and an isochoric part $\overline{\psi }$ dependent on the isochoric elastic invariants ${\overline{I}}_{1e} \text{and} {\overline{I}}_{4e}$, related to the isotropic and anisotropic response of the tissue respectively: (5) $${\overline{I}}_{1e}={\overline{\mathbf{C}}}_e:\mathbf{I}\hspace{1em}\text{and} \hspace{1em}{\overline{I}}_{4e}={\overline{\mathbf{C}}}_e:\mathbf{N}$$(5) with structural tensor $\mathbf{N}={\mathbf{n}}_0\otimes {\mathbf{n}}_0$ with ${\mathbf{n}}_0$ a unit vector in the fiber direction of the reference configuration.

Similarly, the second Piola-Kirchhoff stress was additively decomposed into a volumetric and isochoric part: (6) $${\mathbf{S}}_e=2\frac{\partial \psi }{\partial {\mathbf{C}}_e}={\mathbf{S}}_e^{\text{vol}}+{\mathbf{S}}_e^{\text{iso}}$$(6) where the individual volumetric and isochoric parts are given by: (7) $${\mathbf{S}}_e^{\text{vol}}=2\frac{\partial U}{\partial {\mathbf{C}}_e}=2\frac{\partial U}{\partial J_e}\frac{\partial J_e}{\partial {\mathbf{C}}_e}=p_e{J}_e{\mathbf{C}}_e^{-1}$$(7) (8) $${\mathbf{S}}_e^{\text{iso}}=2\frac{\partial \overline{\psi }}{\partial {\mathbf{C}}_e}={\left(J_e\right)}^{-\mathrm{2/3}}{\mathbb{P}}_e:{\overline{\mathbf{S}}}_e$$(8) with (9) $${\overline{\mathbf{S}}}_e=2\frac{\partial \overline{\psi }}{\partial {\overline{I}}_{1e}}\mathbf{I}+2\frac{\partial \overline{\psi }}{\partial {\overline{I}}_{4e}}\mathbf{N},\hspace{1em}{\mathbb{P}}_e=\mathbb{I}-\frac{1}{3}{\mathbf{C}}_e^{-1}\otimes {\mathbf{C}}_e\hspace{1em}\text{and} \mathbb{I}=\frac{1}{2}\left[\mathbf{I}\overline{\otimes }\mathbf{I}+\mathbf{I}\underset{\_}{\otimes }\mathbf{I}\right]$$(9) with ${\left\{\alpha \overline{\otimes }\beta \right\}}_{\mathit{ijkl}}={\left\{\alpha \right\}}_{\mathit{ik}}{\left\{\alpha \right\}}_{\mathit{jl}}$ and ${\left\{\alpha \underset{\_}{\otimes }\beta \right\}}_{\mathit{ijkl}}={\left\{\alpha \right\}}_{\mathit{il}}{\left\{\alpha \right\}}_{\mathit{jk}}$. Assuming a plane stress condition and incompressibility, an explicit expression for p_e can be derived: (10) $${\mathbf{S}}_{e33}={\mathbf{S}}_{e33}^{\text{vol}}+{\mathbf{S}}_{e33}^{\text{iso}}=J_e{p}_e{\mathbf{C}}_{e33}^{-1}+{\mathbf{S}}_{e33}^{\text{iso}}=0\to p_e=-\frac{1}{J_e}\frac{{\mathbf{S}}_{e33}^{\text{iso}}}{{\mathbf{C}}_{e33}^{-1}}$$(10)

Finally, a push forward operation of the second Piola-Kirchhoff stress into the current configuration gives the Cauchy stress: (11) $$\mathbf{\sigma }=\frac{1}{J_e}{\mathbf{F}}_e·{\mathbf{S}}_e·{\mathbf{F}}_e^T$$(11)

2.3. Material model

To model the constitutive behavior of the ground matrix and fibrous tissue in the mitral valve, a hyperelastic Holzapfel model (Gasser et al. 2006; Rausch et al. 2013) with a single fiber family was used. Two decoupled terms were employed to describe the isotropic and anisotropic tissue response respectively: (12) $$\overline{\psi }=c_0({\overline{I}}_{1e}-3)+\frac{c_1}{2c_2}\left(\hspace{0.17em}\text{exp}\hspace{0.17em}\left(c_2{(\kappa {\overline{I}}_{1e}+(1-3\kappa ){\overline{I}}_{4e}-1)}^2\right)-1\right)$$(12) with constants $c_0,c_1,c_2$ and fiber dispersity κ. The first derivatives of this strain energy function with respect to the invariants (${\overline{\psi }}_i=\frac{\partial \overline{\psi }}{\partial {\overline{I}}_i}$) are given by: (13) $${\overline{\psi }}_1=c_0+c_1\kappa \left(\kappa {\overline{I}}_1+(1-3\kappa ){\overline{I}}_4-1\right)\hspace{0.17em}\text{exp}\hspace{0.17em}\left(c_2{(\kappa {\overline{I}}_1+(1-3\kappa ){\overline{I}}_4-1)}^2\right)$$(13) (14) $${\overline{\psi }}_4=c_1(1-3\kappa )\left(\kappa {\overline{I}}_1+(1-3\kappa ){\overline{I}}_4-1\right)\hspace{0.17em}\text{exp}\hspace{0.17em}\left(c_2{(\kappa {\overline{I}}_1+(1-3\kappa ){\overline{I}}_4-1)}^2\right)$$(14) which are used in Equation (9)(9) $${\overline{\mathbf{S}}}_e=2\frac{\partial \overline{\psi }}{\partial {\overline{I}}_{1e}}\mathbf{I}+2\frac{\partial \overline{\psi }}{\partial {\overline{I}}_{4e}}\mathbf{N},\hspace{1em}{\mathbb{P}}_e=\mathbb{I}-\frac{1}{3}{\mathbf{C}}_e^{-1}\otimes {\mathbf{C}}_e\hspace{1em}\text{and} \mathbb{I}=\frac{1}{2}\left[\mathbf{I}\overline{\otimes }\mathbf{I}+\mathbf{I}\underset{\_}{\otimes }\mathbf{I}\right]$$(9) . The second derivatives (${\overline{\psi }}_{\mathit{ij}}=\frac{{\partial }^2\overline{\psi }}{\partial {\overline{I}}_i\partial {\overline{I}}_j}$) are defined as: (15) $${\overline{\psi }}_{11}=c_1{\kappa }^2\left(1.2c_2{\left(\kappa {\overline{I}}_1+(1-3\kappa ){\overline{I}}_4-1\right)}^2\right)\hspace{0.17em}\text{exp}\hspace{0.17em}\left(c_2{(\kappa {\overline{I}}_1+(1-3\kappa ){\overline{I}}_4-1)}^2\right)$$(15) (16) $${\overline{\psi }}_{14}=c_1\kappa (1-3\kappa )\left(1.2c_2{\left(\kappa {\overline{I}}_1+(1-3\kappa ){\overline{I}}_4-1\right)}^2\right)\hspace{0.17em}\text{exp}\hspace{0.17em}\left(c_2{(\kappa {\overline{I}}_1+(1-3\kappa ){\overline{I}}_4-1)}^2\right)$$(16) (17) $${\overline{\psi }}_{44}=c_1{(1-3\kappa )}^2\left(1.2c_2{\left(\kappa {\overline{I}}_1+(1-3\kappa ){\overline{I}}_4-1\right)}^2\right)\hspace{0.17em}\text{exp}\hspace{0.17em}\left(c_2{(\kappa {\overline{I}}_1+(1-3\kappa ){\overline{I}}_4-1)}^2\right)$$(17) and are used for the tangent modulus, which was implemented according to Prot et al. (2007).

2.4. Cell-mediated prestretch

To implement the local influence of cellular traction forces on prestretch, a recently developed model by Loerakker et al. (Obbink-Huizer et al. 2014; Loerakker et al. 2014) predicting actin stress fiber-generated stresses was included into the model. It was assumed that cellular stress fibers locally apply a stress magnitude ${\sigma }_{\mathit{sf}}^i$ in two directions, one along the main fiber direction and one perpendicular to this direction. The magnitude of this stress is dependent on the Green-Lagrange strain and strain rate (${ϵ}^i$ and ${\dot{ϵ}}^i$ respectively) in each direction, and a scalar σ_max, representing the maximum amount of cellular contraction. (18) $${\sigma }_{\mathit{sf}}^i={\sigma }_{\mathit{max}}{f}_{ϵ}\left({ϵ}^i\right)f_{\dot{ϵ}}\left({\dot{ϵ}}^i\right)$$(18) with ${ϵ}^i$ being dependent on the deformation gradient $\mathbf{F}$ : (19) $${ϵ}^i=\frac{1}{2}({\left(F^i\right)}^2-1)$$(19) with $F^i=\sqrt{n_0·{\mathbf{F}}^T·\mathbf{F}·n_0}$, the stretch in the direction of the fibers. In this particular case, the strain rate dependent function $f_{\dot{ϵ}}({\dot{ϵ}}^i)$ from the original model by Loerakker et al. (2014) was omitted since we only consider stress fiber stresses at the minimum and maximum points of the dynamic load curves (see section 2.5), where $f_{\dot{ϵ}}({\dot{ϵ}}^i)=1$. In Ristori et al. (2016), the same method was adopted, and even in combination with an analytical method to approximate the results of the stress fiber remodeling algorithm based on simplified strain dynamics, the model predictions appeared to be hardly affected by using just a limited part of the strain history.

The function $f_{ϵ}$ consists of two parts, accounting for the active contraction, and passive strain-hardening response of the stress fibers respectively: (20) $$f_{ϵ}=f_{ϵ,a}+f_{ϵ,p}$$(20) (21) $$f_{ϵ,a}\left({ϵ}^i\right)=\hspace{0.17em}\text{exp}\hspace{0.17em}(-{({ϵ}^i/{ϵ}_0)}^2)$$(21) (22) $$f_{ϵ,p}\left({ϵ}^i\right)=\{\begin{array}{c}{({ϵ}^i/{ϵ}_1)}^2,\hspace{1em} \text{if}\hspace{1em}{ϵ}^i\ge 0\hfill \\ 0,\hspace{1em}\text{if}\hspace{1em}{ϵ}^i<0\hfill \end{array}$$(22)

It was assumed that the actin stress fibers compact the tissue (ground matrix and fibers) until the magnitude of the Cauchy stress of the tissue in that direction is equal to that of the stress fibers, if the stress in the tissue is lower than that in the stress fibers). In that case, a preferred value for the tissue stress is obtained: (23) $${\sigma }_{i,\mathit{pref}}={\sigma }_{\mathit{sf}}^i,\hspace{1em}\text{if}\hspace{1em}{\sigma }_i<{\sigma }_{\mathit{sf}}^i$$(23)

To fulfill equation (23)(23) $${\sigma }_{i,\mathit{pref}}={\sigma }_{\mathit{sf}}^i,\hspace{1em}\text{if}\hspace{1em}{\sigma }_i<{\sigma }_{\mathit{sf}}^i$$(23) , a local preferred prestretch magnitude ${\lambda }_{p,\mathit{pref}}^i$ in that direction (i is either circumferential or radial) is found using a bisectioning method (where ${\lambda }_{p,\mathit{pref}}^i$ is restricted to values $\ge$ 1). Next this preferred prestretch was gradually applied and updated in the model by means of a first order rate equation with constant ${\tau }_{\lambda }$ : (24) $$\frac{d{\lambda }_p^i}{\mathit{dt}}=\frac{1}{{\tau }_{\lambda }}({\lambda }_{p,\mathit{pref}}^i-{\lambda }_p^i)$$(24)

At each point in time, the prestretch tensor at each material point ${\mathbf{F}}_p$ (used in equation (1)(1) $${\mathbf{F}}_e=\mathbf{F}·{\mathbf{F}}_p$$(1) ) can be derived by assuming that the applied local prestretch is volume-conserving: $${\mathbf{F}}_p=\left(\begin{array}{ccc}{\lambda }_p^{\mathit{circ}}& 0& 0\\ 0& {\lambda }_p^{\mathit{rad}}& 0\\ 0& 0& \frac{1}{{\lambda }_p^{\mathit{circ}}{\lambda }_p^{\mathit{rad}}}\end{array}\right)$$

Note that the actin stress fiber induced prestretch tensor ${\mathbf{F}}_p$ is defined at a local level, which implies that it can be heterogeneous and vary between the individual integration points. In addition, ${\lambda }_p^{\mathit{circ}}$ is not necessarily equal to ${\lambda }_p^{\mathit{rad}}$ leading to a potentially anisotropic (local) prestretch. Together with the load-induced deformation gradient ($\mathbf{F}$), the total elastic tensor ${\mathbf{F}}_e$ can be determined using Equation (1)(1) $${\mathbf{F}}_e=\mathbf{F}·{\mathbf{F}}_p$$(1) .

2.5. Implementation

2.5.1. Prestretch in the mitral valve

The entire framework was implemented in the finite element package Abaqus (Dassault Systèmes Simulia Corp., Providence, RI, USA, version 6.14-1) using the user-defined subroutine UMAT, following the implementation described in Prot et al. (2007). The anterior mitral valve leaflet was created in Abaqus using 1035 linear S3 shell elements (Göktepe et al. 2010; Rausch et al. 2013), with a shear stiffness of 100 MPa and a thickness of 1 mm (Figure 1A). The leaflet chordae tendinae were modeled as Neo-Hookean rods with a stiffness of 20 MPa and a cross-sectional area of 1 mm². To account for leaflet anisotropy, a discrete field of circumferentially oriented fibers is mapped over the valve (Rausch et al. 2013).

Figure 1. A–C) Meshes of the simulations used, where the fiber orientation is indicated with black arrows. A) The mitral valve including the chordae tendinea with the fibers orientated circumferentially. B) Circular sheet with the fibers orientated uniformly in the horizontal direction. C) Circular sheet with the fibers orientated circumferentially. The displacement magnitude with respect to the starting configuration (A–C) at the maximally loaded state without prestretch in the mitral valve (D) and the circular sheet (F). E, G) The displacement magnitude in the circular sheet and the mitral valve at their prestretched and minimally loaded state.

In the first step of the analysis, the lowest pressure during the ovine cardiac cycle (Rausch et al. 2011) was applied at the ventricular side of the leaflet. Next the maximum pressure was gradually applied, together with displacement boundary conditions to the outer boundary nodes to accommodate for the deformation of the annulus (Figure 1D). In the following step, actin fiber-generated stresses, and consequent valvular prestretches (section 2.4) were allowed to develop. As cell traction-generated prestress tends to be a process which typically stabilizes within hours (Zhao et al. 2013; van Kelle et al. 2019), the time period for this step in the simulation was set at 6 hours, after which an equilibrium in the predicted prestretch was reached. The maximal left ventricular pressure and annulus deformation were maintained throughout this remodeling step. After this the pressure was lowered again to minimum left ventricular pressure. The reasoning behind this is that Amini et al. (2012) used this situation as the benchmark configuration to reference the reported ex vivo stretches (Figure 1E). Next the pressure and displacement boundary conditions, as well as the chordae tendinea were removed from the model while only fixing the valve in one middle node in all three dimensions to determine the unloaded configuration of the leaflet. Since cells are known to actively exert traction forces on their environment after release of boundary constraints (van Vlimmeren et al. 2011), one additional remodeling step was added to the analysis to account for this process. Finally, the leaflet contraction in the circumferential and radial directions were determined. In analogy with Amini et al. (2012), this was done by referencing four nodal coordinates in the midsection at the end of the simulation, to the coordinates of the same nodes in vivo in the prestretched valve at minimum left ventricular pressure.

2.5.2. Prestretch in circular sheets

In the model of the valve, there is a complex interplay between the predicted prestretch on the one hand, and structural and anisotropic properties on the other hand. To gain a mechanistic understanding of how fiber organization influences prestretch development, the proposed framework was also applied to a relatively simple model in the form of a circular sheet. This sheet was constructed using 1177 S3 shell elements, while being fully encastred at the edge. The simulations were performed with two fiber orientations: uniformly horizontal (Figure 1B) and fully circumferential (Figure 1C). The same simulation steps were undertaken as with the valve: first pressurization of the sheet (Figure 1F), followed by prestretch development, application of minimum left ventricular pressure (Figure 1G), release of all boundary constraints while fixing the middle node and finally measurement of the sheet’s contraction.

2.5.3. Material parameters

For both the circular sheets and the valve, the parameters associated with the stress generated by the stress fibers, and time constant ${\tau }_{\lambda }$ (equation 24(24) $$\frac{d{\lambda }_p^i}{\mathit{dt}}=\frac{1}{{\tau }_{\lambda }}({\lambda }_{p,\mathit{pref}}^i-{\lambda }_p^i)$$(24) ) were directly adopted from Obbink-Huizer et al. (2014) and Loerakker et al. (2016) (Table 1).

Table 1. Parameter values used in all simulations.

Model component	Parameter	Value	Reference
Cell-mediated traction	σmax	200 kPa	Obbink-Huizer et al. (2014)
	${\epsilon }_0$	0.12	Obbink-Huizer et al. (2014)
	${\epsilon }_1$	0.17 kPa	Obbink-Huizer et al. (2014)
	${\tau }_{\lambda }$	1h	Loerakker et al. (2016)

For the Holzapfel model, initially the original material parameters for 30% explicitly applied prestrain in Rausch et al. (2013) were used (Table 2). To gain insight in the separate effect of the isotropic and anisotropic parameters on the predicted prestretch, c₀ and c₁ were varied. To possibly obtain a more anisotropic prestretch first the isotropic parameter c₀ was decreased, second the anisotropic parameter c₁ was increased followed by a combinations of these two adjustments (Table 2).

Table 2. Holzapfel Parameters.

Simulation	c0 (kPa)	c1 (kPa)	c2 (–)	κ (–)
Original parameters	21840.0	116.7	17.5	0.05
Lower c0	10920.0	116.7	17.5	0.05
Higher c1	21840.0	11670.0	17.5	0.05
Lower c0 & higher c1	10920.0	11670.0	17.5	0.05

3. Results

In order to elucidate to what extent cell-mediated traction forces form an explanation for mitral valve prestretch, first the framework was applied to the mesh representing the anterior mitral valve leaflet, using the original material parameters by Rausch et al. (2013) (Table 2). After application of the maximum left ventricular pressure and annulus deformation, the elastic stretches are merely caused by the load induced part $\mathbf{F}$ (Figure 2, top row). In the following step the cellular traction stresses were predicted in the model. These clearly lead to the development of significant amounts of prestretch in both the circumferential and radial direction (Figure 2, bottom row), as the total elastic stretches are higher compared to the merely hemodynamically loaded valve (Figure 2, top row). After release of all boundary constraints, and the following final remodeling step in which actively exerted traction forces after valvular excision were modeled, the total local prestretch in the circumferential (${\lambda }_p^{\mathit{circ}}$ ) and radial directions (${\lambda }_p^{\mathit{rad}}$ ) are very comparable (Figure 3A and B). In line with the local prestretch, this leads to a global valvular retraction which is also rather isotropic, being 22% circumferentially, 21% radially (Figure 3D). These values agree well with retractions reported by Amini et al. (2012) (18% circumferentially, 31% radially). What is also important to mention is that there are remaining residual stresses present in the valve after excision (Figure 3C).

Figure 2. Elastic Green-Lagrange stretches in the circumferential (left column) and radial (right column) directions due to hemodynamic loading alone (pressure and annulus deformation) (top row) and after stress fiber-generated prestretch was developed (bottom row) for the mitral valve with the original material parameters (Table 2).

Figure 3. The mitral valve with the original material parameters (Table 2) after release of all boundary constraints, with the final predicted prestretch in the circumferential (A) and radial direction (B). C) The maximal principal residual stresses left in the valve after excision. D) Overlay of the excised valve over the original geometry (maximal principal elastic Green-Lagrange stretch), with black arrows indicating the percentile valvular circumferential (bottom) and radial (right) retraction referenced to the in vivo loaded state at minimum pressure.

Although the predicted valvular contractions seem physiological, the magnitudes in both directions are similar, while Amini et al. (2012) reported rather anisotropic retractions. To elucidate what the relative contribution of the currently used material parameters is on the predicted prestretch, the parameters c₀ and c₁ (accounting for the isotropic and anisotropic stiffness) were decreased and increased respectively. Note that κ, the term directly associated with anisotropy in the model, was not altered during the simulations. The rationale behind this was that κ was already close to 0, hence anisotropic, while the overall material response with the given material parameters was not strongly anisotropic. Therefore, studying relevant and large variations in the anisotropic material behavior, could only be achieved by varying the constants c₀ and c₁.

Before immediately applying these changes to the material parameters to the valve, a virtual biaxial tensile test was performed, in which the constitutive response of the parameter variations was tested. This was done by implementing the material model in MATLAB (Mathworks, Natick, MA, USA) and virtually applying physiological circumferential (1.1) and radial (1.31) stretches (Amini et al. 2012) (Figure 4). The original model parameters yield a rather linear response in both the circumferential and radial direction for the entire physiological stretch range in the mitral valve (Amini et al. 2012). Moreover, the gross material response is isotropic, indicating that the isotropic parameter c₀ is apparently very dominant. Even after lowering c₀, the obtained stresses in both directions are quite comparable, albeit considerably lower compared to the original parameter set. In contrast, a higher anisotropic parameter c₁ gives a non-linear, stiffer response in the circumferential direction, while it has little effect on the constitutive behavior in the radial direction compared to the original model parameters hence, the material behavior is also more anisotropic. Finally, the combined adjustments to c₀ and c₁ naturally also yield a more anisotropic material response, being more non-linear circumferentially, and less stiff and linear radially.

Figure 4. Stress-stretch response in the circumferential (A) and radial direction (B) for different material parameters (Table 2) for a biaxial test with ${\lambda }_{\mathit{circ}}=1.1$ and ${\lambda }_{\mathit{rad}}=1.31$ (Amini et al. 2012).

Besides altering the material parameters, we also sought to investigate the effect of fiber orientation on the predicted prestretch and resulting valvular contraction. Therefore, relatively simple circular meshes with horizontally (Figure 5, left column) and circumferentially oriented fibers (Figure 5, right column) were used, in combination with the previously mentioned adjustments to the material parameters. When using the original model parameters, in line with the valve also the circular sheets contract in an isotropic fashion for both fiber orientations (Figure 5A and B). When lowering the isotropic parameter c₀, the overall resistance against the cellular traction forces decreases, leading to an increased magnitude of contraction increases in both cases (Figure 5C and D). However, with an increased anisotropic parameter c₁, a horizontal fiber orientation shows a clear anisotropic contraction (Figure 5E), whereas in the circumferential case, again a clear isotropic shrinkage of the sheet is observed (Figure 5F). If we finally combine the adjustments to c₀ and c₁, the retraction is once more isotropic in the sheet with the circumferential alignment (Figure 5H), while it is anisotropic (with a higher magnitude compared to Figure 5E) in the horizontally aligned fibers (Figure 5G).

Figure 5. Maximal principal elastic Green-Lagrange stretch after release of all boundary constraints for the circular sheets with a horizontal (left column) and circumferential (right column) fiber orientation. The final configuration of the sheets is depicted over the original geometry, with black arrows indicating the horizontal (bottom) and vertical (right) retraction. Each row corresponds to a different material parameter set (Table 2).

The last parameter set (lower c₀ and higher c₁) yielded either an anisotropic or isotropic contraction in the circular sheets, depending on the fiber orientation. This set of parameters was also applied to the mesh of the valve to see if an anisotropic contraction (as for instance reported in Amini et al. 2012) could also be obtained. Again, after the stress fiber remodeling step, the resulting valvular prestretch leads to higher elastic stretches in both the circumferential and radial direction (Figure 6, bottom row) compared to the merely hemodynamically loaded state (Figure 6, top row). However, contrary to the original parameter set, in this case the predicted prestretch is more anisotropic, being lower in the circumferential compared to the radial direction (Figure 7A and B). Once more, there are residual stresses in the contracted valve (Figure 7C). Finally, the valvular contraction is more anisotropic then the previous case, being 18% in the circumferential direction and 22% in the radial direction.

Figure 6. Elastic Green-Lagrange stretches in the circumferential (left column) and radial (right column) directions due to hemodynamic loading alone (pressure and annulus deformation) (top row) and after stress fiber-generated prestretch was developed (bottom row) for the mitral valve with the adjusted material parameters (lower c₀ and higher c₁, Table 2).

Figure 7. The mitral valve with the adjusted material parameters (lower c₀ and higher c₁, Table 2) after release of all boundary constraints, with the final predicted prestretch in the circumferential (A) and radial direction (B). C) The maximal principal residual stresses left in the valve after excision. D) Overlay of the excised valve over the original geometry (maximal principal elastic Green-Lagrange stretch), with black arrows indicating the percentile valvular circumferential (bottom) and radial (right) retraction referenced to the in vivo loaded state at minimum pressure.

4. Discussion

In this study, we aimed to elucidate to what extend cell-mediated traction forces may explain the development of anisotropic tissue prestretch in the mitral valve. Towards this end, a model predicting actin stress fiber traction forces (Loerakker et al. 2016) was implemented in a finite element framework of mitral valve (Rausch and Kuhl 2013; Rausch et al. 2013). The predicted cellular traction forces on a local level indeed lead to significant amounts of global prestretch in the mitral valve in vivo. Therefore, the actual elastic stretches in the valve are higher compared to the merely hemodynamically loaded state (Figures 2 and 6). This highlights the importance of implementing prestretch when modeling the mechanical response of soft biological tissues such as the mitral valve. Upon release of all boundary constraints, using the original model parameters, the prestretch in the valve was rather isotropic, which led to an isotropic valvular retraction. The overall predicted magnitude of this prestretch induced valvular contraction that was in good agreement with data reported on valvular retraction after excision from the heart (Amini et al. 2012). However, the retraction in the study by Amini et al. (2012) was actually more anisotropic.

In an attempt to explain why the predicted contraction was rather isotropic, the model parameters were varied and the effect of fiber organization was investigated. In a circular sheet with uniformly horizontally oriented fibers, an increased anisotropic material response led indeed to an anisotropic contraction, being larger in the radial direction. In addition, the magnitude of contraction in both directions could be influenced by adjusting the isotropic parameter. These findings show that structural material parameters are important in determining the degree of predicted prestretch. However, also structural tissue organization plays a significant role. In similar simulations with a circumferential fiber organization in the circular sheets, the degree of contraction was always isotropic, regardless of the relative contribution of the anisotropic fibers to the constitutive behavior. Considering that the valve’s fiber organization is somewhat in between the extreme fiber organizations prescribed in the circular sheets, it is not surprising that the valvular leaflet contraction was more anisotropic than before using the new model parameters, albeit not to the same extent as in the sheets with the purely horizontally organized fibers. So even though locally the prestretch is rather anisotropic, on a global level this can lead to a more isotropic valvular contraction due to structural properties.

These results demonstrate that both the material model with associated parameters, and the structural organization are important contributors to the magnitude and anisotropy of the overall tissue prestretch. A well-educated choice for such a material model and parameters is therefore imperative, to accurately predict tissue prestretch. In the past, various material models have been used to model the constitutive response of the mitral valve (May-Newman and Yin 1998; Prot et al. 2007, 2010; Rausch et al. 2013; Lee et al. 2015, 2016). In this particular case, a Holzapfel model was chosen as this allowed for a systematic decoupled variation of the isotropic and anisotropic responses and rendered the most stable results.

The valve’s behavior was modeled using a single family of fibers arranged along the circumferential direction in a saddle horn-like fashion as used earlier by Rausch and Kuhl (2013), while Lee et al. (2016) showed that radially oriented collagen and elastin fibers are also present in the atrialis layers of the valve. A proper representation of the valve’s structural properties is thus crucial for proper prestretch predictions. A future study could implement the different layers and associated fiber organizations into the model, to obtain a more elaborate representation of the mitral valve. The fact that tissue organization is an important determinant of valvular prestretch could explain the fact that a loss of tissue organization during mitral valve repair, induces the potential development of a non-physiological prestretch that could impair normal valvular function (Braunberger et al. 2001; David et al. 1993, 2013; Gillinov et al. 2001).

In conclusion, this study demonstrates that cell traction forces can explain the formation of significant degrees of prestretch in the mitral valve and the resulting valvular retraction magnitude is in good agreement with data reported in literature. Potentially, this model can be used to predict tissue prestretch in a wide variety of other soft biological tissues.

Acknowledgements

We would like to thank Harkamaljot Kandail for his valuable input during discussions about the computational model. The authors furthermore acknowledge Frank Baaijens and Harkamaljot Kandail for reviewing the manuscript.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

We acknowledge the support from the Netherlands Cardiovascular Research Initiative (CVON 2012-01): The Dutch Heart Foundation, Dutch Federation of University Medical Centers, the Netherlands Organization for Health Research and Development and the Royal Netherlands Academy of Sciences.