Is Stress Relaxation in Sea Cucumber Dermis Chemoelastic?

Abstract

Echinoderms, such as sea cucumbers, have the remarkable property of changing the stiffness of their dermis according to the surrounding chemical environments. When sea cucumber dermal specimens are constantly strained, stress decays exponentially with time. Such stress relaxation is a hallmark of visco-elastic mechanical behavior. In this paper, in contrast, we attempted to interpret stress relaxation from the chemoelasticity viewpoint. We used a finite element model for the microstructure of the sea cucumber dermis. We varied stiffness over time and framed such changes against the first-order reactions of the interfibrillar matrix. Within this hypothetical scenario, we found that stress relaxation would then occur primarily due to fast crosslink splitting between the chains and a much slower macro-chain scission, with characteristic reaction times compatible with relaxation times measured experimentally. A byproduct of the model is that the concentration of undamaged macro-chains in the softened state is low, less than $10\%$, which tallies with physical intuition. Although this study is far from being conclusive, we believe it opens an alternative route worthy of further investigation.

1. Introduction

Mutable collagenous tissue (MCT), found in all echinoderms, is a remarkable example of a collagenous tissue that can alter its mechanical state (e.g., stiffness) within a few seconds [1]. Echinoderms are an evolutionarily ancient animal phylum and include starfish, sea cucumbers, and sea urchins. Physiologically, mutability in MCT mechanics manifests in energy-efficient rigid posture maintenance with minimal muscle involvement [2], in irreversible softening processes, including autotomy (or defensive tissue detachment) [3], or in the generation of tensile force in feather stars [4]. MCT has attracted attention as a template for developing mechanically adaptive materials, with applications in biomedical, cosmetics, and bioinspired fields [1,5,6,7]. Understanding the dynamic mechanical and stress-relaxation behavior of MCT can have important implications in the biomedical materials and cosmetics field. The ability to regulate (increase or decrease) the stress relaxation rate in MCT (or its mimics) by controlling chemical bond breakage and reformation could enable the development of reconfigurable or injectable biomaterial implants in soft tissues. Especially in tissues with very low stiffness, such as neural tissue and its surroundings, regulating or increasing the stress relaxation or creep of MCT-based biomaterials would facilitate implantation without injury. In a similar manner, our understanding of the time dependence of mechanical properties in MCT would aid in optimizing the processes of the topical application and absorption of collagen-based cosmetics derived from MCT.

However, there is still a lack of clarity on how mechanical mutability—especially dynamic or time-dependent changes—occurs in MCT, specifically in models that link the extracellular matrix (ECM) to mechanical behavior.

1.1. Structure and Mechanics

Structurally, the major components of the MCT ECM are collagen fibrils with a characteristic axial electron-density periodicity D (≈65–67 nm) [1,8,9] and a heterotrimeric structure at the molecular level, along with fibrillin-rich microfibrils, and proteoglycans in the interfibrillar space. Neural effector cells (called juxtaligamental cells or JLCs), closely interpenetrating the ECM, are likely to play a significant role in modulating MCT mechanics by releasing effector molecules like tensilin or softenin [10,11], though the evidence is still indirect. Biomechanically, MCT stiffness (elastic modulus) ranges from 1–200 MPa on average [1], assessed using tensile tests or cyclic dynamic mechanical analysis. Failure strains of between 20–60% are intermediate between those of aligned fiber composites like tendons and tissues, with a wide range of fiber orientations like skin [1,12]. A change in mechanical properties can be induced in the laboratory by changing the ionic composition of the tissue bathing fluid or by mechanical stimulation. Regarding the modulus, elevations of 4–6x have been reported for immersion in potassium-enriched artificial seawater (KASW) [8] and a factor of 2–3x for mechanical stimulation [13,14]. A range of biological effector molecules isolated from MCT, like tensilin, NSF, softenin, and stiparin, have been implicated in the stiffening and softening of MCT [1], potentially by increasing interfibrillar crosslinking or stiffening the interfibrillar matrix. Local water content has also been known to change (reduce) during MCT stiffening [15], which can cause increased or closer fibril packing, which, in turn, may lead to greater interfibrillar stress transfer and stiffness. In this regard, fiber-composite models have been proposed to explain the variation in mechanics in MCT due to changes in the structure and properties of the ultrastructural constituents—the collagen fibrils and interfibrillar ECM components [8,16]. By using synchrotron X-ray nanomechanical imaging to quantify the fibril strains and orientation under different states of mechanical stimulation and quasi-static deformation, we showed that mutability in MCT arises from changes in the interfibrillar matrix properties and effective interfibrillar cohesion [8] rather than changes in collagen fibril properties.

1.2. Modeling the Time-Dependent Behavior

However, modeling the dynamic mechanical changes of MCT ECM—like stress relaxation, creep, or changes in mechanical state—in terms of the behavior of the ultrastructural building blocks is more challenging. MCT exhibits the visco-elastic behavior seen in most hydrated soft tissues, but the characteristics change with the altered mechanical state [17]. In order to estimate viscosity or flow in the ECM, earlier experiments tested the speed of bending of tissues under set weights and in the presence of stiffening factors [10] or with changes in calcium concentration [18], as the force is proportional to strain rate in Newtonian viscosity, the time for bending was taken as proportional to material viscosity. Structurally, however, the mechanisms by which the ECM mechanics change during typical time-dependent mechanical alterations, like stress relaxation, have been studied little at the ultrastructural level. Recently, we used experimental data on the time-dependent fibril stress-relaxation in chemically stimulated MCT (determined using small-angle X-ray scattering) to develop an ultrastructural model with shear-lag between the fully elastic collagen fibrils embedded in a viscous extra-fibrillar matrix, (which can contain fibrillin, proteoglycans, water, and effector molecules) [17].

In this previous work [17], we assumed a fully elastic behavior for the fibrils embedded into a visco-elastic interfibrillar matrix. A shear-lag model connected the fibrils and matrix. The matrix can only withstand shear stress, $\tau$, and the relationship with the shear strain, $\gamma$, is related to a convolution integral with the derivative of the interfibrillar matrix shear modulus $G_m\left(t\right)$, which was unknown.

$$\tau \left(t\right)=G_m\left(0\right)\gamma \left(t\right)-{\int }_0^t\frac{\mathrm{d}{G}_m(t-\tau )}{\mathrm{d}\tau }\gamma \left(\tau \right)\mathrm{d}\tau$$

(1)

Equation (1) is typical of visco-elastic behaviors [19]. We then obtained $G_m\left(t\right)$ and tried to explain the origin of $G_m\left(t\right)$ visco-elasticity by framing the visco-elastic parameters within an extended Doi-Edwards model. This interpretation agreed with some intuitive findings on the visco-elasticity of sea cucumbers: the increase in crosslinks and interchain friction from a standard-to-stiffening solution and the decrease in crosslinks and interchain friction from a standard-to-softening solution. However, these two results implied an apparent counter-intuitive conclusion on the elastic chain scission of the interfibrillar matrix. Specifically, we found that the molecular weight decreased during stiffening, and it increased during softening.

In order to address this discrepancy, we interpret stress relaxation as a chemical process in this paper, following [20], where elastically active chains and crosslinks break and reform dynamically. Here, by elastically active chains, we mean the chains in tension contributing to the shear modulus, $G_m\left(t\right)$, as in the molecular theory of rubber elasticity [21].

$$G_m\left(t\right)=N\left(t\right)kT$$

(2)

where k is the Boltzmann constant, T is the absolute temperature, and N is the number density of elastically active chains.

In a chemoelastic stress relaxation process, N changes over time due to chemical reactions. The chemoelasticity, which is sometimes referred to as chemorheology, has been applied, for example, to model aging-induced fracturing [22], a high-temperature elastomer response [23,24], and reversible crosslinking in hydrogels [25].

In summary, firstly, we built a staggered finite element numerical model. We extracted an approximate solution of the finite element model to obtain insights into the possible mechanisms behind tissue stress relaxation. From this solution, we inferred that a chemoelastic model with a shear modulus that degrades exponentially with time is very plausible. Based on this assumption, we used an optimization algorithm to derive the mechanical parameters. Finally, we explain these mechanical parameters for fibrils and interfibrillar matrix with first-order reaction kinetics to reveal the mechano-cracking mechanisms.

2. Materials and Methods

The experimental results, re-analyzed here using the chemoelastic model, were collected, reported, and published in previous papers [8,17]. Hence, only a summary of the experimental methods used in those papers is provided here as the experimental work was not part of the current study, and we refer to [8,17] for the details. Black sea cucumber (Holothuria atria) specimens were collected from a commercial wholesaler. After placing them in a $-{20}^{\circ }$ freezer, sections of tissue of $12\mathrm{mm}\times {\left(1.8\mathrm{mm}\right)}^2$ were dissected out. All studies were carried out in accordance with the Animals (Scientific Procedures) Act 1986 of the UK, including revision 2013: invertebrates (except cephalopods) are not considered protected species under the Act. Samples were incubated in artificial (ASW; control), potassium-enriched ASW (KASW), and calcium-free ASW (CaFASW) seawater, following the protocols published previously [8,26]. Synchrotron SAXS measurements were carried out at beamline I22, Diamond Light Source (Harwell, UK), using an X-ray energy of $12.46\mathrm{keV}$, a collimated beam size of $20\times 25\mathsf{\mu }{\mathrm{m}}^2$, and an exposure time of $0.5$ s per SAXS pattern. For the SAXS measurements, the tissues were mounted in a hydrated condition in a custom micromechanical tester in line with the X-ray beam at BL-I22 (developed by us [8]), using sandpaper to ensure adhesion to the grips. Stress-relaxation tests were carried out as described, and the SAXS patterns were azimuthally integrated to obtain radial line plots, from which the fibril strain was calculated as a percentage change in the 5th-order Bragg peak from the meridional fibril scattering, as described previously [8]. By combining the stress, tissue strain, and fibril strain plots as functions of time, the tissue and fibril-level mechanical data re-analyzed in this paper were obtained.

We formulate a shear-lag finite element model, where the interfibrillar matrix can only transmit shear stresses but not carry tensile loads (Figure 1). The implementation and solution details are in Appendix A.

We proceed in two steps. Firstly, we use approximated analytical results (Appendix B) to deduce qualitatively the mechanical properties of the fibrils and matrix. Such an analytical solution can be obtained only for a small number of fibrils; for more fibrils, the models can only be solved numerically.

Secondly, we use full numerical simulations to quantitatively derive the parameters of the elastic models for fibril and matrix for all the cases.

From the qualitative approximate analytical solution (Appendix B), we obtain

$${\sigma }_T\left(t\right)\approx \frac{1}{2}N{\rho }_F^2\sqrt{\frac{\varphi }{\pi }}{G}_m\left(t\right){ϵ}_T$$

(3)

where N is the number of fibrils, $\varphi$ is the fibril volume fraction, ${\rho }_F$ is the fibril aspect ratio, and ${ϵ}_T$ is the constant applied strain.

Equation (3) states that the tissue stress directly relates to the interfibrillar matrix. Experimentally, the tissue stress is captured by a two-term Prony series [8,17], as seen in Figure 2. Therefore, in the following, we assume that $G_m\left(t\right)$ is given by

$$G_m\left(t\right)=G_{\infty }\left(1+{\gamma }_1{\mathrm{e}}^{-t/t_{M_1}}+{\gamma }_2{\mathrm{e}}^{-t/t_{M_2}}\right)$$

(4)

where $G_{\infty }$ is the asymptotic shear modulus, ${\gamma }_1$ and ${\gamma }_2$ are dimensionless parameters of magnitude $\mathcal{O}\left(1\right)$, and $t_{M_1}<t_{M_2}$ are time constants.

Furthermore, experimentally, fibril strains are given simply by a one-term Prony series in the stiffened and softened state and by a two-term Prony series in the artificial seawater [8], as in Figure 3. From the approximate analytical solution in Appendix B,

$${ϵ}_F\left(t\right)\approx {ϵ}_{T}N{\rho }_F^2\sqrt{\varphi /\pi }\frac{G_m\left(t\right)}{E_F\left(t\right)}$$

(5)

where $E_F$ is the fibril Young modulus.

In order to minimize the total number of parameters involved in the analysis, we assume that $E_F$ is expressed by a one-term Prony series:

$$E_F\left(t\right)=E_{\infty }\left(1+e_1{\mathrm{e}}^{-t/t_{K_1}}\right)$$

(6)

with $E_{\infty }$ being the asymptotic fibril Young modulus, $e_1$ being a dimensionless parameter of magnitude $\mathcal{O}\left(1\right)$, and $t_{K_1}$ a relaxation time constant.

From experimental values [8], we know that fibril strain decays rapidly, with a time constant close to $t_{M_1}$. As we can see from Equation (5), the ratio $G_m\left(t\right)/E_F\left(t\right)$ is close to a one-term Prony series with a time constant $t_{M_1}$ if ${\gamma }_2\approx e_1$ and if $t_{K_1}\approx t_{M_2}$, meaning that the time constant for the fibril chemoelastic relaxation must be close to the largest time constant for the interfibrillar matrix chemoelastic relaxation. We use these insights to guide the numerical optimizations to quantitatively obtain the parameters of the elastic models for the fibrils and matrix. Therefore, assuming

$$G_M\left(t\right)=G_{\infty }+G_1{\mathrm{e}}^{-t/t_{M_1}}+G_2{\mathrm{e}}^{-t/t_{M_2}}{E}_F\left(t\right)=E_{\infty }+E_1{\mathrm{e}}^{-t/t_{M_2}}$$

(7)

we aim to calculate $G_{\infty }$, $G_1$, $G_2$, $E_{\infty }$, and $E_1$, and we take $t_{M_1}$ and $t_{M_2}$ from the experimental values in [8].

We calculate the parameters $G_{\infty }$, $G_1$, $G_2$, $E_{\infty }$, and $E_1$ in Equation (7) through a least squares minimization of the difference between the stresses and strains from the finite element model and the experimental values.

In all the calculations, we assumed a tissue specimen length of $L_T=10\mathrm{mm}$, a fibril length of $L_F=100\mathrm{\mu }\mathrm{m}$, a fibril aspect ratio of ${\rho }_F=1000$, and an applied tissue strain of ${ϵ}_T=50\%$. We used $N=100$ fibril elements and a final simulation time of $500\mathrm{s}$ with a time step of $\mathsf{\Delta }t=1\mathrm{s}$.

We have two sets of experimental values for matching the model: tissue stress and the fibril strain. Because the tissue stress depends only on the matrix shear modulus, we can proceed in two stages. In the first one, we minimize the difference between the tissue stress values, obtaining, in this way, $G_{\infty }$, $G_1$, and $G_2$. With these parameters at hand, we move on to a second curve, fitting the fibril strains, where we get $E_{\infty }$ and $E_1$.

3. Results

Figure 2 and Figure 3 display the numerical results of such minimizations for the three seawater conditions, namely the Artificial (ASW, taken as reference solution), the potassium-enriched (KASW), and the calcium ions-deprived seawater (CaFASW).

Table 1. Chemoelastical parameters for the interfibrillar matrix shear modulus and fibril Young modulus for the three seawater solutions.

Seawater	${\mathit{G}}_{\mathit{\infty }}$ (MPa)	${\mathit{G}}_1$ (MPa)	${\mathit{G}}_2$ (MPa)	${\mathit{E}}_{\mathit{\infty }}$ (MPa)	${\mathit{E}}_1$ (MPa)	${\mathit{t}}_{{\mathit{M}}_1}$ (s) [8]	${\mathit{t}}_{{\mathit{M}}_2}$ (s) [8]
ASW	6.6974$·{10}^{-8}$	6.1561$·{10}^{-8}$	9.8484$·{10}^{-8}$	27.59	0.1314	13.35	216
KASW	3.8535$·{10}^{-7}$	2.2968$·{10}^{-7}$	2.7412$·{10}^{-7}$	34.6294	19.5533	20.37	229.4
CaFASW	3.6746$·{10}^{-9}$	3.4478$·{10}^{-8}$	1.1543$·{10}^{-8}$	26.2296	0.0839	11.71	171.7

summarizes the calculated chemoelastic parameters for the interfibrillar matrix ($G_{\infty }$, $G_1$, $G_2$) and for the fibrils ($E_{\infty }$, $E_1$, $t_{K_1}$). Figure 4 and Figure 5 show the corresponding time histories as in Equation (7). Figure 4 shows the softening in CaFASW and stiffening in KASW compared to ASW, which is of about an order of magnitude for the interfibrillar matrix. For the fibrils, as in Figure 5, in the KASW solution, the fibril relaxes chemoelastically, with a time constant equal to the longest tissue relaxation time. It is also about two times stiffer than the fibril in ASW and CaFASW. In contrast, the ASW and CaFASW solutions display no fibril chemical relaxation, with the fibrils in ASW slightly stiffer than in CaFASW.

4. Chemoelastic Parameters in Terms of Kinetic Reactions

We now attempt to explain these findings through kinetic reactions involving the scission and recombination of links [20,27].

For the fibrils, we consider the following reversible reaction, where the chain scission occurs in equilibrium with chain recombination (Figure 6):

$$A\underset{k_R}{\stackrel{k_S}{\rightleftharpoons }}B$$

(8)

where A represents the macro-chains before scission, B represents the chains after scission, $k_S$ is the scission rate, and $k_R$ is the recombination rate. We assume that fibrils can only carry elastic loads if intact; therefore, the elastically active chains will depend on the concentration $\left[A\right]$ of species A.

$$E_F\left(t\right)=3\left[A\right]kT$$

(9)

where k is the Boltzmann constant, T is the absolute temperature, and the prefactor 3 appears because we assumed the fibril was incompressible.

For the interfibrillar matrix, we assume that two concurrent breakage mechanisms occur: the breaking of the intramolecular bond, namely chain scission, and the splitting of the crosslinks connecting different macro-chains (crosslink splitting). Additionally, we allow the possibility of simultaneous chain and crosslink recombination after these two splitting mechanisms (Figure 7). This picture is idealized, but it is a necessary simplification that introduces the least number of parameters (three reaction rates and three initial conditions) to fit the experimental data without too much indetermination.

$$\begin{array}{c}A\stackrel{k_1}{⟶}B\\ B\stackrel{k_2}{⟶}C\\ C\stackrel{k_3}{⟶}A\end{array}$$

(10)

where A represents the crosslinked macro-molecular chains, B represents the macro-chains after chain scission (broken bond inside the chain) that are still crosslinked, C represents the macro-chains after chain scission and crosslink splitting, $k_1=1/t_1$ is the chain scission reaction rate ($t_1$ is a characteristic time), $k_2=1/t_2$ is the crosslink splitting reaction rate and $k_3=1/t_3$ is the recombination rate.

We assume the matrix loses elasticity if it is in configuration C. Therefore, the elastically active chains are those in configurations A and B; even if the chains are scissored, the matrix can still carry loads because the chains are crosslinked.

$$G_m\left(t\right)=\left(\left[A\right]\left(t\right)+\left[B\right]\left(t\right)\right)kT$$

(11)

where k is the Boltzmann constant, and T is the absolute temperature.

4.1. Fibril Chemoelasticity

Solution (A28) of the kinetic equation indicates that the sum of the scission rate and recombination rate is equal to the inverse of the time constant $t_{K_1}$, which, in all three chemical conditions, is close to the largest tissue stress relaxation time constant $=t_{M_2}$ (Table 1). The higher the initial percentage of undamaged collagen chains $a_0$, the higher the proportion of the recombination rate (A32). Additionally, the balance of the recombination rate depends on the ratio $E_{F_{\infty }}/(E_{F_{\infty }}+E_1)$. Table 1 and Figure 5 show that $E_F\left(t\right)$ is constant for ASW and CaFASW. This finding means that fibrils do not show chemoelasticity and are only elastic, which is compatible with the previous results [8,17]. Interestingly, in KASW, a chemoelastic relaxation is possible: the ratio $E_{F_1}/(E_{F_{\infty }}+E_{F_1})$ is equal to 0.36. This means that for KASW

$${\kappa }_R=0.36a_0{\kappa }_S=1-0.36a_0$$

(12)

where ${\kappa }_S$ and ${\kappa }_R$ are the scission and recombination rates’ percentages. We are not able to infer $a_0$ from the experimental values. However, from Equation (12), we see that the recombination rate is always less than the scission rate, meaning that chain scission would occur much faster than chain recombination. The same applies to the interfibrillar matrix, as seen in the next subsection.

4.2. Matrix Chemoelasticity

After solving the kinetic reaction Equation (A36) and comparing the solutions’ time constants with the tissue stress relaxation’s time constants, we obtain three possibilities for the reaction rates, as shown in Figure 8: (i) fast recombination, slow crosslink splitting, and slow chain scission; (ii) fast crosslink splitting, slow chain scission, and slow recombination; (iii) fast chain scission, slow crosslink splitting, and slow recombination.

However, considering the coefficients of the Prony series of the matrix shear modulus in Table 1, we found that only one of the three possibilities is feasible. Indeed, from the coefficients in Equation (7), we can calculate the percentages (Equation (A65)) of the initial reagents in Figure 7. It turns out that the only reaction rates giving physically relevant percentages (from 0 to 100) are slow recombination, fast crosslink splitting, and slow chain scission.

Figure 9 shows the corresponding reaction times, i.e., the inverse of the reaction rates from the thick red line on the right of Figure 8. We see that the recombination times are slow or extremely slow in the order of hours or days, whereas crosslink splitting is, at most, 20 s, and the chain scission times are, at most, 500 s (8 min). When comparing these reaction times with the matrix shear modulus time constants in Table 1, we see that the crosslink splitting times are similar to the smallest time constants, $t_{M_1}$, and the chain scission times are comparable to the largest time constants, $t_{M_2}$, for all three chemical environments. Additionally, we see that crosslink splitting and chain scission in CaFASW (softening solution) are faster than those in KASW (stiffening solution) and ASW (standard solution).

It is impossible to infer from the model and the experimental values in our possession the exact initial conditions of the reagents. However, we can compute the possible percentages $a_0$ and $b_0$ (A65). Figure 10 shows percentages $a_0$ (crosslinked macro-chains, x-axis) and $b_0$ (scissored chains, y-axis) for the given reaction rates: $k_1$ (chain scission), $k_2$ (crosslink splitting), and $k_3$ (recombination).

We notice that for ASW and KASW, the initial conditions are such that the intact chains percentage, $a_0$, is always larger than the ones, $b_0$, with internal damage, in a ratio of ${}_0/b_0=2.33$. This finding could lead to speculation that the crosslinking between chains increases in going from a standard to a stiff state.

A completely opposite trend is seen for CaFASW. The initial conditions are such that $a_0$ is low, less than $10\%$, meaning that, initially, most chains are damaged: either scissored ($b_0$) or with their crosslinks removed ($c_0=1-b_0-a_0$). This result confirms the intuition that the chains undergo significant scission in going from a standard state to a soft one.

5. Conclusions

We have hypothesized a relaxation mechanism for the mutable connective tissue of a sea cucumber. Our explanation is alternative to the commonly believed assumption that visco-elasticity is the only mechanism at play. Indeed, in our previous modeling work, we found some discrepancies when we framed the visco-elastic parameters in a polymer theory, such as the Doi-Edwards model.

We then turned to a completely different point of view. We used a chemoelastic framework, where the stiffness varies because of a change in the concentration of polymer macro-chains participating in the elastic response. Such change is described by first-order kinetic reactions.

These reactions involve cracking mechanisms. For the fibrils, the hypothesized mechanism is the scission of collagenous chains. For the matrix, we hypothesized two competing mechanisms: the chain scission and splitting of crosslinks. In both cases, we incorporated the possibility of chain and crosslink recombination.

Our mathematical model revealed that in standard chemical environments, such as ASW, and in softening conditions (CaFASW), fibrils are only elastic, with no stiffness change over time. Only in a stiffening seawater solution (KASW) did we find chemoelastic behavior. If the hypothesized mechanism is correct, such fibril stress relaxation occurs so that scission is always faster than recombination.

Finally, our model indicated that in all three chemical environments, the interfibrillar matrix relaxes primarily through crosslink splitting, where reaction times are in the order of 10 s, similar to the shortest tissue stress relaxation times measured in the experiments [8]. Chain scission participates more slowly, in the order of 5–10 min. Even slower is recombination, which can take place over hours for barely damaged macro-chains or days for highly damaged macro-chains.

A final finding led us to speculate on the transitions from standard to soft and standard to stiff. For instance, after computing the initial conditions compatible with the experiments, it could happen that the standard-to-stiff form occurs by increasing the crosslinks in the interfibrillar matrix and the recombination of collagen chains in the fibril; in contrast, the standard-to-soft form might occur only through chain scission. This is speculative because our model and experimental measures refer to stress relaxation in fixed rather than changing chemical environments, although it seems intuitive.

We would like to remark that we have proposed hypothetical mechanisms, albeit backed by mathematical modeling. Further experimental tests, for example, measuring stress while the water solution changes from standard to soft or standard to stiff, would provide additional insights for a more advanced mathematical model.

Such a sophisticated model could be, for example, a coupled multi-physics three-dimensional finite element model, where the ultra-structure is represented more accurately and more complex reaction kinetics are considered.