Biaxial Estimation of Biomechanical Constitutive Parameters of Passive Porcine Sclera Soft Tissue

This study assesses the modelling capabilities of four constitutive hyperelastic material models to fit the experimental data of the porcine sclera soft tissue. It further estimates the material parameters and discusses their applicability to a finite element model by examining the statistical dispersion measured through the standard deviation. Fifteen sclera tissues were harvested from porcine' slaughtered at an abattoir and were subjected to equi-biaxial testing. The results show that all the four material models yielded very good correlations at correlations above 96%. The polynomial (anisotropic) model gave the best correlation of 98%. However, the estimated material parameters varied widely from one test to another such that there would be need to normalise the test data to avoid long optimisation processes after applying the average material parameters to finite element models. However, for application of the estimated material parameters to finite element models, there would be need to consider normalising the test data to reduce the search region for the optimisation algorithms. Although the polynomial (anisotropic) model yielded the best correlation, it was found that the Choi-Vito had the least variation in the estimated material parameters, thereby making it an easier option for application of its material parameters to a finite element model and requiring minimum effort in the optimisation procedure. For the porcine sclera tissue, it was found that the anisotropy was more influenced by the fiber-related properties than the background material matrix-related properties.


Introduction
It is reported that majority of the vision impairment may be prevented or treated if caught early, and correct therapies are implemented. Globally, 86.1 million cases leading to the moderate and severe vision impairment were caused by undercorrected refractive error. Due to increase in aging population in many countries across the globe, vision impairment including uncorrected refractive error, glaucoma, and contract were reported to be in the rise. The rise in vision impairment is due to rising sociodemographic status and life expectancy witnessed in number of countries globally [1][2][3][4]. It is reported that 4.8 million and 21.4 million people are blind and visually impaired in 48 countries of sub-Saharan Africa [5]. In addition, Africa has been reported to be carrying approximately 19% of the global blindness [6]. Blindness prevalence rates vary widely but the evidence suggests that approximately 1% of people in sub-Saharan Africa are blind [7]. Visual impairment is a significant health problem in the sub-Saharan Africa. The major eye conditions include cataracts, uncorrected refractive errors, glaucoma, age-related macular degeneration, corneal opacities, diabetic retinopathy, trachoma, and onchocerciasis [8,9]. In addition, glaucoma is labelled as worldwide problem of the eye diseases where 111 million of individuals are projected to be affected by 2040 [10]. Therefore, it is vital that the biomechanical properties and computational models of the eye soft tissues are developed with the view to understand the mechanisms of various diseases of the eye. Biomechanical properties of eye soft tissues are vital in improving insight into how eye diseases progress [11].
Biomechanical properties of the soft tissues have been utilised previously to study mechanical behaviour [12][13][14]. The progression of the diseases was also studied by developing accurate and reliable computational models [15,16]. The detailed computational models of the eye request a thorough study of mechanical properties including the fitting of constitutive hyperelastic models. Previously, different biomechanical properties of the different ocular tissues including sclera [17][18][19], orbital connective tissue and fat [20,21], cornea [22][23][24], and extraocular muscles [25,26]. Sclera is normally associated with near-sightedness (myopia), and therefore, there is considerable amount of work done on the investigation of the biomechanical properties of sclera [27,28]. In particular, the mechanical properties of the sclera are critical for the development of synthetic materials for the replacement of eye tissues. This is because globally, there is a rise of accident where organ such as eye may be lost or damaged permanently. In our previous studies, we presented the biaxial mechanical characterisation of the sheep sclera soft tissue [28]. In this study, the Fung and Choi-Vito hyperelastic constitutive models were fitted on the biaxial tensile data. Finite element modelling (FEM) has been previously utilised in studying how soft tissues behave under mechanical loading [29][30][31][32][33][34][35]. These constitutive parameters generated by fitting biaxial tensile data may be utilised for detailed study of diseases' progression in the eye as well as understanding pain dissemination based on the needle injection in the sclera for treatment [36,37].
In this study, we present the biaxial biomechanical properties of the porcine sclera soft tissue. The porcine anatomy has been reported to be like that of human beings; therefore, it is believed that the biaxial mechanical data presented in this study will be useful for studying human eye diseases. Therefore, the aim of this study is to generate constitutive parameters of the porcine sclera soft tissues subjected to equi-biaxial loading. Fung, Choi-Vito, Holzapfel (2000), and polynomial (anisotropic) hyperelastic constitutive models are fitted in the equi-biaxial tensile data to estimate the constitutive material parameters. To our knowledge, there is no study that looked at utilising and comparing the Fung, Choi-Vito, Holzapfel (2000), and polynomial (anisotropic) hyperelastic constitutive models in estimating the material parameters of the porcine sclera soft tissue.

Methods and Materials
2.1. Sample Preparation. The sample consists of sclera (N = 15) soft tissue harvested from the local abattoir with unknown ages and sexes. To preserve the mechanical properties of the sclera soft tissue, the specimen was collected within 4 to 6 hours after slaughter as discussed [38]. The detailed tissue preparation of the porcine sclera is like what we have previously presented [28]. To avoid possible deterioration of the sclera tissue including dehydration, specimen was placed in ice bag on the way to the Unisa Biomechanics Laboratory. During further processing, the unwanted ocular tissues were removed from the porcine eye using surgical blades. Two square samples of 12 × 12 mm ( Figure 1) were cut from each eye roughly 4 to 6 mm away from the optic nerve head and corneal limbus to reduce thickness variation. Vernier caliper was then utilized to measure the thickness of Circumferential direction e sclera tissue is taken from Porcine eye Figure 1: Porcine sclera tissue samples (a) cut from each eye immediately after the delivery from the abattoir to the Unisa Biomechanics Laboratories with a sketch showing how the tissue was excised from the porcine eye (b). The reference axis were defined in such a way that the longitudinal direction is from the cornea to the base while circumferential direction is 90°with the longitudinal direction. 2 Applied Bionics and Biomechanics the sclera tissue. Randomly, the thickness of the sclera tissue was measured in four different places, and the average thickness was taken for further processing of the data. Specimens were subjected to visual inspection before individually subjected to biaxial testing equipment. Reason for visual inspection was to ensure that there was no physical damage on the tissue that may influence the mechanical data. Samples were dipped and kept in isotonic solution 9.0 g/l at pH of 5.5 for about 15 to 20 minutes before testing.

Mechanical
Testing. Biaxial mechanical properties of the porcine sclera soft tissue were captured using custom-built biaxial testing device, Bio-Tester 5000 from Cellscale (Waterloo, Canada) with load capacity of 23000 mN and accuracy of 5000 μN ( Figure 2). CellScale BioRakes were utilized to clamp the sclera soft tissue on all sides when equi-biaxial strain is applied. Strain rate of between 0.452 mm/s was applied during biaxial testing of each sample. To mimic the body temperature, the water bath was heated to 37°C. The temperature sensor is integrated in the custom-built biaxial device to ensure that the temperate is kept constant through testing. X and Y directions were considered as circumferential and longitudinal directions, respectively. Four (4) CellScale BioRakes were used to mount the specimens to develop mechanical properties. On both x and y directions, all fifteen (15) specimens were subjected to 25% strain using only strain-controlled mode equally in both directions.

Theoretical Framework
3.1. Tissue Stress-Strain Analysis. In this study, the stresses were calculated through the first Piola-Kirchoff stress T in the two directions using the equation: where F is the load vector in directions i = 1, 2, and for the current study, these indices represent longitudinal (crossfiber-direction) and circumferential (fiber-direction) direction. The index 0 denotes the undeformed state. L represents tissue length, and h represents tissue thickness.
The finite strains were calculated by the formula: The stress results in Equation (1) are however noisy; therefore, they were further filtered with an 8-point moving average filter in Excel. The data were resampled and further smoothened using a quadratic function.

Constitutive
Modelling. The passive response of a biological soft tissue is more complex than the response of elastic solids due to the fact that they undergo finite deformations under mechanical loads [39,40]. Thus, the passive phases of tissues' deformations are considered using the nonlinear theory of hyperelasticity. The most useful quantity in deriving the expressions for the passive response of materials that undergo finite deformation is the strain energy function [40]. There is a huge number of variations of its implementation in the literature depending on different cases. In this study, we assume that the tracheal tissue is anisotropic and incompressible; therefore, we apply and study the material parameters from six models, namely, the Fung, Choi-Vito, Holzapfel (2000), and polynomial (anisotropic) models. In the following sections, these models are presented in terms of their strain energy functions.

The Fung
Model. The Fung model is a hyperelastic anisotropic material model for stress-strain description of arterial wall. It is phenomenological in nature and its incompressible form is given by [41]. . The hooks (BioRakes) were utilised in securing the sclera tissue and were punched thoroughly before they were inserted into the warm bath of isotonic solution 9.0 g/l at pH of 5.5 at a temperature of 37°C.

Applied Bionics and Biomechanics
where E ij ði, j = θ, Z, RÞ are the Green strain components and ðθ, Z, RÞ are the cylindrical coordinates in the radial, circumferential, and axial reference directions, respectively. c stress-like material parameter and b i are the nondimensional material parameters.

Choi-Vito Model.
Choi-Vito model is hyperelastic anisotropic material model developed for canine pericardium. The model is fully phenomenological and formulated through the components of Green-Lagrange strain tensor. The model is implemented in an exponential format, and its strain-energy function is expressed as [42] where b i are the material parameters. (2000) Model. This model is hyperelastic anisotropic material model for stress-strain description of arterial layers. It is constituted by forms of the strain energy function that represent isotropy and anisotropy. Thus, its strain energy function is given by [43]. For incompressible formulation:

The Holzapfel
where the constant μ is associated with the noncollagenous matrix of the material, which describes the isotropic part of the overall response of the tissue. The constants c 1 and c 2 are associated with the anisotropic contribution of collagen to the overall response of the tissue. The material parameters are constants and do not depend on the geometry, opening angle, or fiber angle. The model is implemented in an exponential format, and I 1 , I 4 , and I 6 are (pseudo-) invariants of C (G 0 and H 0 ). G 0 , H 0 are structural tensors referenced to individual family of fibers.
3.2.4. The Polynomial (Anisotropic) Model. Like its name, this model is hyperelastic anisotropic material model whose strain energy function is expressed as a polynomial series of isotropic and anisotropic strain invariants given as [44] where a i , b j , are stress-like material parameters referenced to isotropic (matrix) properties. c k , and e m are stress-like material parameters referenced to anisotropic (fiber) properties. I 1 , I 2 , I 4 , and I 6 are (pseudo-)invariants of C (A 0 and B 0 ). A 0 , B 0 are structural tensors referenced to individual family of fibers. C is right Cauchy-Green definition tensor. The four material models are classified further according to the type of material parameters that they represent in Table 1. From the table, the models can be grouped into two main classes: those that model the soft tissue as layers of general material matrix and those that model the soft tissue as composite material with embedment of fiber structures within the layers of material matrix. These models may also be classified into purely phenomenological or structural scenario. Material parameter defining the orientation angle of fibers (measured from axis "1") in the undeformed configuration is included in the analysis. In this analysis, both the x and y directions were plotted simultaneously including the fibre angle term.

Data Analysis.
In this study, we use the constrained optimisation by linear approximation algorithm (COBYLA (3 rd party: SciPy)) implemented in Hyperfit software for fitting the equi-biaxial tensile experimental data of Fung, Choi-Vito, Holzapfel (2000), and polynomial (anisotropic) hyperelastic constitutive models given in Equations (3)-(6). The coefficient of determination, R 2 , which measures the correlation of one variable with another, is evaluated for each model. The coefficient of determination (R 2 ) (also known as Nash-Sutcliffe coefficient) is defined as follows [45][46][47]: Table 1: The classification of the four hyperelastic anisotropic material models according to the material properties they represent. The Fung and Choi-Vito are classified as the models that model the soft tissue as layers of general material matrix, while the polynomial (anisotropic) and Holzapfel (2000) models are those that model the soft tissue as composite material with embedment of fiber structures within the layers of material matrix.

Model
Applied Bionics and Biomechanics where y e is the experimental data, y m is the model predicted data, y e is the average value of the experimental data, the indices i, ⋯, n denote the data points, and R 2 ϵh−∞,1i, where a perfect fit is defined for R 2 = 1.

Results
Stress-strain curves for the porcine sclera soft tissue subjected to equi-biaxial tensile loading are plotted in Figure 3. The stress-strain curves for porcine sclera soft tissue show that the tissue has a well-defined nonlinear behaviour as the tissue gets stretched from the toe region to 25% strain. This study is aimed at estimating the biomechanical material parameters of a porcine sclera using four different constitutive material models, namely, Fung, Choi-Vito, Holzapfel (2000), and polynomial (anisotropic) models. All the above material models are hyperelastic anisotropic models which employ different modelling frameworks as presented in Equations (3)-(6) in Section 3.2. The material parameters for each model for all the fifteen tests are given in Tables 2-4. The tabulated results show that all the models yield very good correlations at greater than 96% with standard deviations of less than 3.2%. This implies that the model predictions were very close to the mean stress-strain curves for all the four models in this study.
However, the polynomial (anisotropic) model yielded the best correlation at an average of 98.3% and an average standard deviation of approximately 1.73%, while the Fung model yielded the lowest average correlation of 96.4% with an average standard deviation of 3.2%. The Choi-Vito and Holzapfel models had similar average correlations of 97.2% with average standard deviations of 1.65% and 2.98%, respectively.
In terms of statistical dispersion in the estimated material parameters by the four models, the polynomial (anisotropic) model had the largest variance with the widest scatter in the material parameters over the fifteen different tests. An evaluation of the standard deviation percentage yielded deviations from 110% to 660% for the polynomial (anisotropic) model. On the other hand, the Choi-Vito model had the smallest variance in the material parameters with deviations between 31% and 91%. This result might not be surprising considering that the Choi-Vito model was originally formulated for the pericardium tissue [42,48,49] which, in much respect, may be similar in functionality to the sclera tissue. The implications of such deviations in terms of applying the estimated average material parameters directly to finite element (FE) models is problematic, and this has been discussed in the next section. From Figure 4, the polynomial (anisotropic) constitutive model yields the highest R 2 of 0.983 while the Fung Table 2: Estimated material parameters of porcine sclera soft tissue loaded biaxially using Fung hyperelastic constitutive model including the coefficient of determination (R 2 ) of each specimen (N = 15). The average constitutive (c, b 1 , b 2 , b 3 , b 4 , b 5 , b 6 , and R 2 ) parameters are also included.    Figure 5 shows that the porcine sclera soft tissue was found to be much stiffer in the circumferential direction with stress values taken at 15% strain.

Discussion
The results of correlations between model predictions and tests presented in Section 4 should never be taken at face value. It should be noted that a model is applied to a given set of experimental data from which the material parameters are estimated. Thus, each material model is optimised for each fresh set of experimental data with totally different target values, and therefore, it is likely to yield material parameters that are different from those obtained from the other tests. This is where some form of calibration may become necessary to eliminate huge variances in the model parameters from one test to another [50][51][52][53]. Dubuis et al. [54] employed a calibration technique in which the FE itself was calibrated using 3D CT scans.
Therefore, in applying the average material parameters obtained here to an FE model, one may need to consider implementing an optimisation procedure in which these average values may be used as seed values [55]. Kollmann et al. [55] applied such methods for optimisation of the parameters used in FE modelling of asphalt mixtures. For models with large numbers of material parameters such as the polynomial (anisotropic) model in this study, it might be prudent to conduct an initial sensitivity analysis to determine the most influential material parameters. Biaxial tensile testing has been utilised in understanding the mechanical properties of soft tissues. Additionally, hyperelastic constitutive models have also been presented by fitting the constitutive models of the various soft tissues including heart myocardium [56][57][58] and omasum [59].
Despite the highest correlation obtained for the polynomial (anisotropic) model, the direct application of its parameters to an FE model can yield huge errors due to the wide scatter of its material parameters. One may need to prescribe S2 S3 S4 S5 S6 S7 S8 S9 S10 S11 S12 S13 S14 S15

Stress (kPa)
Circumferential direction 7 Applied Bionics and Biomechanics huge range of its parameter space within which the optimisation algorithm can search for the optimal material parameters to give accurate correlation for any given random test. Among the four material models in this study, the Choi-Vito model offers the best chance for attaining accurate correlations when directly applied to an FE model with minimal effort of optimisation. For the Choi-Vito model, only the stress-like material parameter may take more time in searching for optimal material parameter at 91% deviation from the mean value. Apart from the average value given in Table 3, the other best value from which to start the search for the stress-like material parameter for the Choi-Vito may be the gradient of the linear elastic portion of the stress-strain curve.
In terms of general material matrix anisotropy, both the Fung and Choi-Vito models do not show that the porcine sclera tissue is heavily anisotropic between the circumferential and longitudinal direction. However, the tissue seems to be different in material matrix in the radial direction from the Fung model. The values of a i and b j which also represent the material matrix properties in the polynomial (aniso-tropic) model show that there is no significant anisotropy between the circumferential and longitudinal directions. On the other hand, the c k material parameters in the polynomial (anisotropic) model and the k i parameters in the Holzapfel (2000) model, all of which represent material parameters related with fiber properties, exhibit significant anisotropy in Tables 5 and 4. This implies that the anisotropy in the porcine sclera tissue is mainly due to different fiber properties between the circumferential and longitudinal directions and not the background material matrix properties. It is difficult to establish in this study if the fiber orientation angle had any influence on the accuracy of the correlations in the polynomial (anisotropic) and Holzapfel (2000) models. This can be better investigated through the sensitivity analysis. Table 6 shows the box plot properties of the coefficient of determination (R 2 ) having Q1, median, and Q3. It is clear from the table that the four models have approximately equal coefficients of determination with the Fung and Holzapfel (2000) models exhibiting the largest and lowest ranges in the R 2 values as represented by the IQR values. Table 5: Estimated material parameters of porcine sclera soft tissue loaded biaxially using polynomial (anisotropic) hyperelastic constitutive model including the coefficient of determination (R 2 ) of each specimen (N = 15). The average constitutive (a 1 , a 2 , a 3 , b 1 , b 2 , b 3 , c 2 , c 3 , c 4 , c 5 , c 6 , φ, and R 2 ) parameters are also included.

Conclusion
The study applies four different hyperelastic material models to estimate the material parameters for the computational modelling of porcine sclera tissue. It is found that all the four material models produce very good correlations higher than 96% with the polynomial (anisotropic) model yielding the highest correlation of 98%. However, the variation of the material parameters from one test to another is a challenge as far as their direction application to FE models is concerned. For example, the polynomial (anisotropic) model has the minimum standard deviation of its material parameters from one test to another as 110%. Thus, the average material parameters obtained in this paper may result in inaccurate correlations in the FE predicted results. One way to deal with this challenge is to implement an optimisation routine which uses these material parameters as seed values. However, with such large spread in the parameter space, the search algorithm may take too long to converge especially with material models that have large numbers of estimated material parameters. Therefore, it is proposed that the test data should also be normalised or calibrated before application of the constitutive material models to reduce the variation in the estimated parameters thereby reducing the search region. In this study it is also found that the Choi-Vito model has the least variation in its material parameters. Thus, it might offer much better chance than the other models for convergence to material parameters that can yield very good correlations for the FE model. In terms of the models that are based on fiber-property formulation as identified in Table 1, the Holzapfel (2000) model offers a better option as compared to the polynomial (anisotropic) model for faster convergence as its material variation is much lower between the two models. As far as anisotropy in the porcine sclera tissue is concerned, both the polynomial (anisotropic) model and the Holzapfel (2000) model show that the fiber properties have more influence than the background material matrix. It is not clear if the angle of fiber orientation has any influence in the results. This might require more investigation.

Limitation
The limitation of this experimental study is that it did not segment porcine sclera per region. The thickness of porcine sclera samples used in this work was unknown. This may lead to unclear results since it was reported that some regions in the sclera exhibit anisotropy while others do not. However, it has been reported that porcine and human sclera behave in the same way despite a difference in tissue thickness. The two tissues are similar in terms of structure, histology, and collagen fiber architecture and should have similar material properties and will be useful for studying human eye diseases [60]. No imaging/microscopy was done in this study; hence, the fiber direction cannot be determined with certainty. It has previously been observed that the porcine sclera has a blood vessel that travels parallel to the ora serrata [61]. This vessel could introduce stiffening of one fiber direction and can result in different material properties compared to another fiber direction.
Further studies should be conducted to confirm human and porcine similarities with regard to material properties. Porcine sclera is easily available and can be an efficient substitute for human sclera material properties, in computational modelling.

Data Availability
All data is available in the main manuscript.

Additional Points
Correspondence and requests for materials should be addressed to F.N. Reprints and permissions information is available at http://www.nature.com/reprints.