# Residual strains estimation in the annulus fibrosus through digital image correlation

Gilles Dusfour, Dominique Ambard, Patrick Cañadas, Simon Lefloch

► To cite this version:

Gilles Dusfour, Dominique Ambard, Patrick Cañadas, Simon Lefloch. Residual strains estimation in the annulus fibrosus through digital image correlation. Journal of Theoretical, Computational and Applied Mechanics, 2021, pp.1-15. 10.46298/jtcam.6971 . hal-03058316v4

HAL Id: hal-03058316

https://hal.science/hal-03058316v4

Submitted on 25 Nov 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
**Identifiers**

DOI 10.46298/jtcam.6971

OAI hal-03058316v4

**History**

Received Dec 14, 2020

Accepted Sep 2, 2021

Published Nov 18, 2021

**Associate Editor**

Anna Pandolfi

**Reviewers**

Anonymous

Anonymous

**Open Review**

OAI hal-03406052

**Supplementary Material**

Data permalink
DOI 10.5281/zenodo.5658136

**Licence**

CC BY 4.0
©The Authors

# Residual strains estimation in the annulus fibrosus through digital image correlation

Gilles DUSFOUR, Dominique AMBARD, Patrick CAÑADAS, and Simon LE FLOC'H

LMGC, Université de Montpellier, CNRS, Montpellier, France

Back pain is now among the most prevalent health issues and generally originates from lesions in the annulus fibrosus. While the annulus fibrosus *ex vivo* mechanical properties are increasingly well understood, *in vivo* data are still missing. In particular, very few studies have precisely measured the residual strains within the annulus fibrosus and thus the *in vivo* deformation state of the annulus fibrosus is still misinterpreted and misevaluated. In this work, we propose an original and robust method for quantification of annulus fibrosus residual strains using digital image correlation techniques. Ten pig annuli fibrosi were extracted from adjacent vertebrae followed by a radial incision to release the residual strains. The operations were filmed and then analyzed with custom digital image correlation software in order to quantify the circumferential, radial, and shear residual deformations. Our results show that residual strains are of the same order of magnitude as those *in vivo*. The average circumferential strains are in tension on the outer periphery [3.32; 5.94]% and in compression on the inner periphery [-6.4; -1.69]%. The mean radial residual strains are essentially in compression (-10.4; 2.29]%). Locally, radial and circumferential residual strains can reach very high values of compression, up to 40 %. The mean shear strains remain very small (-0.04; 2.88]%). This study also shows that circumferential and radial residual strains evolve linearly along the radius and nonlinearly along the angle. We propose a simple model to predict their spatial variations. Our results and methods will allow the quantification of more realistic *in vivo* strains and stresses within the human intervertebral disc.

**Keywords:** annulus fibrosus; soft tissue; intervertebral disc; residual strain; mechanical characterization; growth; digital image correlation

## 1 Introduction

Back pains is the most prevalent health issue according to the work of Breivik et al. (2006) and might lead to cessation of professional activity and the emergence of psychological conditions like depression and anxiety (McWilliams et al. 2004). Back pain is mainly due to deteriorated or herniated intervertebral discs (IVDs) where the fribocartilaginous cells struggle to recreate annulus fibrosus (AF) tissue. However, numerical models of the IVD that take into account mechanical and nutritional factors highlight several possible scenarios of annulus fibrosus degeneration (Malandrino et al. 2011). While the *ex-vivo* mechanical behavior of the AF is well documented, very little is known about the *in-vivo* mechanical state. One of the more important factors in realistic models of the *in-vivo* state is accurate identification of the residual strain fields, as has already been done in arteries (Taber and Humphrey 2001).

The most commonly used technique to highlight residual strains is to make one or several incisions within the tissue and to observe its opening. Once the tissue is fully open, the residual strains are estimated through *opening angle method* (Fung 1983). Broad approximations are made as to the roundness of the studied tissue, which lead to global and sometimes poor estimations of the residual strain fields. Recently, Sigaeva et al. (2019) made improvements of the *opening angle method* in order to apply it to different parts of the same tissue. Such improvements allow several tissue geometries to be taken into account, including many cylindrical shapes, enabling the modeling of residual strains in asymmetric geometries. Estimations of residual strains and stresses
via radial incisions have been densely reported in the literature for cardiovascular tissues such as
arteries and heart (Omens et al. 2003; Labrosse et al. 2013), but little is known about the residual
strains in conjunctive tissues like the AF. To date, only Michalek et al. (2012), Duclos and Michalek
(2017) and Mengoni et al. (2017) have estimated the residual strains in the AF of bovine caudal
IVDs. These studies are based on the *opening angle method* and the *arc length method* to obtain
the circumferential strains, whereas the radial strains are estimated by following the evolution of
lamellae thickness during the stress relaxation experiment. However, the estimation of global
residual strains in these studies is based on the roundness of the tissue, which complicates
reproduction of this technique on the more complex geometry of human AF. Noteworthy,
Duclos and Michalek (2017) also used polarization birefringence microscopy to estimate the
fibre crimped period before and after the radial incision. This method shows circumferential
residual strains qualitatively equivalent to the arc length method but quantitatively different.
Therefore, we decided to realize a full field estimation of the residual strains through digital
image correlation (DIC) techniques and to apply this to the pig models more commonly used in
the IVD biomechanics community (Busscher et al. 2010; Yingling et al. 1999).

Residual strains are present in a large variety of biological tissues in both animal and plant
species (Goriely 2017). However, their genesis and purpose have not been clearly identified. Many
authors consider differential growth within organs to be responsible for residual strains (Goriely
2017; Zahn and Balzani 2018). Others hypothesize that gradients in mechanical properties or in
negatively charged molecules, such as proteoglycans, play a major role in the generation of
residual strains and stresses (Azeloglu et al. 2008; Yang and O'Connell 2018). According to Lanir
(2009), residual strains and stresses are the result of all of the above mechanisms combined.
Interestingly, residual strains seem to play a structural role in load-bearing tissues. Indeed, when
biological tissues undergoing internal pressure, such as arteries, heart, or AF, are cut in the radial
direction they tend to open up (Fung 1991; Michalek et al. 2012; Duclos and Michalek 2017). Such
residual strain or stress fields are assumed to create a homeostatic state of strain or stress in the
loaded organs. This hypothesis led to the development of strain-induced growth laws to explain
the pattern of residual strains in biological tissues (Taber and Humphrey 2001).

Several attempts have been made to model the development of biological tissues (Goriely 2017; Ambrosi and Mollica 2002; Lubarda and Hoger 2002). Interestingly, most of these models have been applied to cardiovascular tissues. To this date no growth model of the intervertebral disc has been proposed. Moreover, very few authors have compared experimental residual strain fields to those studied using simulations, excepting the work by Taber and Humphrey (2001). In order to completely model such a complex phenomenon, anatomical, nutritional, mechanical, and biological data are required and the results of the simulation finally need to be validated. Many of the incoming data are available in the literature through the work by Taylor (1975); Maroudas et al. (1975); Peacock (1951); Antoniou et al. (1996). Therefore, it is of crucial importance to accurately evaluate the residual strain fields to both validate and drive future IVD growth laws.

In this study, samples of pig AF were extracted from the adjacent vertebrae and cut in the
radial direction. We propose a new method to measure residual strains regardless of the AF shape.
The AF opening was recorded with a camera and analyzed through a dedicated DIC software to
evaluate the residual strain fields in two dimensions. Novel data on residual strain patterns
are reported. We hypothesized that the AFs display heterogeneous and anisotropic residual
strain patterns dependent on both radius and angle. Moreover, it is assumed that the residual
strain patterns can be modelled in terms of a linear relationship with radius and a nonlinear
relationship with angle.

2 Materials and methods

2.1 Specimen preparation

A total of ten AF sheet specimens were harvested from the spine of a three-month-old white pig weighing approximately 30 kg. Cervical (C2-C3, C3-C4, C5-C6, C6-C7), thoracic (T1-T2, T2-T3, T3-4, T10-T11,T11-T12), and lumbar (L1-L2, L6-S) AFs were extracted. The spine was obtained from a local slaughter house and immediately frozen at –20 °C in a physiological solution of NaCl (9 g/l). Prior to extracting the disc, the spinous and transverse processes were cut using a
saw. The discs were separated from the adjacent vertebral bodies with a homemade double blade cutter to obtain a planar AF sheet of 1 mm thickness. Following extraction, the AF was immediately cleaned with a physiological solution and graphite powder was deposited to enhance contrast. The samples were then immersed in a physiological solution of NaCl (9 g/l) at 37 °C in an aluminum container with a glass window on the bottom. This glass allowed images to be captured from below while making a radial incision with a scalpel on the upper lateral side of the AF, see Figure 1. The images were captured with a camera (IDS-UI-3360-CP-M-GL Rev2, Fujifilm

<figure><img src="image_2.png" /><figcaption>**Figure 1** Sheet specimen extraction protocol. Spinous and transverse processes are cut with a saw, the sheet of AF is obtained with a homemade double blade cutter, and the outer and inner specimens are extracted with a surgical knife. The DIC windows size is depicted in stage 4.</figcaption></figure>

lens Fujinon-1:1.4/16 mm) at an acquisition frequency of 1 Hz during the radial incision and 0.1 Hz for four hours during the residual strain relaxation period. The entire experimental setup is shown in Figure 2.

<figure><img src="image_3.png" /><figcaption>**Figure 2** Image acquisition setup used to acquire the images during the relaxation process.</figcaption></figure>

## 2.2 Numerical sample generation

For each experiment, a numerical geometry and a hexaedric mesh of the open sample were created with the GMSH mesh generator (Geuzaine and Remacle 2009). The boundaries of the numerical sample were defined manually to avoid any side effects. For each sample, the center point of the nucleus pulposus was determined with contour detection using the open cv2 python library. This point allowed us to compute an angle for each node of the mesh and to identify whether the node belonged to the anterior [0° – 60°], lateral [60° – 120°], or posterior [120° – 180°] quadrant as illustrated in Figure 3(b). This calculation can be extended for any division of the AF between two angles [α_min; α_max]. Such a division allows a refinement of the
analysis. Finally, a thermal calculation allowed us to compute a normalized radius as well as radial and circumferential direction for each node of the mesh, see Figure 3(a), Figure 3(c) and Figure 3(d).

## 2.3 Image analysis

### 2.3.1 Digital Image correlation

In order to estimate the residual strain fields, an iterative DIC technique was implemented in python 2.7. For each sample, a total of 30 images were analyzed with a mean number of 3480±2427 correlation points. For each point, a subset dimension of 85 px×85 px (4.6 mm×4.6 mm) was searched in a 150 px × 150 px (8.1 mm²) research area to identify both the displacement and the rotation of the subset. The research area in the deformed image was centered on the coordinates of the correlation point of the reference image. An initial angle interval of ±30° was chosen. For each correlation point, if the correlation criterion, in Equation (1), was higher than 0.9 the reference and deformed coordinates as well as the correlation criterion and the rotation angle were saved in a results array. Otherwise, another step of correlation was launched and an initial guess for the displacement and the rotation of the remaining points was computed using the mean displacement of the neighbouring points already correlated. If needed, the angle interval was enlarged during the computation to maximize the correlation criterion.

### 2.3.2 Correlation criterion

The zero-normalized cross-correlation criterion (ZNCC), as proposed by Pan et al. (2009b), was used to estimate the displacement of each point between the reference image f and the deformed image g,

$$r_{ij} = \frac{\sum_{m \in M} \sum_{n \in N} (f(m, n) - \bar{f})(g(m + i, n + j) - \bar{g})}{\sqrt{\sum_{m \in M} \sum_{n \in N} (f(m, n) - \bar{f})^2} \sqrt{\sum_{m \in M} \sum_{n \in N} (g(m + i, n + j) - \bar{g})^2}} \quad (1)$$

where M and N are the x and y intervals of the subset domain, respectively. The mean grayscale value of the reference and deformed image are represented by $\bar{f}$ and $\bar{g}$, respectively. The displacement in the x and y directions are determined by

$$(\Delta x, \Delta y) = \arg \max_{(ij)} \{r_{ij}\}. \quad (2)$$

### 2.3.3 Pointwise least square for strain estimation

In order to obtain full-field strain distribution from a noisy full-field displacement estimation, a pointwise local least-squares fitting technique was used as proposed by Wattrisse et al. (2001) and Pan et al. (2009a). The implementation of this technique requires, for each node of the mesh, the creation of a circular window (arbitrary radius of 40 px) and to extract all the correlation points in the neighbouring pixels. The displacement distribution in this window is then approximated as the linear plane

$$\begin{aligned} u(i, j) &= a_0 + a_1 x_i + a_2 y_j \\ v(i, j) &= b_0 + b_1 x_i + b_2 y_j \end{aligned} \quad (3)$$

where $u(i, j)$ and $v(i, j)$ are the original horizontal and vertical displacements at pixel location $(i, j)$. The least-squares method is then used to estimate the local unknown polynomial coefficients. Finally, the imposed displacement for each node of the hexaedric mesh is computed with identified polynomial coefficients and node coordinates. The strains are then computed at each Gauss Point element, through the derivative of the form function from the finite element technique, and extrapolated at each node location.

## 2.4 Radial and circumferential strain field estimation

For each node of the final closed mesh, the Eulerian strain tensor $\underline{e}$ is projected on both radial $\underline{e}_{RR}$ and circumferential $\underline{e}_{\theta\theta}$ directions (see Figure 3(c) and Figure 3(d)) to express the radial $\underline{e}_{RR}$,
circumferential $e_{\theta\theta}$ and shear $e_{\theta R}$ strains as

$$e_{\theta\theta} = \underline{e} : (\underline{e}_{\theta\theta} \otimes \underline{e}_{\theta\theta}), \qquad (4)$$

$$e_{RR} = \underline{e} : (\underline{e}_{RR} \otimes \underline{e}_{RR}), \qquad (5)$$

$$e_{\theta R} = \underline{e} : (\underline{e}_{\theta\theta} \otimes \underline{e}_{RR}). \qquad (6)$$

In order to facilitate the reading of residual strain fields, we created graphics to visualize the

<figure><img src="image_4.png" /><figcaption>**Figure 3** Numerical model parameters. The normalized radius, radial and circumferential vectors were estimated using a thermal flux, see Section 2.2. Anterior, lateral and posterior quadrants were determined by their specific angles from the antero-posterior central line.</figcaption></figure>

average evolution of residual strains according to the normalized radius and area of study. Therefore, the radius was divided into twelve equal intervals with a $\frac{1}{12}$ step. For each radius interval $[r_{min}; r_{max}]$ and each zone interval $[\alpha_{min}; \alpha_{max}]$, the median and interquartile range of residual strain were computed.

## 2.5 Residual strain model fitting

In order to take advantage of such rich residual strain fields, we derive a simple model to fully describe the residual strain state as a function of the local coordinates. The heterogeneous residual strain pattern of the AF was modelled by the following function linear in the normalized radius:

$$e_{\theta\theta}(r, \alpha) = K_{\theta\theta}(\alpha) \times r + O_{\theta\theta}(\alpha), \qquad (7)$$

$$e_{RR}(r, \alpha) = K_{RR}(\alpha) \times r + O_{RR}(\alpha), \qquad (8)$$

$$e_{\theta R}(r, \alpha) = K_{\theta R}(\alpha) \times r + O_{\theta R}(\alpha). \qquad (9)$$

Above, $K_{ii}$ is the gradient of residual strain along the normalized radius and $O_{ii}$, the offset of the residual strain for the zero normalized radius. In order to respect the strain continuity at 0° and 180°, the gradient $K_{\theta\theta}$ and the initial offset $O_{\theta\theta}$ of the circumferential strain were modelled as the
sum of cosine functions with respect to the normalized angle:

$$K_{\theta\theta}(\alpha) = k^K + \sum_{i=1}^{n} A_i^K \cos i\alpha, \qquad (10)$$

$$O_{\theta\theta}(\alpha) = k^O + \sum_{i=1}^{n} A_i^O \cos i\alpha. \qquad (11)$$

Fitting was performed using the `curve_fit` function from the `scipy/optimize` library of Python 2.7. The solution of the `curve_fit` function minimizes the squared error between the predicted and the actual data. The goodness of the fit between the predicted and measured residual strain e was evaluated with the normalized root mean square error (NRMSE):

$$\text{NRMSE} = \frac{\sqrt{\frac{\sum_{i=1}^{n} (e_i^{\text{num}} - e_i^{\text{exp}})^2}{n}}}{e_{\text{max}}^{\text{exp}} - e_{\text{min}}^{\text{exp}}} \times 100\%. \qquad (12)$$

## 2.6 Statistical analysis

The numerical data are presented as a median ± interquartile range (IQR). All statistical analyses were performed at 95 % confidence level using the statistical environment R. Nonparametric tests had to be applied because normality tests (Shapiro-Wilk) showed that none of the data groups has a normal distribution. The studied values are the radial strain e<sub>RR</sub>, the circumferential strain e<sub>θθ</sub>, the shear strain e<sub>θR</sub>, the gradient of each strain along the radius K<sub>ii</sub> and the initial offset O<sub>ii</sub>. We performed a Kruskal-Wallis test to assess the dependence of the radial strain e<sub>RR</sub>, the circumferential strain e<sub>θθ</sub> and the shear strain e<sub>θR</sub> on the normalized radius. Kruskal-Wallis tests were also performed on the gradient of residual strain K<sub>ii</sub> and on the initial offset O<sub>ii</sub> to assess their dependence on the normalized angle.

# 3 Results

## 3.1 Qualitative method results

After the incision, all the samples open up to release the residual strains. These results highlight circumferential tension in the outer AF and compression in the inner AF. Figure 4 shows residual strains fields for three different AFs. Qualitatively, these results suggest that the strain fields are widely heterogeneous and anisotropic. Indeed the circumferential residual Eulerian strains e<sub>θθ</sub> can exceed 15 % in both tension and compression while the radial residual strains e<sub>RR</sub> can reach 40 % in compression but rarely go into tension. Clearly, the radial e<sub>RR</sub>, circumferential e<sub>θθ</sub> and shear e<sub>θR</sub> residual strains present different pattern highlighting the anisotropic residual strain fields. Moreover, we note that both circumferential e<sub>θθ</sub> and radial e<sub>RR</sub> strains vary smoothly with radius and angle but also with disc level. Finally, the shape of the final AF presented in Figure 4 highlights the impaired ability of the DIC software to fully close the incisions.

## 3.2 Comparison to literature

In order to facilitate analysis and to compare our results with the literature (Duclos and Michalek 2017), the regional (anterior, lateral and posterior) variation along the normalized radius was averaged and is presented in Figure 5.

Figure 6 presents the quantitative data of the circumferential e<sub>θθ</sub> residual strains along the normalized radius for the anterior, lateral and posterior regions. The exterior periphery (r = 1) is in tension with a more pronounced tensile state on the anterior and posterior quadrant than in the lateral quadrant (respectively 5.94 % ± 1.54 %, 5.48 % ± 3.06 % and 3.32 % ± 6.65 %). These strains gradually decrease along the radius, getting closer to the zero strain state for the median radius (r = 0.5, respectively 2.26 % ± 4.94 %, -1.79 % ± 5.56 % and -0.12 % ± 5.25 %) and finally reach a global large compression state at the inner periphery (r = 0, respectively -6.4 % ± 13.71 %, -6.33 % ± 11.09 %, -1.69 % ± 9.24 %). We note the decreasing IQR with increasing radius. Importantly, the Kruskal-Wallis test shows that the circumferential residual strains e<sub>θθ</sub> are dependent on the normalized radius (p < 0.001).
<figure><img src="image_5.png" /><figcaption>**Figure 4** Circumferential, radial and shear residual strain fields for AFs.</figcaption></figure>

Quantitative data for the radial $e_{RR}$ and shear $e_{R\theta}$ residual strains along the normalized radius are presented in Figure 5(b) and Figure 5(c), respectively.

<figure><img src="image_6.png" /><figcaption>**Figure 5** Residual strains in the annulus fibrosus along the normalized radius for the three anterior, lateral and posterior quadrants. Circumferential (a), radial (b) and shear (c) strains. Black circles represent boxplot outliers.</figcaption></figure>

Like the circumferential strains $e_{\theta\theta}$, the Kruskal-Wallis test showed that the radial strains $e_{RR}$ are dependent on the normalized radius ($p < 0.05$) but not the shear residual strains $e_{R\theta}$. The radial residual strains $e_{RR}$ gradually increase in the anterior, lateral and posterior portions from the inner periphery ($r = 0$, respectively -10.4% ± 9.79%, -6.24% ± 10.75% and -4.99% ± 5.55%) to the outer periphery ($r = 1$, respectively -4.76% ± 9.07%, -5.77% ± 5.33% and 2.29% ± 4.97%). Finally, the shear residual strains $e_{R\theta}$ are globally close to the zero strain state (-0.04% ± 2.88%).

## 3.3 New local estimation of residual strain

### 3.3.1 Raw data

Nevertheless, the DIC method provides local estimates allowing us to increase the number of studied portions even with such a large amount of data. This technique allows us to identify significant differences between 15 different portions compartmentalized in 12° angular intervals, starting from the anterior side (0° – 12°) and ending at the posterior side (168° – 180°).

### 3.3.2 Linear model of residual strain

Figure 6(a), Figure 6(c) and Figure 6(e) display the average of all AF residual circumferential, radial and shear strains fields on an idealized pig intervertebral disc.
<figure><img src="image_7.png" /><figcaption>(a) Average circumferential experimental strains field</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(c) Average radial experimental strains field</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(e) Average shear experimental strains field</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(b) Circumferential strains</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(d) Radial strains</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(f) Shear strains</figcaption></figure>

**Figure 6** Projection of the average experimental strain fields on an idealized AF and plot of residual strains along the normalized radius for all zones of the AF. In (b), (d) and (f), average experimental data represented by solid lines with markers and linear model fits as dashed lines.

In addition, Figure 6(b), Figure 6(d) and Figure 6(f) show overlays of experimental and model-predicted residual strains $e_{\theta\theta}$, $e_{RR}$ and $e_{R\theta}$ versus the normalized radius $r$ for 15 different portions and 10 radius values of the idealized pig AF. We note that the IQRs are not presented for reasons of readability. For all the studied discs and zones, the predictions of the circumferential residual strains $e_{\theta\theta}$ agree extremely well with experimental data (NRMSE = 12.54% ± 6.58%) while the predictions of the radial $e_{RR}$ and shear $e_{R\theta}$ strains agree less well (respectively NRMSE = 20.9% ± 7.74% and NRMSE = 21.2% ± 7.31%).

The first order polynomials model allows us to obtain a gradient and an initial offset of residual strains along the normalized radius for 15 angular portions for all the studied discs. The gradients $K_{ii}$ and initial offsets $O_{ii}$ with respect to the angle $\alpha$ for the circumferential, radial and shear residual strains are presented in Figure 7.

While the circumferential gradient $K_{\theta\theta}$ and initial offset $O_{\theta\theta}$ evolve in a nonlinear way, both showing a peak centred on the [48° – 72°] interval, the radial and shear gradients, respectively $K_{RR}$ and $K_{R\theta}$, and corresponding initial offsets, $O_{RR}$ and $O_{R\theta}$, seem more stable with respect to the angle $\alpha$.
<figure><img src="image_13.png" /><figcaption>**Figure 7** Boxplots of the residual strain gradients and initial offsets along the normalized radius in the circumferential, radial and shear directions. The results of the cosine models are plotted in (a) and (b), showing good agreement with the experimental data. Circumferential residual strains: (a) gradient, (b) offset; Radial residual strains: (c) gradient, (d) offset; Shear residual strains: (e) gradient, (f) offset. For each angular value, the results of the cosine models for each sample are shown as gray circles. Error bars are calculated as follows: the lower value corresponds to the minimum value in the interval [Q1-1.5*IQR; Q1], the upper value corresponds to the maximum value in the interval [Q3; Q3+1.5*IQR]. IQR is the interquartile range (Q3-Q1). Q3 is the 75th percentile and Q1, the 25th percentile of the dataset.</figcaption></figure>

The results of the Kruskal-Wallis test show that the circumferential strains gradient $K_{\theta\theta}$ and the initial offset $O_{\theta\theta}$ are dependent on the angular portion (respectively $p < 0.05$ and $p < 0.005$). The Wilcoxon test with a Bonferroni correction highlight that the gradient of the [48° – 60°] and [60° – 72°] portions are significantly higher than the gradient of the [144° – 156°]; [156° – 168°] and [168° – 180°] portions ($p < 0.05$). On the other hand the Kruskal-Wallis test shows that the gradient and initial offset of the radial ($K_{RR}$, $O_{RR}$) and shear ($K_{R\theta}$, $O_{R\theta}$) residual strains are not dependent on the angular portion ($p > 0.05$).

Therefore, the evolution of the circumferential strains gradient $K_{\theta\theta}$ and initial offset $O_{\theta\theta}$ are modelled with a triple and double sum of cosine functions, respectively. The parameters of both fit functions are summarized in Table 1. Predictions of circumferential gradient and initial offset with respect to the angle agree well with the experimental data (respectively NRMSE = 14.19% and NRMSE = 12.85%). Predictions of circumferential gradient and initial offset are plotted in Figure 7(a) and Figure 7(b). Regarding the independence of radial and shear strains gradient and initial offset on the angle, all of the data have been approximated with a single scalar also
presented in Table 1.

<table><caption>**Table 1** Cosine function parameters to describe the evolution of residual strains with respect to the normalized radius and angle.</caption><thead><tr><th>Parameters</th><th>A<sub>1</sub></th><th>A<sub>2</sub></th><th>A<sub>3</sub></th><th>k</th><th>NRMSE</th></tr></thead><tbody><tr><td>K<sub>θθ</sub></td><td>0.045</td><td>-0.024</td><td>-0.033</td><td>0.096</td><td>14.19%</td></tr><tr><td>O<sub>θθ</sub></td><td>-0.031</td><td>0.051</td><td>-</td><td>-0.057</td><td>12.85%</td></tr><tr><td>K<sub>RR</sub> = 0.012</td><td>-</td><td>-</td><td>-</td><td>-</td><td>36.88%</td></tr><tr><td>O<sub>RR</sub> = -0.063</td><td>-</td><td>-</td><td>-</td><td>-</td><td>29.62%</td></tr><tr><td>K<sub>Rθ</sub> = 0.002</td><td>-</td><td>-</td><td>-</td><td>-</td><td>23.12%</td></tr><tr><td>O<sub>Rθ</sub> = -0.003</td><td>-</td><td>-</td><td>-</td><td>-</td><td>28.39%</td></tr></tbody></table>

Finally, Figure 8 shows the linear correlation between the residual strain gradients and the initial offset. Each point represents the residual strain gradient and the initial offset for a portion

<figure><img src="image_14.png" /><figcaption>**Figure 8** Gradient of residual strains $K_{ii}$ along the normalized radius with respect to the initial offset $O_{ii}$ for each portion of AF of each disc. Circumferential (a), radial (b) and shear (c) strains. These figures show the linear relationship between the strain gradient and the intensity of the strains. The more the strains tend towards a negative value on the inner periphery of the AF, the greater the strain gradient along the normalized radius.</figcaption></figure>

of a disc. These results establish that the greater the compression of the inner periphery the greater the gradient of residual strain along the normalized radius (NRMSE: $\theta\theta = 7.47\%$; NRMSE: $RR = 6.88\%$; NRMSE: $R\theta = 7.47\%$).

## 4 Discussions

The aim of this study is to develop an original method to locally evaluate residual strains of the pig AF using a digital image correlation technique allowing full field evaluation. In a similar spirit to Fung (1991) and Duclos and Michalek (2017) who respectively cut arteries and bovine AF in the radial direction to release the residual strains, we analyzed residual strains in pig AF and to the best of our knowledge are the first to use DIC to locally estimate the circumferential $e_{\theta\theta}$, radial $e_{RR}$ and shear $e_{R\theta}$ residual strains in this tissue.

The studies by Michalek et al. (2012) and Duclos and Michalek (2017), realized on bovine AF, reveal similar trends in the distribution of circumferential strains along the normalized radius except for the posterior region. Indeed, whereas the work by Duclos and Michalek (2017) presents a minimum of the circumferential strains in the middle zone of the posterior region, our study shows a steady increase from the inner to the outer zone of the AF, in all regions. We find that the amplitude of compressive circumferential strain on the inner zone is half that reported by Duclos and Michalek (2017). Moreover, our results reveal a generalized state of compression in the radial direction which is not in accordance with the state of tension in the inner and middle zones presented by Duclos and Michalek (2017). These large differences in both amplitude and state of deformation are probably due to the species of the extracted DIV.
While our study was realized on pig cervical, thoracic and lumbar AF, the studies by Michalek et al. (2012) and Duclos and Michalek (2017) used bovine tail IVD. Which lends well to this type of experiment because of its roundness and availability, although its morphology is different from the kidney shape of human and pig AF. However, despite showing similar morphology and mechanical behavior, human and pig IVD do not have equivalent dimensions and *in vivo* loading (Busscher et al. 2010; Beckstein et al. 2008). This highlights the necessity to evaluate the residual strain fields in human AF tissue. Therefore, our robust methodology for the estimation of residual strains is of crucial importance regarding the kidney shape of the human AF.

Furthermore, to the best of our knowledge this is the first time that shear residual strains of the AF are reported. Contrary to the circumferential and radial residual strains, the amplitude range of the shear residual strains is very small. The residual strain fields presented in this article, as well as the evolution of residual strains along the radius highlight the relevance of studying the residual strains via the material directions of the AF. Also, radial and tangential strains evolve in a smooth way along the normalized radius. On the contrary, this does not seem to be the case of shear strains. In addition, it has been shown that the residual strains evolve in a linear fashion along the normalized radius with good NRMSE [12% - 21%]. Our analysis highlights significant differences between the various angular portions studied, both in terms of gradient and initial offset of the residual strains. In the case of circumferential strains, the evolution of the gradient and the initial offset are significantly dependent on the angle. Two lateral portions ([48° - 60°] and [60° - 72°]) show gradients that are significantly greater than those of the posterior portions ([144° - 156°], [156° - 168°] and [168° - 180°]), the gradients of which are the lowest in AF. The gradients and the offsets of the radial and shear strains also evolve in a nonlinear way but more experimental data are needed in order to make robust conclusions about the differences between these portions. In a remarkable way, there are locations where the radial and circumferential strains are both positive or both negative, especially in the anterior area. This feature highlights the fact that AF can locally have a negative Poisson's ratio, as already shown by (Dusfour et al. 2020).

There have been many numerical studies in which the authors have attempted to predict the degeneration of the IVD by investigating the coupling between mechanical loading and biological properties (metabolic and nutrient transport, and growth and remodelling through mechanostransduction) (Huang and Gu 2008) and (Malandrino et al. 2011). While most of these studies used mechanical properties found experimentally, none of them considered the residual strain and stress fields as initial conditions of the numerical models. Our results reveal the large amplitude of the medial residual strain distribution along the normalized radius (circumferential: [-6.4%; 6%], radial: [-10.4%; 2.29%] and shear: [-0.04%; 2.88%]) which could greatly influence the computed global mechanical response as well as the nutrient transport. Such identified residual strain fields would reduce the hoop strains and stresses in the inner periphery of the loaded AF. As a consequence, circumferential strains are likely to be homogenized along the normalized radius during normal loading. On the other hand, the experimental studies in (O'Connell et al. 2007; O'Connell et al. 2011; Tsantrizos et al. 2005) found that the AF, regardless of the region, is in a radial tensile state when the entire IVD is loaded under compression, flexion or lateral bending. Moreover, the numerical study of (Yang and O'Connell 2018) shows large radial strains due to the swelling properties developed by the high glycosaminoglycan content. Therefore, the radial residual strains would tend to decrease the global radial strain state generated by both mechanical and chemical loads. In regards to all these features, it would be interesting to integrate an initial residual strain field in a numerical study by using the multiplicative decomposition of the deformation gradient, as proposed by (Jha et al. 2019; Peña et al. 2006; Rausch and Kuhl 2013; Balzani et al. 2006), and to compute its influence on the mechanical and nutrient transport fields. To this end, we propose a simple *cosinus* sum function model based on normalized radius and angle to fully describe the residual strain state within the IVD. Finally, in addition to the homogenization of the *in vivo* strains and stresses of the AF, the residual strains could also contribute to increase the global stiffness of the vertebral column.

The origins of such residual strain fields are still a matter of debate, and it is unclear whether these come from a gradient in glycosaminoglycan or a gradient in growth? While the circumferential strains that we see in the circumferential direction are equivalent to those induced
by the osmotic pressure in the study of (Yang and O'Connell 2018), our analyses of residual
strains in the radial direction lead to completely different results. These large discrepancies
highlight the fact that residual strains may not arise only from the osmotic pressure but also from
other phenomenon such as growth and remodelling. However, it seems difficult to segregate the
impact of osmotic pressure from the growth and remodelling process. Indeed, during the growth
and remodelling process both collagen proteins and glycosaminoglycan are secreted. While
the new created matter will exert forces on the old matter, the new glycosaminoglycan will
increase the negative charge density of the AF and thus increase the osmotic pressure. One
possible way to determine the impact of glycosaminoglycan distribution would be to perform
residual strain relaxation in hypertonic and hypotonic conditions as has already been done
on arteries by (Azeloglu et al. 2008). We have carried out some preliminary experiments by
modifying the NaCl bath from an isotonic saline (9 g/l) to a hypertonic saline (100 g/l). We
observed contraction by less of 2% in the radial direction and dilatation by up to 2% in the
circumferential direction. These preliminary results suggest that the effect of negative charges is
of minor importance in residual strain genesis but further studies are needed to fully understand
this phenomenon. It is now of crucial importance to develop growth models that take into
account the multiplicative decomposition of the gradient of the deformation but also the creation
of new negative charge (Armstrong et al. 2016). Finally, the strong linear correlation between the
level of residual strains at the inner periphery and the gradient of these same strains along the
normalized radius (Figure 8) highlights the spatial interdependence of the residual strains. Such a
phenomenon could be key point to the understanding of growth pattern and the development of
new IVD growth laws.

One major limitation of our work is the impossibility to determine residual strains in the axial direction. To this end, it may be possible to develop an experimental setup with a stereo-correlation technique. It could also be relevant to develop an adapted white and black powder to optimize the speckle pattern and thus greatly refine the digital image correlation results. The use of such a powder could also allow a more detailed study of the tensile test sample as well as the analysis of failure mechanisms of the AF. In addition, better experimental conditions would allow us to minimize the use of the pointwise least square technique to smooth the results. Indeed, a sensitivity study shows us that we probably underestimate the circumferential and radial strains in the lateral area because of our circular window of 40 pixels (5 % to 30 % of the radial thickness of the samples). The objective was to find a balance between softening the noisy experimental results and accounting for the large values of residual strains. Finally, as pointed out by Lee et al. (2020), most of the studies on residual strains focus on the continuous scale whereas the microscopic scale seems essential. Therefore, future studies should combine the relaxation of residual strains by the radial incision technique and the microscopic study of the *annulus fibrosus*.

5 Conclusion

By combining different methods, namely radial incision and digital image correlation, this study highlights the complex residual strain fields in the annulus fibrosus. In particular, we identify a strong and zone dependent circumferential residual strain gradient along the normalized radius, as well as stable and strong radial compressive strains. Our results show that our digital image correlation technique is a robust tool for the characterization of residual strains on the annulus fibrosus and highlights its potential for application to human samples. Our study shows that numerical models of the annulus fibrosus should integrate residual strains to better mimic the *in vivo* behavior of the herniated intervertebral discs. We hope that this work will lead to a better understanding of the origin of residual strains and allow a more detailed characterization. Future work could focus on assessing the impact of the anatomical level on the pattern and intensity of residual strains.

References

Ambrosi, D. and F. Mollica (2002). On the mechanics of a growing tumor. *International Journal of Engineering Science* 40(12):1297-1316. [DOI].

Antoniou, J., T. Steffen, F. Nelson, N. Winterbottom, A. Hollander, R. Poole, M. Aebi, and M. Alini
(1996). The human lumbar intervertebral disc: Evidence for changes in the biosynthesis and denaturation of the extracellular matrix with growth, maturation, ageing, and degeneration. *Journal of Clinical Investigation* 98(4):996-1003. [DOI], [OA].

Armstrong, M. H., A. B. Tepole, E. Kuhl, B. Simon, and J. Vande Geest (2016). A finite element model for mixed porohyperelasticity with transport, swelling, and growth. *PloS One* 11(4):e0152806. [DOI].

Azeloglu, E., M. Albro, V. Thimmappa, G. Ateshian, and K. Costa (2008). Heterogeneous transmural proteoglycan distribution provides a mechanism for regulating residual stresses in the aorta. *American Journal of Physiology-Heart and Circulatory Physiology* 294(3):H1197-H1205. [DOI].

Balzani, D., J. Schröder, and D. Gross (2006). Simulation of discontinuous damage incorporating residual stresses in circumferentially overstretched atherosclerotic arteries. *Acta Biomaterialia* 2(6):609-618. [DOI].

Beckstein, J., S. Sen, T. Schaer, E. Vresilovic, and D. Elliott (2008). Comparison of animal discs used in disc research to human lumbar disc: Axial compression mechanics and glycosaminoglycan content. *Spine* 33(6):E166-E173. [DOI].

Breivik, H., B. Collett, V. Ventafridda, R. Cohen, and D. Gallacher (2006). Survey of chronic pain in Europe: Prevalence, impact on daily life, and treatment. *European Journal of Pain* 10(4):287-333. [DOI].

Busscher, I., J. J. W. Ploegmakers, G. Verkerke, and A. Veldhuizen (2010). Comparative anatomical dimensions of the complete human and porcine spine. *European Spine Journal* 19(7):1104-1114. [DOI], [OA].

Duclos, S. and A. Michalek (2017). Residual strains in the intervertebral disc annulus fibrosus suggest complex tissue remodeling in response to in-vivo loading. *Journal of the Mechanical Behavior of Biomedical Materials* 68:232-238. [DOI].

Dusfour, G., S. LeFloc'h, P. Cañadas, and D. Ambard (2020). Heterogeneous mechanical hyperelastic behavior in the porcine annulus fibrosus explained by fiber orientation: An experimental and numerical approach. *Journal of the Mechanical Behavior of Biomedical Materials* 104:103672. [DOI], [HAL].

Fung, Y. C. (1983). What principle governs the stress distribution in living organs? *Biomechanics in China, Japan and USA* (Wuhan), pp 1-13.

Fung, Y. C. (1991). What are the residual stresses doing in our blood vessels? *Annals of Biomedical Engineering* 19(3):237-249. [DOI].

Geuzaine, C. and J.-F. Remacle (2009). Gmsh: A 3D finite element mesh generator with built-in pre-and post-processing facilities. *International Journal for Numerical Methods in Engineering* 79(11):1309-1331. [DOI], [HAL].

Goriely, A. (2017). *The Mathematics and Mechanics of Biological Growth*. Vol. 45. Interdisciplinary Applied Mathematics. Springer. [DOI].

Huang, C.-Y. and W. Y. Gu (2008). Effects of mechanical compression on metabolism and distribution of oxygen and lactate in intervertebral disc. *Journal of Biomechanics* 41(6):1184-1196. [DOI].

Jha, N.-K., J. Reinoso, H. Dehghani, and J. Merodio (2019). A computational model for fiber-reinforced composites: Hyperelastic constitutive formulation including residual stresses and damage. *Computational Mechanics* 63(5):931-948. [DOI], [OA].

Labrosse, M., E. Gerson, J. Veinot, and C. Beller (2013). Mechanical characterization of human aortas from pressurization testing and a paradigm shift for circumferential residual stress. *Journal of the Mechanical Behavior of Biomedical Materials* 17:44-55. [DOI].

Lanir, Y. (2009). Mechanisms of residual stress in soft tissues. *Journal of Biomechanical Engineering* 131(4):044506-1-5. [DOI].

Lee, T., M. Holland, J. Weickenmeier, A. Gosain, and A. B. Tepole (2020). The geometry of incompatibility in growing soft tissues: Theory and numerical characterization. *Journal of the Mechanics and Physics of Solids*, 104177. [DOI].

Lubarda, V. and A. Hoger (2002). On the mechanics of solids with a growing mass. *International Journal of Solids and Structures* 39(18):4627-4664. [DOI].

Malandrino, A., J. Noailly, and D. Lacroix (2011). The effect of sustained compression on oxygen metabolic transport in the intervertebral disc decreases with degenerative changes. *PLoS
*Computational Biology* 7(8):e1002112. [DOI], [OA].

Maroudas, A., R. Stockwell, A. Nachemson, and J. Urban (1975). Factors involved in the nutrition of the human lumbar intervertebral disc: Cellularity and diffusion of glucose in vitro. *Journal of Anatomy* 120(Pt 1):113. [OA].

McWilliams, L., R. Goodwin, and B. Cox (2004). Depression and anxiety associated with three pain conditions: Results from a nationally representative sample. *Pain* 111(1):77-83. [DOI].

Mengoni, M., O. Kayode, S. Sikora, F. Zapata-Cornelio, D. Gregory, and R. Wilcox (2017). Annulus fibrosus functional extrafibrillar and fibrous mechanical behaviour: Experimental and computational characterisation. *Royal Society Open Science* 4(8):170807. [DOI], [OA].

Michalek, A., M. Gardner-Morse, and J. Iatridis (2012). Large residual strains are present in the intervertebral disc annulus fibrosus in the unloaded state. *Journal of Biomechanics* 45(7):1227-1231. [DOI].

O'Connell, G., W. Johannessen, E. Vresilovic, and D. Elliott (2007). Human internal disc strains in axial compression measured noninvasively using magnetic resonance imaging. *Spine* 32(25):2860-2868. [DOI].

O'Connell, G., E. Vresilovic, and D. Elliott (2011). Human intervertebral disc internal strain in compression: The effect of disc region, loading position, and degeneration. *Journal of Orthopaedic Research* 29(4):547-555. [DOI], [OA].

Omens, J., A. McCulloch, and J. Criscione (2003). Complex distributions of residual stress and strain in the mouse left ventricle: Experimental and theoretical models. *Biomechanics and Modeling in Mechanobiology* 1(4):267-277. [DOI].

Pan, B., A. Asundi, H. Xie, and J. Gao (2009a). Digital image correlation using iterative least squares and pointwise least squares for displacement field and strain field measurements. *Optics and Lasers in Engineering* 47(7-8):865-874. [DOI].

Pan, B., K. Qian, H. Xie, and A. Asundi (2009b). Two-dimensional digital image correlation for in-plane displacement and strain measurement: a review. *Measurement Science and Technology* 20(6):062001. [DOI].

Peacock, A. (1951). Observations on the prenatal development of the intervertebral disc in man. *Journal of Anatomy* 257(6649):260-274. [OA].

Peña, E., M. A. Martinez, B. Calvo, and M. Doblaré (2006). On the numerical treatment of initial strains in biological soft tissues. *International Journal for Numerical Methods in Engineering* 68(8):836-860. [DOI].

Rausch, M. and E. Kuhl (2013). On the effect of prestrain and residual stress in thin biological membranes. *Journal of the Mechanics and Physics of Solids* 61(9):1955-1969. [DOI].

Sigaeva, T., M. Destrade, and E. Di Martino (2019). Multisector method for arteries: The residual stresses of circumferential rings with non-trivial openings. *Journal of the Royal Society Interface* 16(156):20190023. [DOI], [arXiv].

Taber, L. and J. Humphrey (2001). Stress-modulated growth, residual stress, and vascular heterogeneity. *Journal of Biomechanical Engineering* 123(6):528-535. [DOI].

Taylor, J. (1975). Growth of human intervertebral discs and vertebral bodies. *Journal of Anatomy* 120(Pt 1):49-68. [OA].

Tsanrizos, A., K. Ito, M. Aebi, and T. Steffen (2005). Internal strains in healthy and degenerated lumbar intervertebral discs. *Spine* 30(19):2129-2137. [DOI].

Wattrisse, B., A. Chrysochoos, J.-M. Muracciole, and M. Némoz-Gaillard (2001). Analysis of strain localization during tensile tests by digital image correlation. *Experimental Mechanics* 41(1):29-39. [DOI], [HAL].

Yang, B. and G. O'Connell (2018). GAG content, fiber stiffness, and fiber angle affect swelling-based residual stress in the intact annulus fibrosus. *Biomechanics and Modeling in Mechanobiology* 18(3):617-630. [DOI].

Yingling, V., J. Callaghan, and S. McGill (1999). The porcine cervical spine as a model of the human lumbar spine: An anatomical, geometric, and functional comparison. *Journal of Spinal Disorders* 12(5):415-423. [DOI].

Zahn, A. and D. Balzani (2018). A combined growth and remodeling framework for the approximation of residual stresses in arterial walls. *ZAMM-Journal of Applied Mathematics and Mechanics* 98(12):2072-2100. [DOI].
**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the authors—the copyright holder. To view a copy of this license, visit creativecommons.org/licenses/by/4.0.

<figure><img src="image_15.png" /></figure>

**Authors' contributions** G.D.: Investigation, Conceptualization, Formal analysis, Writing original draft, Software. D.A.: Methodology, Conceptualization. S.L.: Methodology, Conceptualization, Supervision, Writing, Review and Editing. P.C.: Project administration, Writing, Review and Editing.

**Supplementary Material** Supplementary material is available at the permalink 10.5281/zenodo.5658136. It corresponds to the experimental data and all the Python and LMGC software pipeline to process the data. The raw data correspond to images of the stress/strain relaxation process within the annulus fibrosus following a radial cut.

**Acknowledgements** The authors would like to thank Jonathan Bares for his help and advice on digital image correlation techniques.

**Ethics approval and consent to participate** Not applicable.

**Consent for publication** Not applicable.

**Competing interests** The authors declare that they have no competing interests.

**Journal's Note** JTCAM remains neutral with regard to the content of the publication and institutional affiliations.