Determining hyperelastic properties of the constituents of the mussel byssus system

Abstract

The mussel byssus system, comprising the adhesive plaque, distal thread, and proximal thread, plays a crucial role in the survival of marine mussels amongst ocean waves. Whilst recent research has explored the stress–strain behaviour of the distal thread and proximal thread through experimental approaches, little attention has been paid to the potential analytical or modelling methods within the current literature. In this work, analytical and finite element (FE) inverse methods were employed for the first time to identify the hyperelastic mechanical properties of both the plaque portion and the proximal thread. The results have demonstrated the feasibility of applied inverse methods in determining the mechanical properties of the constituents of the mussel byssus system, with the residual sum of squares of 0.0004 (N²) and 0.01 (mm²) for the proximal thread and the plaque portion, respectively. By leveraging mechanical and optical tests, this inverse methodology offers a simple and powerful means to anticipate the material properties for different portions of the mussel byssus system, thus providing insights into mimetic applications in engineering and material design.

---

1 Introduction

The mussel plaque within the mussel byssus system (Fig. 1b) can achieve a strong adhesive strength under wet conditions.^1–3 Therefore, scientists have conducted numerous studies on the chemical composition of mussel plaques,^4–6 inspiring the development of advanced adhesive materials in medical and engineering fields.¹ For example, six main mussel foot proteins (MFPs) were identified in mussel plaques, which play an essential role in maintaining adhesion at the plaque–substrate interface and establishing connections between threads and plaques.^4–6 These MFPs contain high levels of amino acid 3,4-dihydroxyphenylalanine (DOPA), which forms solid and reversible bonds with various surfaces and is an essential component of adhesive chemistry.⁴ Meanwhile, studies have also been carried out to investigate the mechanical behaviour of mussel plaques.^7–12 Desmond et al. (2015) analysed the adhesive performance of mussel plaques subjected to directional tension and summarised three different failure modes: adhesive, interior cohesive, and exterior failure.¹² Based on this fundamental research, Pang et al. (2023) further demonstrated that the failure modes of mussel plaques could be changed by the pull angle or the stiffness of the substrate.¹³

Fig. 1 (a) Mussel byssus system showing the threads and plaques, (b) a schematic illustration of a single mussel byssus structure and its protein composition, and (c) a schematic diagram of sample cutting positions: cutting position ① for collecting the whole mussel byssus system samples, ① and ③ for collecting thread samples, ② and ③ for collecting distal portion samples.

Furthermore, as researchers delved deeper into the mussel byssus system, the material properties of the constituents of the mussel threads were also revealed. The mussel thread comprises a distal and a proximal portion. Each portion has different compositions and it exhibits distinct mechanical behaviours that enable mussels to resist dynamic loads and effectively dissipate energy.^9,14 Specifically, preCol-D and preCol-P, two pre-pepsinized collagens, exhibit a prominent presence in the distal and proximal portions, respectively, potentially leading to different mechanical characteristics.^15–17 Gosline and Carrington (2004), revealed that the distal portion of M. californianus mussels can extend approximately half the length of the proximal portion, while demonstrating nearly twice the strength and fifty times the stiffness.¹⁰ Furthermore, when subjected to tension, the proximal portion does not exhibit a prominent point of failure, whereas the distal portion experiences four distinct stress regimes: elastic, plateau, hardening, and ultimate failure.^9,10 Additionally, researchers have also performed cyclic loading tests to analyse the self-healing mechanism of mussel threads.^7,18,19 Reversible sacrificial bonds between proteins and metals are widely recognised to play a crucial role in the self-healing process of mussel threads.^20–23

In order to understand the complex mechanics of the mussel byssus system, recent studies have established robust experimental methods to characterise its tensile properties, focusing primarily on the entire thread-plaque system, the mussel thread or solely on its distal portion.^10,24,25 Nevertheless, mechanical properties of proximal portion or plaque portion are yet to be solely studied due to geometrical limits and fragility when subjected to tensile loads.⁹ Notably, a fundamental nonlinear spring model was introduced by Waite et al. (2019) to obtain the spring constant and Youngs modulus for the plaque.²⁶ However, these material parameters were obtained by simplifying the stress–strain curve of the combined distal and plaque system into three linear regions.²⁶ Given the similarity between the nonlinear stress–strain behaviour exhibited by mussel threads and the hyperelastic relationship observed in rubber materials,²⁷ the errors in calculating the modulus of the distal and plaque portions by simplifying the overall response to linear are non-negligible. Furthermore, the study has not analysed stress–strain curves for the plaque and distal portion separately.²⁶ Therefore, further investigations of the individual contributions of the distal portion and the plaque to the overall stress–strain response is necessary.

In the present study, we first investigated the microstructure of three different portions within the mussel byssus system and established correlations between the microstructure and mechanical behaviour. Micro-tensile tests were then performed to obtain tensile force–displacement curves for the distal portion, the entire thread (proximal and distal portions), and the entire plaque-thread system. Meanwhile, the side view of plaque deformation was recorded to verify the inverse finite element (FE) approach. Finally, the analytical and FE inverse methods were applied to identify the material properties of the proximal and plaque portions. We validated the application of analytical and FE inverse methods in obtaining material properties for each of the three portions within the mussel byssus system and determined material constitutive models of those portions. The application of the FE inverse approach to obtain the potential plaque bottom surface detachment radius was also discussed in this research. These studies offer demonstrable evidence of the applicability of inverse methods in the analysis of multi-phase biomaterials and provide insight into the mechanical behaviour of each portion of the mussel byssus system for future bio-mimicking.

2 Materials and methods

2.1 Sample preparation procedure

2.1.1 Micro-tensile test sample preparation. Blue mussels (Mytilus edulis) were collected from the coast of Hunstanton in England and kept in a water tank filled with flowing and oxygen-rich seawater at temperatures between 8 °C and 10 °C. The collected mussels of 3–5 cm in length were tied to 50 × 50 mm transparent acrylic sheets to allow plaque deposition (see the ESI,† for additional details). Three different sample types were obtained for micro-tensile testing. The first type comprised whole mussel byssus systems, including an adhesive plaque, distal and proximal portions. These were collected by severing mussel threads close to the mussel body using dissecting scissors (Fig. 1c). The second type of sample is mussel threads, including both distal and proximal portions, separated from the adhesive plaques using blades (Fig. 1c). The third type exclusively encompassed the distal portions of mussel threads, collected using dissecting scissors (Fig. 1c).

2.1.2 Scanning electron microscopy (SEM). The SEM method was conducted to investigate the microstructures of the mussel byssus system. The threads and plaques were cut using dissecting scissors to prepare the SEM samples separately. The samples of the byssal threads and plaques were fixed in a 3.7% formaldehyde and 2.5% glutaraldehyde solution for approximately four hours and stored in Milli-Q water for three days at a temperature of 4 °C. The samples were first dehydrated in ethanol at increasing concentrations of 50% (2 times), 70% (2 times), 90% (2 times) and 100% (3 times) and then further dried in hexamethyldisilane (HDMS) 2 times for 5 min each.²⁸ After dehydration and drying, samples were sputter coated with 10 nm thick iridium for 120 s using a Quorum Q150V Plus and imaged using a JEOL 7100F FEG-SEM using an accelerating voltage of 5 kV.

2.2 Micro-tensile tests

Two types of micro-tensile tests: plaque tension (Fig. 2b) and thread tension (distal threads and whole threads) (Fig. 2c), were conducted on a customised micro-tensile experimental setup (Fig. 2a). For thread tests, mussel threads were held between two tweezers. The lower tweezer was replaced by an acrylic sheet to which mussel plaques were attached for plaque tests. The upper tweezer is connected to a load cell (Honeywell Model 34, 1000g-10 LBS) to record the loading force. The linear actuator (Thomson MLA11A05) connected to the load cell applies directional tension to the mussel thread at a slow loading speed of 0.25 mm s⁻¹ to measure the quasi-static tensile responses which are not sensitive to further reduction of the loading speed. The linear actuator is powered by a micro-stepping drive (Kollmorgen P5000) connected to the power supply. A monochrome camera (Thorlabs DCC1545M, 1280 × 1024 pixels, 25 fps), lens system (Navitar 12x Macro Zoom) and mirror were positioned on the acrylic plate side to record the side view of the plaque deformation process under tension.

Fig. 2 (a) Schematic illustration of a customised micro-tensile experimental setup, (b) focused plaque test setup, and (c) focused thread test setup.

2.3 Finite element (FE) forward modelling for mussel plaque

A three dimensional (3D) FE model of the mussel plaque was prepared using ABAQUS/Standard based on the actual measurement from the experiment (Fig. 3a and c). The diameter of the upper surface and the height of the mussel plaque were approximately 100 μm and 335 μm, respectively (Fig. 3a). Tensile loads perpendicular to the surface were applied to the upper surface of the plaque using displacement control. The plaque was modelled with 4-node 3D tetrahedral elements in ABAQUS notation (C3D4H elements). The boundary condition of the bottom surface for the plaque model was assumed first to be fully fixed. The Neo–Hookean model was applied to present the material properties of the mussel plaque. The more detailed equations of this material model can be seen in Section 3.1. The material properties of the mussel plaque model were iterated in the inverse optimisation algorithm until the tensile force–displacement curve and the deformed shapes derived from the FE model synchronised with the results obtained experimentally.

Fig. 3 (a) Side view of the mussel plaque, (b) side view of the FE model of the mussel plaque, (c) bottom view of the mussel plaque, and (d) bottom view of the FE model of the mussel plaque, all under undeformed configuration.

2.4 Fundamental of inverse approaches

It was assumed that within a composite material comprising a series of sub-materials, there exists an unidentified sub-material with unknown material parameters. Therefore, the inverse method can be applied to identify the material properties of the unknown sub-material by simplifying the composite as a series spring system (Fig. 4a). The total number of sub-materials was assumed to be p, where the material parameters of the sub-materials with number p − 1 are known (Fig. 4a). Therefore, three main inputs in the analytical inverse approach were as follows:

Fig. 4 (a) Schematic diagram of the simplified spring system, (b) schematic representation of the principle of the inverse method using displacement as a metric, (c) schematic representation of the principle of the inverse method using force as a metric, and (d) schematic of the deformed shape optimisation approach.

1. The experimental force–displacement data of the overall spring system (f^e,d^e), where f^e = [f₁^e, f₂^e, f₃^e,…, f_k^e]^T, d^e = [d₁^e, d₂^e, d₃^e,…, d_k^e]^T, and k is the number of data points obtained from the force–displacement curves (Fig. 4b and c).

2. Material parameters of known sub-materials ω_i, where i = 1, 2, 3,…, p − 1.

3. Guessing material parameters of the unknown sub-material m = [m₁, m₂, m₃,…, m_q]^T, where q is the number of material parameters for the unknown sub-material.

As the spring system is compounded in series, the tensile force of the overall spring system should be equal to the tensile force of any sub-material, and the tensile displacement of the overall spring system should equal to the sum of the elongations of all sub-materials:

f^e = f^1e = f^2e = ⋯ = f^pe

d^e = d^1e + d^2e +⋯+ d^pe (1)

Because the material parameters ω_i, material constitutive equations and sizes of the p − 1 sub-materials were assumed to be known, the elongations of these sub-materials can be calculated based on eqn (1). Meanwhile, the elongation of the unknown sub-material d^u(m) can also be calculated using the guessing material parameters m, assumed material constitutive equation, measured sizes and the force of unknown sub-materials f^u(m) = f^e. Notably, d^u(m) and f^u(m) denote their values are functions of m, and the rule applies to subsequent arguments containing m in parentheses.

Therefore, the displacement of the overall spring system from the analytical inverse approach d(m) will be:

d(m) = d¹ + d² +…+ d^p−1 + d^u(m) (2)

where d¹,d²,…,d^p−1 are the elongation vectors for known sub-materials calculated from the known material parameters ω_i, material constitutive models, measured sizes and experimental forces.

Thus, the residual vector F(m) based on the guessing material parameters m can be defined as (Fig. 4b):

F(m) = d^e − d(m) = [F₁(m), F₂(m), F₃(m),…, F_k(m)]^T (3)

If the force is used as the optimisation metrics, the residual vector F(m) can be determined as (Fig. 4c):

F(m) = f^e − f(m) (4)

where f(m) is the force vector of the overall spring system based on the guessing material parameters m from the analytical inverse approach.

The optimisation process is achieved by applying the nonlinear least squares estimation function lsqnonlin within the Matlab Optimisation Toolbox.²⁹ This method effectively addresses curve-fitting problems characterised by using the following form of nonlinear least-squares:²⁹

The lower and upper bounds lb and ub were set to limit the value ranges of material parameters m and thus accelerate the inverse approach. The functional tolerance θ is initially set to 0.01 in the ‘options’ and can be changed to the desired degree of precision.²⁹ When the final change in the sum of squares relative to its initial value is less than the value of the function tolerance, lsqnonlin will stop and the optimised material parameters m* will be obtained (Fig. 5):

Fig. 5 Flow chart of both analytical inverse and FE-based inverse procedures.

However, if the shape of the material is irregular, the analytical inverse approach is no longer applicable, and the FE-based inverse approach is needed. The “.inp” file for ABAQUS was first written, including the material geometry, known and guessed material parameters, material constitutive models and boundary conditions. A post-processing Python script was then created to extract force–displacement results based on guessing material parameters from the ABAQUS “.odb” file. The remaining optimisation process is consistent with the analytical inverse method by applying lsqnonlin, meaning that when eqn (6) is satisfied, the optimised material parameters m* will be obtained (Fig. 5).

3 Mussel byssus system material constitutive models and inverse procedure

3.1 Hyperelastic models

Based on the characteristics of the stress–strain curves of the distal and proximal threads and the force–displacement curves of the entire byssus system in the previous literature,^9,10,12,30 different hyperelastic material constitutive models were assumed to present plaques, distal threads and proximal threads, respectively. The incompressible Ogdens model was applied to the distal threads because the shape of the strain–stress curves for the distal threads has characteristics similar to those of the Ogden's model curve. Similarly, the constitutive model of both the proximal portion and the plaque was assumed to be the incompressible Neo–Hookean material model due to the resemblance of the curve shape.

The strain energy potential W of the three-order Ogden model is expressed in terms of principal stretches as follows³¹:

where N is the order of strain energy potential (N = 3), μ_β and α_β (β = 1, 2, 3) are six material constants. [small lambda, Greek, macron]₁, [small lambda, Greek, macron]₂ and [small lambda, Greek, macron]₃ are the distortional principal stretches which can be calculated using principal stretches , J is the total volume ratio.

Meanwhile, the strain energy potential of the incompressible Neo–Hookean material model is given below:³²

W = C₁([small lambda, Greek, macron]₁² + [small lambda, Greek, macron]₂² + [small lambda, Greek, macron]₃² − 3) (8)

where C₁ is a material constant.

Tensile tests were assumed to be uniaxial. Therefore, the uniaxial stress σ and stretch λ for the Ogdens model can be expressed as

The uniaxial stress σ and stretch λ for the Neo–Hookean hyperelastic model are defined as shown below:

3.2 Application of inverse approaches to the mussel byssus system

3.2.1 Problem definition. The shorter length of the proximal portion compared to the distal portion of threads and the irregular shape of mussel plaques increase the difficulty of the micro-tensile tests. Therefore, analytical and FE inverse approaches were employed to obtain their material properties.

3.2.2 Analytical inverse approach for the proximal thread. The three portions of the mussel byssus system were simplified as a spring system in series. Tensile force and displacement data for the distal portion (f_d^e,d_d^e) was obtained from the micro-tensile tests for the distal portion, and the stress–strain data (σ_d,ε_d) can be further calculated by
where r_d and l_d are the radius and length of the distal portion before tension.

Based on the experimental data and the hyperelastic Ogden model (eqn (9)), six material parameters of the distal portion (μ_β, α_β), (β = 1, 2, 3) can be determined using a simple analytical or FE-based inverse approach shown in Fig. 5 and were used as input for the further analytical inverse approach.

Next, the tensile force and displacement data of the whole thread (distal and proximal portions) (f_t^e,d_t^e) was also obtained from micro-tensile tests for the whole thread. Based on eqn (1), the following equations can be obtained:

f_t^e = f_p^e = f_d^e

d_t^e = d_p^e + d_d^e (12)

where f_p^e = [(f_p^e)₁, (f_p^e)₂,…,(f_p^e)_k]^T and f_t^e = [(f_t^e)₁, (f_t^e)₂,…,(f_t^e)_k]^T are the tensile force vectors of the proximal portion and the distal portion, respectively. d_p^e = [(d_p^e)₁,(d_p^e)₂,…,(d_p^e)_k]^T and d_t^e = [(d_t^e)₁,(d_t^e)₂,…,(d_t^e)_k]^T are the elongation vectors of the proximal portion and the distal portion, respectively.

Therefore, the elongations of the distal portion can be calculated based on eqn (9), (11) and (12), and the six material parameters obtained from the previous distal portion tensile tests:

where j = 1, 2, 3,…,k

Meanwhile, the elongation of the proximal portion d(C₁^p)_p based on the given material constant C₁^p can be calculated in the analytical inverse approach based on eqn (10) and (12):

where r_p and l_p are the radius and length of the proximal portion before tension. (d(C₁^p)_p)_j is one data point within d(C₁^p)_p and a function of C₁^p. (The rule also applies to subsequent arguments containing C₁^p in parentheses)

As a result, the tensile displacement of the whole thread based on the given material constant C₁^p can be expressed as

d(C₁^p)_t = d(C₁^p)_p + d_d^e (15)

The residual vector F(C₁^p) based on the given material constant C₁^p then can be defined as

F(C₁^p) = d_t^e − d(C₁^p)_t = [F₁(C₁^p), F₂(C₁^p), F₃(C₁^p),…, F_k(C₁^p)]^T (16)

If the tensile force was used as metrics, the residual vector F(C₁^p) can be expressed as

F(C₁^p) = f_t^e − f(C₁^p)_t (17)

lsqnonlin was applied for parameter optimisation in section 2.4, and the functional tolerance θ was initially set as 0.01. lsqnonlin will stop when the following convergence criterion is satisfied:

3.2.3 FE-based inverse approach for the mussel plaque. The FE inverse approach was employed as the irregular shape of mussel plaques makes the above analytical inverse approach inapplicable. Micro-tensile tests were first performed on the combined distal and plaque portions because it was difficult to perform experiments on the plaques themselves. Based on eqn (1), the force of the mussel plaque f_pl^e can be obtained from the experiment:

f_pl^e = f_d^e = f_d+pl^e (19)

where f_d^e and f_d+pl^e are the tensile force vectors of the distal portion and distal-plaque system, respectively.

Meanwhile, the tensile displacement of the upper surface of the mussel plaque d_pl^e was also measured from the side-view video (Fig. 11a) to generate the experimental force–displacement data (f_pl^e,d_pl^e) as input in the FE inverse approach. The guessing material parameter C₁^pl of the mussel plaque was then written on the ABAQUS CAE Interface and updated in lsqnonlin. A post-processing Python script was created to extract the force–displacement data (f(C₁^pl)_pl,d(C₁^pl)_pl) based on the guessing material paramter C₁^pl from the ABAQUS odb file.

The residual vector F(C₁^pl) and the sum of squared residuals from the force–displacement data based on the guessing material paramter C₁^pl can be defined as:

F(C₁^pl) = d_pl^e − d(C₁^pl)_pl = [F₁(C₁^pl),F₂(C₁^pl),…,F_k(C₁^pl)]^T

To further optimise the FE inverse approach, the deformed side views of the mussel plaque measured from the experiment were compared with the FE results. It was assumed that the deformed shapes with a total number of k₁ were analysed. For each deformed shape, the leftmost node was set as (0,0) (Fig. 3a), and nodes with a total number of k₂ were selected on each shape. Graph Grabber v 2.0.1 software was used to obtain the coordinates (y,z) of selected nodes (https://www.quintessa.org/software). Meanwhile, the coordinates (y,z(C₁^pl)) of corresponding nodes based on the guessing material parameter C₁^pl were exported from the FE model using a post-processing Python script.

The residual vector F^s(C₁^pl) of each deformed shape can be expressed as (Fig. 4d)

F^s(C₁^pl) = z − z(C₁^pl) = [F₁^s(C₁^pl),F₂^s(C₁^pl),…,F_k2^s(C₁^pl)]^T (21)

Thus, the sum of squared residuals for all the deformed shapes were defined as:

The total sum of squared residuals R(C₁^pl) for both force–displacement curves and deformed shapes is defined below:

where W₁ and W₂ are the non-negative weights, W₁,W₂ ∈ [0,1], W₁ + W₂ = 1. W₁ and W₂ were set to 0.8 and 0.2, respectively, indicating that the accuracy of the force–displacement data is more critical in this research. Their values can be adjusted according to study interest.^33,34

lsqnonlin was applied for parameter optimisation in section 2.4, and the functional tolerance θ was initially set as 0.01. The value of θ can be adjusted based on the optimisation conditions. lsqnonlin will stop when the following convergence criteria are satisfied:

The analytical inverse and FE-based inverse procedures specifically for the proximal portion and the plaque portion are summarised in Fig. 6.

Fig. 6 Flow chart of determining the hyperelastic properties of the plaque and the proximal portion by using FE-based inverse approach and analytical inverse approach, respecticvely.

4 Results

4.1 Microstructures of the mussel byssus system

Four different microstructures from the mussel byssus system were identified by SEM analysis (Fig. 7c–f). For the proximal portion, it can be seen that numerous folds of various sizes cover almost the entire surface of the proximal portion, and the furrows between the folds are parallel or perpendicular to the longitudinal axis (Fig. 7c). The surface of the distal portion is smooth, but longitudinal patterns of approximately 1 μm in width can still be seen on the surface (Fig. 7e). Meanwhile, the surface of the transition part between these two portions has the characteristics of both, with a smoother surface and a reduced number of folds compared to the proximal portion, and longitudinal patterns begin to appear (Fig. 7d). Due to the difficulty in identifying the location of the transition portion during the experiments, it was assumed that mussel threads were mainly composed of proximal and distal portions in this research. Furthermore, the mussel plaque with a porous core structure that is distinctly different from that of the thread is shown in Fig. 7f. Different microstructural characteristics and protein composition^35–37 may be the reasons for the difference in material properties for different portions and provide evidence for the feasibility of the spring system assumption, which is the main basis of inverse approaches (Fig. 7a and b).

Fig. 7 (a) The idealization of the thread-plaque mechanics model, (b) the overall thread–plaque structure, and (c)–(f) the micro-structures of the proximal portion, transition portion, distal portion and plaque core, respectively.

4.2 Micro-tensile test results

It is evident from Fig. 8a and b that the force–displacement and strain–stress curves obtained from the tensile tests in the distal portions exhibit different failure points, reflecting the unique characteristics of each thread.^24,38,39 Nevertheless, the distal portion tests show consistent failure mechanism patterns, suggesting the applicability of a single material constitutive model to represent the distal portion.

Fig. 8 (a) The tested distal portion force–displacement curves, (b) the converted nominal stress–strain curves, (c) the tested whole thread (proximal and distal) portion force–displacement curves, and (d) the tested distal–plaque portion force–displacement curves.

Meanwhile, the force–displacement curves of the combined distal and proximal portions and the combined distal and plaque portions are shown in Fig. 8c and d, respectively. An inflection point can be found in some of the force–displacement curves of the combined distal and plaque portions (Fig. 8c), which may indicate that the ultimate failure of the combined portions is likely to occur after the yield point of the distal portion. In contrast, the force–displacement curves of the combined distal and plaque portions did not show an obvious inflection point, suggesting that plaque detachment may occur before the yield point of the distal section. The comparison of Fig. 8c and d also demonstrates that the failure displacement of the combined distal and proximal portions is larger than the displacement of the combined distal and plaque portions, and the failure force varies from 0.1 N–0.2 N during these tests.

4.3 Material properties of individual constituents

4.3.1 Material properties of the distal portion. The “fminsearch” method in MATLAB was first used to generate a target curve that closely matched the experimental data, which was then used as input to the inverse approach to speed up the inverse calculation. The optimised curves of the analytical and FE inverse approaches with the assumed hyperelastic Odgen model show close agreement with the target curve (Fig. 9). The material properties obtained from both inverse approaches are shown in Table 1 and were used as input for the latter proximal and plaque inverse approaches.

Fig. 9 Nominal stress–strain curves of the distal portion obtained by experiment, fitting and inverse.

Table 1 Material properties of the distal portion

Material constants	Analytical inverse solution	FE inverse solution
μ 1 (MPa)	−529	−528
α 1	7.19	6.69
μ 2 (MPa)	425	426
α 2	7.12	6.79
μ 3 (MPa)	229	229
α 3	−14.5	−12.4

---

4.3.2 Material properties of the proximal portion. The process of optimising the material parameters for the proximal portion followed the method shown in Fig. 6. Different thread dimensions can cause a significant variation in the force–displacement curve. Therefore, instead of using the average experimental value, test 2 was selected from the three thread (proximal and distal) portion force–displacement curves in Fig. 8a to remove the noise and set as the target curve. The force–displacement curves from inverse approaches on the basis of the assumed Neo–Hookean model are in good agreement with the experimental result (Fig. 10). Therefore, the material parameter of the proximal portion C₁^p and residuals from both inverse approaches are given in Table 2.

Fig. 10 Nominal force–displacement curves of the thread (proximal and distal) portion obtained by experiment, fitting and inverse.

Table 2 Proximal portion material parameter C₁^p(MPa) and the min of the sum of square residuals (N²)

Inverse methods	C 1 ^p(MPa)
Analytical inverse method	9.16	0.0004
FE-based inverse method	11.2	0.0005

---

4.3.3 Material properties of the plaque portion.
4.3.3.1 Material parameter C₁^pl of the plaque portion. Due to the complex and challenging nature of conducting tensile tests on the small and fragile mussel plaques, the distal–plaque tensile test results were analysed instead of conducting plaque-only experiments. From the experimental video, the displacement of the plaque portion prior to failure was found to be 0.25 mm. Subsequently, by applying eqn (19) to the force data obtained from the distal–plaque test 1 (Fig. 8d), the force–displacement curve for the individual plaque can be constructed, as shown in Fig. 11b and d. Furthermore, three deformation shapes at different displacements were obtained from test 1 (Fig. 11a). The leftmost points of the deformed shapes were set as the coordinate origins, and the shapes were further converted into three shape curves using Graph Grabber v 2.0.1 (Fig. 11c and e).

Fig. 11 (a) The side view of plaque deformation at d_pl = 0 mm, d_pl = 0.16 mm, and d_pl = 0.25 mm. (b) The force–displacement curves without detachment at d_pl = 0 mm, 0.16 mm, 0.25 mm, (c) the deformed shape comparison of the plaque portion without detachment at d_pl = 0 mm, 0.16 mm, 0.25 mm, (d) the force–displacement curves with detachment at d_pl = 0 mm, 0.16 mm, 0.25 mm, and (e) the deformed shape comparison of the plaque portion with detachment at d_pl = 0 mm, 0.16 mm, 0.25 mm.

The experimental force–displacement and shape curves (Fig. 11b–e) were used as input for the FE inverse approach, then the material parameter for the plaque portion C₁^pl was output as 1.2 MPa with a minimum sum of squared residuals 0.137 mm² (Table 3). The general agreement between the force–displacement curves and the shape curves obtained from the experimental results and the inverse results can be seen in Fig. 11b and c. However, by carefully comparing the shape curves of the FE and the experiment, it can be observed that the shapes of the FE model with a fully fixed bottom surface showed obvious deformation only in the upper part (Fig. 11c). On the other hand, the shapes of the experiment showed a more monolithic deformation, which also occurred in the area of the plaque that is far away from the stretched thread (Fig. 11c). Meanwhile, after d_pl = 0.16 mm, the experimental curves were slightly lower than the simulated curves (Fig. 11b). Therefore, the results could indicate a possible detachment on the bottom surface at d_pl = 0.16 mm during the experiment.

Table 3 Plaque portion material parameter C₁^pl(MPa) and the minimum value of the sum of square residuals (mm²) from the FE inverse approach

	Inverse output without detachment
C 1 ^pl	1.2
min[R(C1^pl)]	0.137

---

4.3.3.2 Potential bottom detachment radius optimisation. In order to further analyse the bottom surface detachment radius at d_pl = 0.16 mm, it was assumed that the detachment area was a circle and that the projection of the centre point on the upper surface to the lower surface was used as the detachment centre point. Furthermore, the material parameter C₁^pl was assumed to be the same as 1.2 MPa and the detachment radius from d_pl = 0.16 mm to d_pl = 0.25 mm was assumed to be constant. The updated target force–displacement curve shown in Fig. 11d was applied as the new input. Meanwhile, the radius of bottom surface detachment was optimised inside the FE inverse approach instead of the material parameter C₁^pl. The 3D FE model used within this FE inverse approach was the same as in section 2.3. The only difference was that the fixed boundary conditions at the bottom points within the detachment area were removed at d_pl = 0.16 mm to simulate the bottom surface detachment during the optimisation process. Specifically, the detachment radius had to be less than the distance from the centre of the detachment to the outer boundary of the bottom surface.

The updated inverse results with detachment can be seen in Fig. 11d and e with a detachment radius = 0.6 mm and a less minimal sum of squared residuals 0.01 mm² (Table 4). The FE and experimentally deformed shapes at d_pl = 0.16 mm and 0.25 mm were more consistent than the previous results (Fig. 11d and e). The discrepancy between the FE and the experimentally deformed shapes at d_pl = 0.25 mm suggests a possible detachment range greater than 0.6 mm. However, it is important to note that 0.6 mm represents the maximum detachment radius without touching the edge of the bottom surface. This observation suggests that the detachment might be a dynamic process in which the detachment range gradually expands and is not static. Overall, the implementation of shape optimisation in the mussel plaque's FE inverse approach leads to a more accurate result.

Table 4 Bottom surface detachment radius r_b(mm) and the minimum value of the sum of square residuals (mm²) from the FE inverse approach

	FE inverse approach output
r b	0.6
min[R(C1^pl)]	0.01

---

5 Discussion

In recent decades, several reports have investigated the tensile stress–strain behaviour of the mussel byssus system as a whole and its distal and proximal portions individually.^9,10,12,30 However, due to the size of the mussel byssus system and the complexity of the biomaterial composition, there are still some challenges to comprehensively understand its mechanism using experimental or numerical methods. In our study, the analytical and FE inverse approaches were applied for the first time to obtain the hyperelastic properties of the constituents of the mussel byssus system in this investigation, which could simplify the process of obtaining experimentally challenging material properties of the mussel byssus system or any other compound material system. The results show that the force–displacement curves based on the identified material properties from both inverse approaches are in good agreement with the micro-tensile test results, which demonstrates the feasibility of inverse methods for identifying the material properties of biomaterials. Notably, the biological characteristics of mussel plaques contribute to the variability in plaque underside shapes. Even plaques from the same mussel exhibit different structures depending on the temperature or PH levels.^28,40 This investigation focuses on the application of the finite element (FE) inverse method, where the selected plaque had a nearly heart-shaped underside. However, in the study of Desmond et al. (2015) or Pang and Liu (2023), the undersides of the plaques showed an oval or square shape.^12,13 Therefore, the integration of the inverse method with basic mechanical and optical tests is essential for a comprehensive analysis. Furthermore, it is pertinent to highlight that our experimental protocol only applied a stretching angle of 90° rather than aligning with the natural growth direction of mussels, as the research focus is to assess the feasibility of inverse methods to obtain the material properties of the mussel byssus system. Pang and Liu (2023) have conducted a detailed analysis to explain why the natural growth angle of mussel plaques is close to 15° using directional tensile experiments and FE simulations.¹³ Additionally, the deformed side shapes and the potential detachment of the plaque bottom surface were also included in the FE inverse approach to further increase the accuracy of the results. Specifically, in the study conducted by Waite et al. (2019), the spring constant and Young's modulus of the mussel plaque were obtained by simplifying the stress–strain curve of the combined distal and plaque system into three linear regions,²⁶ while the consistency between our experimental and inverse curves suggests that the third-order Ogden model may be a more accurate option for characterising the material properties of distal threads, and the Neo–Hookean model may be a better choice for describing the material properties of both proximal threads and mussel plaques.

In the present study, the tensile behaviour of the different portions of the mussel byssus system, as well as the microstructures of the different portions derived from micro-tensile tests, are also in agreement with the results of previous studies.^9,10,28,41 The unique folds and furrows found in the proximal portion might be one of the reasons why the proximal portion is more ductile than the distal portion, as they allow the proximal portions more room to stretch and expand, helping to dissipate the strain energy (Fig. 7c and e). Meanwhile, the variety in protein composition could be another reason to explain the differences in the tensile properties of different portions. The biochemical composition of the mussel byssus system has been well researched, particularly the protein composition and metals.^{6,35,36,42,43} The distal portion is enriched in PreCol-D, whereas the proximal part is enriched in PreCol-P. PreCol-P with elastin-like folding sequences in the flanking domains is embedded in the matrix proteins of the proximal portion to provide greater elasticity. While preCol-D with spider silk-like sequences in the flanking domains could be reinforced with a β-sheet crystalline particle matrix to provide greater strength and toughness in the distal portion.⁴⁴ Additionally, the protein composition of the plaques differs from that of the threads, which is mainly composed of MFPs that help in the adhesion of the mussel plaque. The material properties, chemical composition, and microstructure of these three portions contribute to a better understanding of the mechanism and energy consumption underlying the composite characteristics of the mussel byssus system.

Although identifying the material properties of each portion in the mussel byssus system is a significant first step in analysing the mechanism of how mussels remain stable on the seafloor against persistent waves, future investigation is needed once the material properties have been obtained. For example, mussel threads are a composite of an external cuticle and a fibrous core, whereas mussel plaques are a composite of a cuticle and a porous core. It is not yet clear how the cuticle layer and core affect tensile behaviour, respectively. Meanwhile, several questions have not been addressed, such as the transition mechanism between the proximal portion and the distal portion, the dynamic plaque detachment process within the inverse approach and the fracture properties of these different portions. In addition, the study by Wilhelm et al. (2017) found that due to various forms of microstructural damage, the force–displacement curves of the mussel byssus system under multiple loading cycles showed a hysteresis behaviour, but interestingly the material strength remained comparable to that of the monotonically loaded mussel byssus system.⁴⁵ In our research, the mussel byssus system was assumed to be hyperelastic materials without considering any microstructural damage before the final failure. Therefore, how to further simulate and analyse the damage and self-healing of the mussel byssus systems under multiple loading cycles using inverse methods is also an interesting future research direction.

In our study, we introduced analytical and FE inverse methods to obtain the material properties of the constituents of the mussel byssus system and validated the feasibility of these methods, which provide a new approach to analyse the composite behaviour of the mussel byssus system and any other biocomposite material for future biomaterial mechanism research. Meanwhile, the diversity of mussel byssus system mechanisms also provides many insights for mimetic applications in engineering and material design, such as the design of anchor systems using the feature of the mussel byssus system to reduce fatigue in the mooring lines.

6 Conclusions

In this research, both analytical and FE inverse methods have been applied for the first time to determine the individual material properties of the mussel plaque and the proximal portion. Those inverse methods were critical due to the inherent difficulty in obtaining these properties through conventional tests. The inverse methods utilised experimental force–displacement curves from different constituents of the mussel byssus system and side-deformed shapes of mussel plaques as optimisation metrics.

The findings indicate that these inverse techniques effectively identify the individual material properties of the proximal thread and the plaque using reasonably assumed constitutive material models with the sum of square residuals equal to 0.0004 (N²) and 0.01 (mm²), respectively. Meanwhile, the results support the accurate representation of the material properties of the distal portion with the Ogden hyperelastic model and the proximal and plaque portions with the Neo-Hookean hyperelastic model. Additionally, this study has also conducted a comprehensive characterisation of the microstructure and material properties of the three different portions within the mussel byssus system. The distal portion exhibits elevated strength but reduced ductility in contrast to the proximal portion. The variations in the tensile behaviour among these portions could be attributed to differences in their microstructure or chemical composition. Finally, the study has also discussed the application of the FE inverse method for the calculation of the potential plaque bottom surface detachment radius.

The material parameters obtained from the inverse methods offer an opportunity for in-depth analysis of the material behaviour of the mussel byssus system. Future studies could focus on further improving the precision of inverse approaches by assessing the transition mechanisms between different portions and monitoring the dynamic detachment processes of the plaque bottom surface. Additionally, it will be necessary to determine the fracture characteristics of these three portions in the future and simulate the entire fracture tensile process using a comprehensive 3D finite element model of the mussel byssus system.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work was funded by the Leverhulme Trust Research Grant Scheme, UK (No. RPG-2020-235) and NanoCAT Research Grant, the University of Nottingham.

Notes and references

1. B. P. Lee, P. B. Messersmith, J. N. Israelachvili and J. H. Waite, Annu. Rev. Mater. Res., 2011, 41, 99–132.
2. P. V. Perepelizin and D. Boltovskoy, Biofouling, 2011, 27, 667–674.
3. J. D. Witman and T. H. Suchanek, Mar. Ecol.: Prog. Ser., 1984, 16, 259–268.
4. Y. Li and Y. Cao, Nanoscale Adv., 2019, 1, 4246–4257.
5. Q. Lin, D. Gourdon, C. Sun, N. Holten-Andersen, T. H. Anderson, J. H. Waite and J. N. Israelachvili, Proc. Natl. Acad. Sci. U. S. A., 2007, 104, 3782–3786.
6. D. S. Hwang, H. Zeng, A. Masic, M. J. Harrington, J. N. Israelachvili and J. H. Waite, J. Biol. Chem., 2010, 285, 25850–25858.
7. J. Allen, M. Cook, D. Jackson, S. Preston and E. Worth, J. Molluscan Studies, 1976, 42, 279–289.
8. D. Crisp, G. Walker, G. Young and A. Yule, J. Colloid Interface Sci., 1985, 104, 40–50.
9. E. C. Bell and J. M. Gosline, J. Exp. Biol., 1996, 199, 1005–1017.
10. E. Carrington and J. M. Gosline, American Malacological Bulletin, 2004, 18, 135–142.
11. J. R. Burkett, J. L. Wojtas, J. L. Cloud and J. J. Wilker, J. Adhesion, 2009, 85, 601–615.
12. K. W. Desmond, N. A. Zacchia, J. H. Waite and M. T. Valentine, Soft Matter, 2015, 11, 6832–6839.
13. Y. Pang and T. Liu, arXiv, 2023, preprint, arXiv:2305.15129 DOI:10.48550/arXiv.2305.15129.
14. N. Holten-Andersen and J. Waite, J. Dent. Res., 2008, 87, 701–709.
15. J. H. Waite, X.-X. Qin and K. J. Coyne, Matrix Biol., 1998, 17, 93–106.
16. J. Gosline, M. Lillie, E. Carrington, P. Guerette, C. Ortlepp and K. Savage, Philos. Trans. R. Soc. London, Ser. B, 2002, 357, 121–132.
17. A. Hagenau, H. A. Scheidt, L. Serpell, D. Huster and T. Scheibel, Macromol. Biosci., 2009, 9, 162–168.
18. M. Areyano, E. Valois, I. Sanchez Carvajal, I. Rajkovic, W. R. Wonderly, A. Kossa, R. M. McMeeking and J. H. Waite, J. R. Soc., Interface, 2022, 19, 20210828.
19. M. J. Harrington, H. S. Gupta, P. Fratzl and J. H. Waite, J. Struct. Biol., 2009, 167, 47–54.
20. C. Sun, E. Vaccaro and J. H. Waite, Biophys. J., 2001, 81, 3590–3595.
21. M. J. Harrington and J. H. Waite, J. Exp. Biol., 2007, 210, 4307–4318.
22. Q. Xu, M. Xu, C.-Y. Lin, Q. Zhao, R. Zhang, X. Dong, Y. Zhang, S. Tian, Y. Tian and Z. Xia, Adv. Sci., 2019, 6, 1902043.
23. S. Zechel, M. D. Hager, T. Priemel and M. J. Harrington, Biomimetics, 2019, 4, 20.
24. N. Aldred, T. Wills, D. Williams and A. Clare, J. R. Soc., Interface, 2007, 4, 1159–1167.
25. J. Smeathers and J. Vincent, J. Molluscan Studies, 1979, 45, 219–230.
26. N. Cohen, J. H. Waite, R. M. McMeeking and M. T. Valentine, Philos. Trans. R. Soc., B, 2019, 374, 20190202.
27. E. Vaccaro and J. H. Waite, Biomacromolecules, 2001, 2, 906–911.
28. E. Filippidi, D. G. DeMartini, P. Malo de Molina, E. W. Danner, J. Kim, M. E. Helgeson, J. H. Waite and M. T. Valentine, J. R. Soc., Interface, 2015, 12, 20150827.
29. M. U. Guide, Inc., Natick, MA, 1998, 5, p. 4.
30. J. H. Waite, E. Vaccaro, C. Sun and J. M. Lucas, Philos. Trans. R. Soc. London, Ser. B, 2002, 357, 143–153.
31. R. W. Ogden, Philos. Trans. R. Soc. London, Ser. A, 1972, 326, 565–584.
32. R. S. Rivlin, Philos. Trans. R. Soc. London, Ser. A, 1948, 241, 379–397.
33. D. S. Schnur and N. Zabaras, Int. J. Numer. Meth. Eng., 1992, 33, 2039–2057.
34. W. Sun, Z. Yue, Z. Wen and M. Li, J. Mater. Res. Technol., 2023, 22, 3132–3145.
35. K. J. Coyne, X.-X. Qin and J. H. Waite, Science, 1997, 277, 1830–1832.
36. X.-X. Qin, K. J. Coyne and J. H. Waite, J. Biol. Chem., 1997, 272, 32623–32627.
37. J. Wang and T. Scheibel, Biotechnol. J., 2018, 13, 1800146.
38. H. A. Price, J. Zoology, 1981, 194, 245–255.
39. G. M. Moeser and E. Carrington, J. Exp. Biol., 2006, 209, 1996–2003.
40. J. H. Bernstein, E. Filippidi, J. H. Waite and M. T. Valentine, Soft Matter, 2020, 16, 9339–9346.
41. E. Vekhova, Biol. Bull., 2021, 48, 1443–1451.
42. M. H. Suhre, M. Gertz, C. Steegborn and T. Scheibel, Nat. Commun., 2014, 5, 3392.
43. C. N. Schmitt, Y. Politi, A. Reinecke and M. J. Harrington, Biomacromolecules, 2015, 16, 2852–2861.
44. A. V. Golser and T. Scheibel, Fibers, 2018, 6, 21.
45. M. H. Wilhelm, E. Filippidi, J. H. Waite and M. T. Valentine, Soft Matter, 2017, 13, 7381–7388.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d3sm01677c

---