Deterministic and Explicit: A Quantitative Characterization of the Matrix and Collagen Influence on the Stiffening of Peripheral Nerves Under Stretch

Abstract

The structural organization of peripheral nerves enables them to adapt to different body postures and movements by varying their stiffness. Indeed, they could become either compliant or stiff in response to the amount of external solicitation. In this work, the global response of nerves to axial stretch was deterministically derived from the interplay between the main structural constituents of the nerve connective tissue. In particular, a theoretical framework was provided to explicitly decouple the action of the ground matrix and the contribution of the collagen fibrils on the macroscopic stiffening of stretched nerves. To test the overall suitability of this approach, as a matter of principle, the change of the shape of relevant curves was investigated for changes of numerical parameters, while a further sensitivity study was performed to better understand the dependence on them. In addition, dimensionless stress and curvature were used to quantitatively account for both the matrix and the fibril actions. Finally, the proposed framework was used to investigate the stiffening phenomenon in different nerve specimens. More specifically, the proposed approach was able to explicitly and deterministically model the nerve stiffening of porcine peroneal and canine vagus nerves, closely reproducing ($R^2>0.997$) the experimental data.

1. Introduction

Peripheral nerves are crucial to connect sensory and motor organs to the central nervous system [1,2]. Motor nerves carry information towards organs, muscles, or other peripheral effectors, and also, sensory nerves link sensory receptors to the brain [3]. These complex structures are hierarchically organized: internal axons, together with Schwann cells, are protected by a loose connective tissue (i.e., the endoneurium), which was found to be mainly composed by I and II type collagen, macrophages, and mast cells and filled by endoneurial fluid [2,4,5]. Endoneurial components are, in turn, enveloped in fascicles by the epineurial sheet, which is formed by layers of perineural cells, I and II type collagen, together with elastic fibers [2,4,5]. Finally, fascicles are grouped together by a further layer of connective tissue (i.e., the epineurium) mainly composed by I and III type collagen, elastic fibers, and mast and fat cells [2]. Despite this richness of internal substructures, the main constitutive components of the connective tissues of nerves are elastin and collagen, which were found to be mainly axially oriented [4].

In addition, peripheral nerves are physical objects, which obey physical laws [6]. Indeed, they react to external mechanical stimuli trough a specific mechanical response, increasing their stiffness to keep axons’ integrity and to preserve endoneurial structures from longitudinal overstretch [7]. The knowledge and the prediction of this kind of response are important in medicine [8,9] and physical therapy [2,10,11], while they are strategic in modern and high-tech bioengineering fields, as in neural interfaces design [12,13,14].

Although the anatomic path of nerves involves branching points, their mechanical response to longitudinal stretch was closely related to the response of the rectilinear nerve trunks [7]. As a consequence, and for the sake of simplicity, different kinds of models were implemented to predict the rectilinear nerve trunk response to external stimuli. Indeed, from one side, peripheral nerves were modeled as isotropic materials, providing a deterministic elastic [15,16] or hyperelastic response [17,18,19,20,21]. From the other side, a stochastic iterative fibril-scale mechanical model was implemented to reproduce the straightening of wavy fibrils and to account for the effects of interfibrillar crosslinks on the overall properties of the tissue [22,23]. Furthermore, recent studies aimed at investigating the relationship between macroscopic tensile loading of nerves and micro-scale deformations [24].

Nevertheless, a deterministic and explicit continuum model to explore the influence of the ground matrix and collagen fibrils on the nerve stiffening is apparently not available in literature. As a consequence, in this work, such a theoretical framework is provided to quantify the action of both the ground matrix and collagen fibrils on the nerve stiffening.

In particular, the text is organized as follows: First, a closed-form strain energy function (SEF) is proposed to account for the actions of the matrix and fibrils. Then, the change of the shape of the stress/stretch curve is studied for variations of relevant numerical parameters, while a further sensitivity study is performed. Finally, the proposed framework is used to investigate the stiffening phenomenon in different nerve specimens coming from different animal models [19,25].

2. Methods

The peripheral nerve was modeled as a three-dimensional continuum body (elliptic cylinder), which was deformed into a final configuration $\mathsf{\Gamma }$ starting from a reference configuration ${\mathsf{\Gamma }}_r$. The position of each point of this solid was identified through the vectors $\mathbf{X}$ and $\mathbf{x}$, respectively related to the reference and the deformed configuration. As shown in Figure 1 below.

The deformation of the body was described through a time-dependent map $\mathbf{x}=\beta (\mathbf{X},t)$, where $\beta$ is an invertible and regular function. Unlike previous models [17,18,19,20,21,22,23], a unit vector field $\mathbf{M}\left(\mathbf{X}\right)$ was used to explicitly account for the presence of the collagen fibrils in the reference configuration ${\mathsf{\Gamma }}_r$ [26,27]. For the sake of simplicity, the collagen fibrils were modeled as axially oriented with no dispersion [4].

The response of the nerve was described through a strain energy function $\mathsf{\Psi }$ (per unit of volume), which was affected by the initial fibrils’ direction, as well as by the strain measure. Therefore, the mean Cauchy stress tensor within the nerve is written as:

$$\begin{array}{c}\hfill \sigma =J^{-1}\mathbf{F}\frac{\partial \mathsf{\Psi }(\mathbf{C},\mathbf{M}\otimes \mathbf{M})}{\partial \mathbf{F}}\end{array}$$

(1)

where $\mathbf{F}=Grad\mathbf{x}$ is the deformation gradient, $J=det\left(\mathbf{F}\right)$, and $\mathbf{C}={\mathbf{F}}^{\mathbf{T}}\mathbf{F}$ is the right Cauchy–Green strain tensor.

In particular, an invariant-based form of this energy is used:

$$\begin{array}{c}\hfill \mathsf{\Psi }=\mathsf{\Psi }(I_1,I_2,I_3,I_4,I_5)\end{array}$$

(2)

where, $I_1=tr\left(\mathbf{C}\right)$, $I_2=\frac{1}{2}[I_1^2-tr\left({\mathbf{C}}^{\mathbf{2}}\right)]$, $I_3=det\left(\mathbf{C}\right)$, $I_4=\mathbf{M}·\left(\mathbf{CM}\right)$, $I_5=\mathbf{M}·\left({\mathbf{C}}^{\mathbf{2}}\mathbf{M}\right)$. The Cauchy stress tensor is, then, better specified as:

$$\sigma =J^{-1}[2{\mathsf{\Psi }}_1\mathbf{B}+2{\mathsf{\Psi }}_2(I_1\mathbf{B}-{\mathbf{B}}^{\mathbf{2}})+2I_3{\mathsf{\Psi }}_3\mathbf{I}+2{\mathsf{\Psi }}_4\mathbf{m}\otimes \mathbf{m}+2{\mathsf{\Psi }}_5(\mathbf{m}\otimes \mathbf{Bm}+\mathbf{Bm}\otimes \mathbf{m})]$$

(3)

where $\mathbf{B}={\mathbf{FF}}^{\mathbf{T}}$ is the left Cauchy-Green tensor, ${\mathsf{\Psi }}_i=\frac{\partial \mathsf{\Psi }}{\partial I_i}$, $i=1,2,3,4,5$, $\mathbf{m}=\mathbf{FM}$. As in previous literature [7,19,20,21], the nerve was considered incompressible (i.e., the nerve volume remained the same, when the specimen underwent deformation), then $I_3=det\left(\mathbf{C}\right)=1$, and the pressure contribution within the Cauchy stress tensor was made explicit. Thus, the stress is written as:

$$\sigma =-p\mathbf{I}+2{\mathsf{\Psi }}_1\mathbf{B}+2{\mathsf{\Psi }}_2(I_1\mathbf{B}-{\mathbf{B}}^{\mathbf{2}})+2{\mathsf{\Psi }}_4\mathbf{m}\otimes \mathbf{m}+2{\mathsf{\Psi }}_5(\mathbf{m}\otimes \mathbf{Bm}+\mathbf{Bm}\otimes \mathbf{m})$$

(4)

where p is the hydrostatic pressure, while $\mathbf{I}$ is the unit tensor. The physical behavior of the nerve was assumed to result from the contribution of the main structural components of its connective tissue, that is the elastin matrix and the collagen fibers. Therefore, to account for both their action, the strain energy function is written as:

$$\begin{array}{c}\hfill \mathsf{\Psi }(I_1,I_4)={\mathsf{\Psi }}_{mat}\left(I_1\right)+{\mathsf{\Psi }}_{fib}\left(I_4\right)\end{array}$$

(5)

where ${\mathsf{\Psi }}_{mat}\left(I_1\right)$ is the energetic term due to the ground elastin matrix, while ${\mathsf{\Psi }}_{fib}\left(I_4\right)$ is the energetic term describing the collagen fibrils.

As a consequence, ${\mathsf{\Psi }}_2={\mathsf{\Psi }}_5=0$, since $\mathsf{\Psi }(I_1,I_4)$ is independent of $I_2$ and $I_5$, and the Cauchy stress tensor is further simplified as:

$$\sigma =-p\mathbf{I}+2{\mathsf{\Psi }}_1\mathbf{B}+2{\mathsf{\Psi }}_4\mathbf{m}\otimes \mathbf{m}$$

(6)

More specifically, the energetic contribution of the matrix is written as [28]:

$$\begin{array}{c}\hfill {\mathsf{\Psi }}_{mat}\left(I_1\right)=\sum _{i=1}^U{K}_i{(I_1-3)}^i\end{array}$$

(7)

where $U\in \mathbb{N}$. In contrast, the energetic contribution of the collagen fibrils is derived from:

$$\begin{array}{c}\hfill {\mathsf{\Psi }}_{fib}(I_4,A,D)=A[\mu (I_4,D)-1]\end{array}$$

(8)

where $\mu (I_4,D)=\mathsf{\Theta }\left[D{\xi }^2\left(I_4\right)\right]$ and $\xi \left(I_4\right)=(I_4-1)$. In addition, $K_1,A,D\in \mathbb{R}$ are numerical parameters, while $\mathsf{\Theta }\in C^{\left(n\right)}$, and ($n\ge 2$) is a real function, a continuum with its n-th order derivatives. Here, for the sake of simplicity, $\mathsf{\Theta }(\ast )=exp(\ast )$ [29,30], and $U=1$.

Since the nerve was axially stretched, its lateral surface was stress-free, and ${\sigma }_{xx}={\sigma }_{yy}=0$. Therefore, the contributions of the axial Cauchy stress due to the matrix and to the fibrils are written as ${\sigma }_{z{z}_{mat}}\left(\lambda \right)$ and ${\sigma }_{z{z}_{fib}}\left(\lambda \right)$, respectively. However, in the following, these contributions will be simply written as ${\sigma }_{mat}\left(\lambda \right)$ and ${\sigma }_{fib}\left(\lambda \right)$, dropping the $zz$ index related to the direction of stretch (z axis). Similarly, in the following, the longitudinal stretch is $\lambda =\frac{L_f}{L_0}$, where $L_0$ is the initial length of the specimen, while $L_f$ is the final length of the nerve after the axial stretching. Recalling this notation, the axial stress due to the action of the matrix is:

$$\begin{array}{c}\hfill {\sigma }_{mat}\left(\lambda \right)=2K_1({\lambda }^2-\frac{1}{\lambda })\end{array}$$

(9)

while the axial stress due to the fibrils’ action is:

$$\begin{array}{c}\hfill {\sigma }_{fib}\left(\lambda \right)=4AD\chi \left(\lambda \right)[\omega (\lambda ,D)-1]\end{array}$$

(10)

where $\chi \left(\lambda \right)={\lambda }^2-1$ and $\omega (\lambda ,D)=exp\left[D{\chi }^2\left(\lambda \right)\right]$. Similarly, their rates of change with stretch are:

$$\begin{array}{c}\hfill {\partial }_{\lambda }{\sigma }_{mat}\left(\lambda \right)=2K_1\left(2\lambda +\frac{1}{{\lambda }^2}\right)\end{array}$$

(11)

$$\begin{array}{c}\hfill {\partial }_{\lambda }{\sigma }_{fib}\left(\lambda \right)=8AD\lambda \{[1+2D{\chi }^2\left(\lambda \right)]\omega (\lambda ,D)-1\}\end{array}$$

(12)

while their second derivatives with respect to $\lambda$ are:

$$\begin{array}{c}\hfill {\partial }_{\lambda }^2{\sigma }_{mat}\left(\lambda \right)=4K_1\left(1-\frac{1}{{\lambda }^3}\right)\end{array}$$

(13)

$$\begin{array}{c}\hfill {\partial }_{\lambda }^2{\sigma }_{fib}\left(\lambda \right)=8AD\{[2D\chi \left(\lambda \right)P(\lambda ,\chi \left(\lambda \right),D)+1]\omega (\lambda ,D)-1\}\end{array}$$

(14)

where $P(\lambda ,\chi \left(\lambda \right),D)=4D{\lambda }^2{\chi }^2\left(\lambda \right)+\chi \left(\lambda \right)+6{\lambda }^2$. As a consequence, the total amount of stress within the nerve, together with its derivatives, are written as:

$$\begin{array}{c}\hfill {\sigma }_{tot}\left(\lambda \right)={\sigma }_{mat}\left(\lambda \right)+{\sigma }_{fib}\left(\lambda \right)\end{array}$$

(15)

$$\begin{array}{c}\hfill {\partial }_{\lambda }{\sigma }_{tot}\left(\lambda \right)={\partial }_{\lambda }{\sigma }_{mat}\left(\lambda \right)+{\partial }_{\lambda }{\sigma }_{fib}\left(\lambda \right)\end{array}$$

(16)

$$\begin{array}{c}\hfill {\partial }_{\lambda }^2{\sigma }_{tot}\left(\lambda \right)={\partial }_{\lambda }^2{\sigma }_{mat}\left(\lambda \right)+{\partial }_{\lambda }^2{\sigma }_{fib}\left(\lambda \right)\end{array}$$

(17)

Furthermore, the partial amount of nerve stiffening due to each component (matrix and fibrils), as well as the total amount of nerve stiffening were measured through the extrinsic curvature of each stress/stretch curve. In particular, the dimensionless extrinsic curvature due to the matrix-derived stress/stretch curve, accounting for the stiffening due to the ground matrix action, is written as:

$$\begin{array}{c}\hfill {\widehat{\kappa }}_{mat}\left(\lambda \right)=\frac{{\partial }_{\lambda }^2{\widehat{\sigma }}_{mat}}{{[1+{\left({\partial }_{\lambda }{\widehat{\sigma }}_{mat}\right)}^2]}^{3/2}}\end{array}$$

(18)

where:

$$\begin{array}{c}\hfill {\widehat{\sigma }}_{mat}=\frac{{\sigma }_{mat}}{6K_1}.\end{array}$$

(19)

In contrast, the extrinsic curvature of the fibril-derived stress/stretch curve was studied through two different formulas. The first one is written as:

$$\begin{array}{c}\hfill {\tilde{\kappa }}_{fib}\left(\lambda \right)=\frac{{\partial }_{\lambda }^2{\tilde{\sigma }}_{fib}}{{[1+{\left({\partial }_{\lambda }{\tilde{\sigma }}_{fib}\right)}^2]}^{3/2}}\end{array}$$

(20)

where:

$$\begin{array}{c}\hfill {\tilde{\sigma }}_{fib}=\frac{{\sigma }_{fib}}{A}\end{array}$$

(21)

and accounts for the nerve stiffening due to the collagen fibrils’ action, while the second one is written as:

$$\begin{array}{c}\hfill {\widehat{\kappa }}_{fib}\left(\lambda \right)=\frac{{\partial }_{\lambda }^2{\widehat{\sigma }}_{fib}}{{[1+{\left({\partial }_{\lambda }{\widehat{\sigma }}_{fib}\right)}^2]}^{3/2}}\end{array}$$

(22)

where:

$$\begin{array}{c}\hfill {\widehat{\sigma }}_{fib}=\frac{{\sigma }_{fib}}{6K_1}\end{array}$$

(23)

and measures the amount of the nerve stiffening due to the collagen action with reference to the stiffening due to the matrix action. In addition, the total dimensionless stress in the nerve is written as:

$$\begin{array}{c}\hfill {\widehat{\sigma }}_{tot}=\frac{{\sigma }_{mat}+{\sigma }_{fib}}{6K_1}\end{array}$$

(24)

while the total dimensionless curvature of the stress/stretch curve, accounting for the coupled action of matrix and fibrils, is defined as:

$$\begin{array}{c}\hfill {\widehat{\kappa }}_{tot}\left(\lambda \right)=\frac{{\partial }_{\lambda }^2{\widehat{\sigma }}_{tot}}{{[1+{\left({\partial }_{\lambda }{\widehat{\sigma }}_{tot}\right)}^2]}^{3/2}}\end{array}$$

(25)

Once having defined the previous theoretical framework, the range of variation of Equations (9) and (10) was explored by changing their main numerical parameters (i.e., $K_1,A,D$) through the auxiliary set $Q=\{0,0.1,0.5,1,5,10,50\}$, to test their suitability for reproducing the classic evolution of the stress/stretch curve [2,7]. Moreover, to assess their quantitative dependence on the numerical parameters, a sensitivity study was performed through the index [31,32]:

$$\begin{array}{c}\hfill s[q\left(\lambda \right),p_i]=\frac{q_{out}\left(\lambda \right)-q_{ref}\left(\lambda \right)}{q_{out}\left(\lambda \right)}\end{array}$$

(26)

where q is the studied quantity, $q_{out}$ is the output value of q derived from the variation of the chosen parameter $p_i$, while $q_{ref}$ is the reference value, when all the parameters are set to 1. The further the resulting value of Equation (26) was far from zero, the more the effects of the studied parameter affected the quantity q.

Finally, theoretical predictions were compared to experimental data collected from five consecutive extensions of a porcine peroneal nerve [19] and for one extension of a canine left vagus nerve [25]. More specifically, Equation (15) was used to reproduce the experiments, while numerical parameters were optimized through a non-linear procedure (quasi-Newton algorithm, Scilab; Scilab Enterprises S.A.S 2015), allowing the $R^2$ statistic to be maximized and residual plots minimized for each extension. Once having obtained the relevant parameters ($K_1,A,D$) for the real experiments, Equations (18)–(25) were used to explore the action of the ground matrix and the collagen fibrils on the stiffening of the nerve during axial stretch.

3. Results

To explore whether, as a matter of principle, the previous framework was able to reproduce the classic evolution of the nerve stiffening phenomenon, the range of variation of the main theoretical curves was studied. First, the range of variability of Equation (9), accounting for the stress due to the matrix contribution (Figure 2a), was investigated for increasing values of $K_1$. In particular, the values of ${\sigma }_{mat}(\lambda ,1+\mathsf{\Delta }{K}_1)$ were computed for $\mathsf{\Delta }{K}_1\in Q$. For a small stretch (i.e., $\lambda \to 1$), Equation (9) reduced to ${\sigma }_{mat}=6K_1(\lambda -1)$, so the initial slope of the curve was ${\partial }_{\lambda }{\sigma }_{mat}(1,K_1)=6K_1$.

Then, Equation (10), accounting for the stress derived from the collagen fibrils’ contribution, was investigated for $D=1$ and $A=1+\mathsf{\Delta }A$, $\mathsf{\Delta }A\in Q$. Larger values of A resulted in larger values of ${\sigma }_{fib}(\lambda ,A,1)$ (as shown in logarithmic scale in Figure 2b), while for $\lambda \to 1$, its slope was ${\partial }_{\lambda }{\sigma }_{fib}(\lambda ,A,1)=0,\forall A\in \mathbb{R}$. The contribution of fibrils was also investigated for $A=1$ and for different values of $D=1+\mathsf{\Delta }D$, $\mathsf{\Delta }D\in Q$, through the evolution of the natural logarithm of ${\sigma }_{fib}(\lambda ,1,D)$ (Figure 2c). More specifically, the larger the value of D was, the more the value of fibril stress increased, while for $\lambda \to 1$, its slope was ${\partial }_{\lambda }{\sigma }_{fib}(\lambda ,1,D)=0,\forall D\in \mathbb{R}$.

In addition, Equation (18), accounting for the stiffening due to the matrix action ${\widehat{\kappa }}_{mat}\left(\lambda \right)$, monotonically increased for $1\le \lambda \le 1.54$, while after a maximum ${\widehat{\kappa }}_{mat}\left(1.54\right)=0.133$, it monotonically decreased for $1.54\le \lambda \le 2$ (Figure 2d). For $\lambda \to 1$, its slope was ${\partial }_{\lambda }{\widehat{\kappa }}_{mat}\left(1\right)=2^{-1/2}$.

Similarly, Equation (22), accounting for the stiffening due to the action of fibrils, for $D=1$ (i.e., ${\tilde{\kappa }}_{fib}(\lambda ,1)$), had a fixed bell-like shape (Figure 2e), with a maximum ${\tilde{\kappa }}_{fib}(1.06,1)=11.317$. Its full width at half maximum was FWHM$\left({\tilde{\kappa }}_{fib}\right)\simeq 0.082$, while for $\lambda \to 1$, its slope was ${\partial }_{\lambda }{\tilde{\kappa }}_{fib}\left(1\right)=192$. On the contrary, the same Equation (22) for different values of $D=1+\mathsf{\Delta }D$, $\mathsf{\Delta }D\in Q$ (i.e., ${\tilde{\kappa }}_{fib}(\lambda ,D)$), had a variable bell-like shape (Figure 2f). The larger the value of D was, the more the value of the peak increased, while the more the value of the FWHM$\left({\tilde{\kappa }}_{fib}\right)$ decreased. In other words, for $D\to \infty$, $max\left\{{\tilde{\kappa }}_{fib}(\lambda ,D)\right\}\to \infty$, while FWHM$\left({\tilde{\kappa }}_{fib}\right)\to 0$. However, these curves had the following characteristic: ${lim}_{{\lambda }_{max}\to \infty }{\int }_1^{{\lambda }_{max}}{\tilde{\kappa }}_{fib}(\lambda ,D)d\lambda =1$, $\forall D\in \mathbb{R}$. Furthermore, for $\lambda \to 1$, its slope quadratically increased with D, since it was ${\partial }_{\lambda }{\tilde{\kappa }}_{fib}\left(D\right)=192D^2$.

To better quantify the influence of each parameter, a sensitivity study was also performed on each contribution. More specifically, the sensitivity of ${\sigma }_{mat}$ with respect to changes in $K_1$ was tested in Figure 3a. Since it was found that $s({\sigma }_{mat},K_1)=\frac{K_1-1}{K1}$, the larger the value of $K_1=1+\mathsf{\Delta }{K}_1$ ($\mathsf{\Delta }{K}_1\in Q$) was, the more the value of its value converged towards one.

Again, the effects of the changes in A or D were investigated on the stress functions due to the collagen fibrils. In this case as well, it was found that $s({\sigma }_{fib},A)=\frac{A-1}{A}$, (Figure 3b), for different values of $A=1+\mathsf{\Delta }A$, $\mathsf{\Delta }A\in Q$. Thus, the larger the value of A was, the more the value of $s({\sigma }_{fib},A)$ converged towards one. Similarly, the evolution of $s({\sigma }_{fib},D)$ is shown in Figure 3c for different values of $D=1+\mathsf{\Delta }D$, $\mathsf{\Delta }D\in Q$. The larger the value of D was, the quicker the value of $s({\sigma }_{fib},D)$ converged towards one. On the contrary, the larger the value of D was, the quicker the function $ln\left[s({\tilde{\kappa }}_{fib},D)\right]$ diverged, as shown in Figure 3d.

Once having tested the suitability of the provided framework, Equation (15) was used to study the experimental response of a porcine peroneal nerve (five consecutive elongations up to $\lambda =1.08$) [19]. Stress data are reported as a function of stretch (Figure 4) and quantitatively compared to theoretical predictions through the $R^2$ statistic and the computation of residuals ($\mathsf{\Delta }\sigma$). This comparison resulted in $R^2=0.9979$ for the first elongation (Figure 4a), $R^2=0.9976$ for the second one (Figure 4b), $R^2=0.9987,0.9989,0.9989$ (Figure 4c–e). Similarly, residuals ranged between $-0.056$ kPa and $0.091$ kPa for the first elongation (Figure 4f), while they ranged between $-0.105$ kPa and $0.086$ kPa for the second one (Figure 4g). Moreover, residuals spanned from $-0.078$ kPa to $0.078$ kPa (Figure 4h), from $-0.077$ kPa to $0.083$ kPa (Figure 4i), and from $-0.074$ kPa to $0.108$ kPa (Figure 4j), for the third, fourth, and fifth elongation, respectively. Finally, relevant quantities (mean ± standard deviation) in Equation (15) were quantified and resulted in $K_1=3.018\pm 0.130$ kPa, $A=0.014\pm 0.001$ kPa, and $D=66.893\pm 1.527$.

Moreover, the dimensionless stress and curvature were investigated. In particular, the dimensionless stress due to the matrix contribution (${\widehat{\sigma }}_{mat}$) monotonically increased from zero ($\lambda =1$) to $0.080$ ($\lambda =1.08$, Figure 5a). The slope of ${\widehat{\sigma }}_{fib}$ was (mean ± standard deviation) $13.011\pm 1.797$ for $\lambda =1.08$, while it was zero by definition at $\lambda =1$ (Figure 5b). Similarly, the slope (mean ± standard deviation) of the total dimensionless stress ${\widehat{\sigma }}_{tot}$ was $14.017\pm 1.797$, at $\lambda =1.08$, while, by definition, the slope was one at $\lambda =1$ (Figure 5c). The dimensionless curvature ${\widehat{\kappa }}_{mat}$ monotonically increased from zero ($\lambda =1$) up to $0.048$ ($\lambda =1.08$) with a slope of about $0.59$ (Figure 5d). On the contrary, the dimensionless curvature due to the fibrils’ action ${\widehat{\kappa }}_{fib}$ had a bell-like shape (Figure 5e). In particular, this curve reached its maximum value (mean± standard deviation) $28.771\pm 0.727$ at $\lambda =1.035\pm 0.001$. The median, maximum, and minimum values of $max\left\{{\widehat{\kappa }}_{fib}\right\}$ were, respectively $29.085,29.141,27.471$, as shown in Figure 5f, together with the shape of their distribution. Finally, ${\widehat{\kappa }}_{tot}$ is studied in Figure 5g. In this case as well, the curve had a bell-like shape and reached the maximum value of $6.840\pm 0.176$ at $\lambda =1.035\pm 0.002$. The median, maximum, and minimum values of $max\left\{{\widehat{\kappa }}_{tot}\right\}$ were respectively $6.953,6.970,6.568$, as shown in Figure 5h, together with their distribution shape.

The proposed framework was also used to explore experimental data derived from the axial stretching of a canine left vagus nerve [25]. In this case, Equation (15) with $K_1=9.731$ kPa, $A=5.277$ kPa, $D=10.171$, was able to closely reproduce the experiments ($R^2=0.9998$, Figure 6a), while residuals ($\mathsf{\Delta }\sigma$) spanned from $-0.316$ to $0.165$ (Figure 6b). The dimensionless stress due to the action of the matrix linearly increased from zero ($\lambda =1$) to $0.088$ ($\lambda =1.089$), and the slope of this curve was $58.39$ for $\lambda \to 1$, as shown in Figure 6c. On the contrary, the dimensionless stress due to fibrils increased in a non-linear way from zero to $0.089$, and the slope of this curve increased from zero ($\lambda =1$) to $4.676$ ($\lambda =1.089$) (Figure 6d). Similarly, the total dimensionless stress increased in a non-linear way from zero to $0.177$, and its slope spanned from $1.000$ ($\lambda =1$) to $12.384$ ($\lambda =1.089$), as shown in Figure 6e. Furthermore, the dimensionless curvature ${\widehat{\kappa }}_{mat}$ monotonically increased from zero ($\lambda =1$) up to $0.052$ ($\lambda =1.089$) (Figure 6f), while the dimensionless curvature due to the fibrils’ action ${\widehat{\kappa }}_{fib}$ had a bell-like shape (Figure 6g). In this case, the maximum value of the curve was $max\left\{{\widehat{\kappa }}_{fib}\right\}=25.732$, which was reached at $\lambda =1.0565$. Finally, ${\widehat{\kappa }}_{tot}$ was studied. This curve had a bell-like shape, and the maximum value was $max\left\{{\widehat{\kappa }}_{tot}\right\}=6.183$ at $\lambda =1.058$, as shown in Figure 6h.

4. Discussion

To tackle nerve neuropathies deriving from different kinds of pathologies (i.e., diabetes, scleroderma) [33], dietary deficiency, or exposure to ionizing radiation (e.g., neuromalacia) [34], an explicit understanding of the main features of the peripheral nerve response to stretch is needed. In particular, a computationally light and physically-based theoretical approach, linking in a deterministic and explicit way the main physical constituents of the connective tissue to the macroscopic response of the peripheral nervous tissue, could be strategic to make quantitative and reliable predictions.

Since peripheral nerves are anatomically complex (i.e., they have branching points), but their overall response is closely related to the response of their main rectilinear tract (which is the main structural segment), in this work, they were modeled as continuum rectilinear transverse isotropic bodies. Indeed, although peripheral nerves are hierarchically organized, ultrastructural studies revealed that the main constituents of the connective tissue are a ground matrix of elastin and a longitudinally-oriented fibrous network of collagen [4]. As a consequence, the strain energy function, ruling the nerve behavior, was dominated by ${\mathsf{\Psi }}_{mat}$ and ${\mathsf{\Psi }}_{fib}$, that is by the energetic contributions of the matrix and fibrils, respectively.

In particular, the energetic contribution of the elastic ground matrix was modeled through a first order Yeoh-like strain energy function [28], while the energetic contribution of the collagen fiber network was modeled through Equation (10), where, for the sake of simplicity, the function $\omega (\lambda ,D)$ was a convex energy function, traditionally used to model the behavior of soft tissues [29,30]. The superimposition of these energetic contributions was able to closely reproduce the experimental evolution of the mean Cauchy stress.

First, to qualitatively test the suitability of the provided approach, changes in numerical parameters $(K_1,A,D)$ were investigated with respect to their effectiveness in reshaping both the ${\sigma }_{mat}$ and ${\sigma }_{fib}$ curves. In particular, with reference to ${\sigma }_{mat}$, the more $K_1$ increased, the steeper the stress curve was, due to the action of the ground matrix. Similarly, the steepness of ${\sigma }_{fib}$ increased in a linear way by increasing the A value, since it resulted in a proportional magnification of the curve. On the contrary, D affected ${\sigma }_{fib}$ in a highly non-linear way.

In the same way, the dimensionless extrinsic curvature ${\widehat{\kappa }}_{mat}$ showed that the contribution of the ground matrix to the nerve stiffening was limited with stretch, while ${\tilde{\kappa }}_{fib}$ heavily affected the stiffening process. As a consequence, three parameters were enough to provide a variety of shapes, allowing experimental trends to be correctly reproduced by theoretical predictions, while in the literature, more parameters were needed to account for the collagen action through an iterative and probabilistic approach [22,23].

In addition, the sensitivity of ${\sigma }_{mat}$ and ${\sigma }_{fib}$ and ${\tilde{\kappa }}_{fib}$ towards changes in numerical parameters was explored. In particular, the sensitivity of the matrix stress function to the variation of $K_1$ and the sensitivity of the fiber stress function to the variation of the A parameter (with respect to the unit values of the other constants) converged towards zero for increasing values of $K_1$ and A. On the contrary, the sensitivity of the fibril stress function (${\sigma }_{fib})$ to the variation of the D parameter quickly reached the steady value of one. The higher the value of D was, the faster the steady state was reached. As a consequence, the parameter D affected ${\sigma }_{fib}$ more extensively than the parameter A. In contrast, A affected ${\sigma }_{fib}$ exactly as $K_1$ affected ${\sigma }_{mat}$. The sensitivity of ${\tilde{\kappa }}_{fib}$ to variations of the parameter D resulted in a highly non-linear and diverging line. The more the value of D increased, the more the line diverged. In other words, the value of D heavily affected the curvature of the stress due to the collagen action.

The parameter $K_1$ was, then, physically related to the apparent stiffness of the nerve for small strain (in this case, the Young modulus was $6K_1$), while the parameter A controlled the linear part of stress increment due to the collagen action. Furthermore, the parameter D affected the curve steepness and controlled both the delay (with respect to the reference configuration) and the velocity of the nerve stiffening process. Therefore, it was shown, as a matter of principle, that the proposed deterministic and explicit framework, exploiting closed-form equations with three main parameters, was able to reproduce the classic nerve stiffening behavior [2,7].

As a consequence, first, this framework was used to explore the stiffening behavior of real specimens. Indeed, with reference to the porcine peroneal nerve, theoretical predictions were able to closely reproduce the experiments ($R^2\ge 0.997$), while residual plots were analyzed to avoid numerical bias (e.g., over-fitting) keeping the value of the $R^2$ statistic artificially high. All residuals ($\mathsf{\Delta }\sigma$), between theoretical predictions and experimental values, were within $0.1$ kPa; that is, the error was about $2.5\%$ in the worst case, confirming that experiments were reproduced with confidence. To further explore the response of the porcine nerve, dimensionless stress and curvature were used, since the dimensionless stress was able to quantify the evolution of stress with respect to the stiffness of nerve for small strain, while the dimensionless curvature was used to measure the nerve stiffening, since it mapped where and how much the dimensionless stress evolved non-linearly because of the action of collagen fibrils.

More specifically, the slope of ${\widehat{\sigma }}_{tot}$ at $\lambda =1.08$ was on average 14 times greater than the slope at the reference configuration (i.e., at $\lambda =1$), while the slope of ${\widehat{\sigma }}_{mat}$ was practically constant. In addition, ${\widehat{\kappa }}_{mat}$ was small with respect to ${\widehat{\kappa }}_{fib}$ (i.e., $max\left\{{\widehat{\kappa }}_{fib}\right\}/max\left\{{\widehat{\kappa }}_{mat}\right\}\simeq 599.39$). These results not only suggested that the stiffening process was dominated by the action of the collagen fibrils, confirming previous literature [2,7,24], but also provided a quantification of the amount of the stiffening phenomenon (i.e., through the change of the slope by 14 times), as well as a quantification of how much the fibrils’ action was able to overcome the contribution of the ground matrix (i.e., about two orders of magnitude). On the contrary, the maximum dimensionless curvature of ${\widehat{\sigma }}_{tot}$ was lower than the maximum curvature of ${\widehat{\sigma }}_{fib}$ (i.e, $max\left\{{\widehat{\kappa }}_{fib}\right\}/max\left\{{\widehat{\kappa }}_{tot}\right\}\simeq 4.21$). This result suggested a possible lowering action of the ground matrix on the huge stiffening action of the collagen fibrils.

Furthermore, since the proposed framework was not intrinsically connected to any kind of animal or human model, it is, as a matter of principle, suitable to explore any kind of nerve. To show this, a different specimen, coming from a different animal model, was investigated. Indeed, with reference to a canine left vagus nerve, theoretical predictions were able to closely reproduce the experiments ($R^2\ge 0.999$), while all residuals were within $0.32$ kPa; that is, the error was less than $1.5\%$. Furthermore, the slope of ${\widehat{\sigma }}_{tot}$ at $\lambda =1.08$ was on average 12 times greater than the slope at the reference state (i.e., at $\lambda =1$), while the slope of ${\widehat{\sigma }}_{mat}$ was practically constant. In addition, ${\widehat{\kappa }}_{mat}$ was small with respect to ${\widehat{\kappa }}_{fib}$ (i.e., $max\{{\widehat{\kappa }}_{fib}/\}max\left\{{\widehat{\kappa }}_{mat}\right\}\simeq 494.84$). Therefore, again, these results were able to quantify the stiffening phenomenon in a canine vagus nerve (change of the dimensionless total stress slope of about 12 times) and showed that the action of collagen fibrils dominated the nerve stiffening phenomenon (of about two order of magnitude). Moreover, since ${\widehat{\kappa }}_{tot}$ was lower than ${\widehat{\kappa }}_{fib}$ and there was $max\left\{{\widehat{\kappa }}_{fib}\right\}/max\left\{{\widehat{\kappa }}_{tot}\right\}\simeq 4.16$, also in this case, a lowering action of the matrix with respect to the stiffening action of the collagen fibrils was suggested.

Curiously, for both specimens $max\left\{{\widehat{\kappa }}_{fib}\right\}/max\left\{{\widehat{\kappa }}_{tot}\right\}\simeq 4$, indicating the existence of common features between different species. More specifically, the analysis of dimensionless stress suggested that collagen fibrils were the main common structural features supporting the strain stiffening behavior. On the contrary, the analysis of the dimensionless curvature suggested a kind of common lowering action of the ground matrix on the huge stiffening due to the collagen fibrils’ action. In addition, although the model made abstraction of the real location of the collagen inside the nerve, it supported the hypothesis that the main structural contribution was performed by connective tissues in epineurium and perineural sheets [7].

Finally, since the provided model was general, it was able to investigate the mechanical response of different nerves coming from different animal models (i.e., suine peroneal and canine vagus nerves). In both cases, the experimental evolution of stretch was closely reproduced ($R^2\ge 0.99$) with small residuals. In other words, the provided framework was predictive, by definition, for the tested specimens, since $R^2$ statistics were high and residuals small. However, it seems to be predictive also for human cadaveric specimens [7] and for nerves coming from other animal species (i.e., rabbit [35], lobster, aplysia, [19,36]), since the evolution of stress with stretch (or strain) is similar to the evolution of the studied specimens. In addition, the main theory was derived from biomechanical considerations; thus, the model seems to have also explanatory power with reference to a possible mechanism of action of elastin and collagen fibers.

5. Conclusions

In this work, an explicit and deterministic approach was provided to investigate the stiffening phenomenon in axially stretched nerves. Although the complex nature of nerves was simplified using continuum, transverse isotropic, incompressible, elliptic cylinders, stretching experiments, involving porcine peroneal and canine vagus nerves, were closely reproduced. Moreover, since the presented approach was computationally light, explicit, and deterministic, it may be easily used as a further tool to investigate the variations of the nerve response due to any change in the content or nature of elastin and collagen. As a consequence, like previous implicit and probabilistic studies [22,23], it could be used to study pathologic states of peripheral nerves (e.g., due to diabetes). Furthermore, this framework, coupled to in vitro, in vivo, and human studies, could be used to better explore the effects of environmental cues on the nerve regeneration in peripheral nerve-gap lesions [37,38].