{'fraction_non_alphanumeric': 0.04211747969001535, 'fraction_numerical': 0.01591104788289469, 'mean_word_length': 4.05143721633888, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
J-shaped stress-strain diagram of collagen fibers: Frame tension of triangulated surfaces with fixed boundaries

We present Monte Carlo data of the stress-strain diagrams obtained using two different triangulated surface models. The first is the canonical surface model of Helfrich and Polyakov (HP), and the second is a Finsler geometry (FG) model. The shape of the experimentally observed stress-strain diagram is called J-shaped. Indeed, the diagram has a plateau for the small strain region and becomes linear in the relatively large strain region. Because of this highly non-linear behavior, the J-shaped diagram is far beyond the scope of the ordinary theory of elasticity. Therefore, the mechanism behind the J-shaped diagram still remains to be clarified, although it is commonly believed that the collagen degrees of freedom play an essential role. We find that the FG modeling technique provides a coarse-grained picture for the interaction between the collagen and the bulk material. The role of the directional degrees of freedom of collagen molecules or fibers can be understood in the context of FG modeling. We also discuss the reason for why the J-shaped diagram cannot (can) be explained by the HP (FG) model.


I. INTRODUCTION
The mechanical properties of macroscopic membranes, such as human skin, have been extensively studied experimentally for a long time [1][2][3]. One interesting mechanical property is the stress-strain diagram. This diagram is called "J-shaped" because of its plateau (linear behavior) in the small (large) strain region [4][5][6][7][8]. This J-shaped curve is quite different from the curve expected from the theory of elasticity, and it is also different from the curve observed in rubber elasticity [9]. Moreover, from an engineering perspective, this non-linear behavior attracts a considerable amount of attention for biomaterial functional technology [10,11]. For these reasons, many efforts have been devoted to understanding the origin of such a specific and unusual response to external forces. However, the mechanism still remains unclear, although it is widely accepted that the internal structure such as the collagen degrees of freedom [12][13][14], the notion of collagen network [15][16][17][18], and the notion of fibers [19,20] play essential roles in the J-shaped behavior.
In this paper, we use surface models for membranes, such as the Helfrich and Polyakov (HP) model [21][22][23][24][25][26] and a Finsler geometry (FG) model [27,28], to calculate the stress-strain diagram. The purpose of this study is to clarify the J-shaped behavior from the perspective of the theory of two-dimensional surfaces, which undergo thermal fluctuations. In such twodimensional surface models, the stress τ can be obtained as the frame tension of the surface that spans the fixed boundaries [29][30][31]. We will show that the J-shaped curve of τ can be obtained in the context of the HP model. However, a linear response of τ against the strain is also detected in an intermediate region of the bending rigidity κ. This linear behavior at the low strain region contradicts with the existing experimental data [2][3][4]. In contrast, the J-shaped curve of τ can be obtained independently of κ in the FG model. From this result, we confirm that the stress-strain diagrams obtained using the FG model are consistent with the existing experimental data.
The FG model is an extension of the HP model; hence, the Hamiltonian is composed of the Gaussian bond potential S 1 and the bending energy S 2 [27]. Moreover, the Hamiltonian includes the sigma model Hamiltonian S 0 for variable σ, which represents directional degrees of freedom of collagen or some internal molecular structures such as liquid crystals (LCs).
This variable σ plays an important role in the J-shaped curve of the stress-strain diagram, just like the polymeric degrees of freedom in the aforementioned collagen network and fiber models. We also note that the variable σ in this paper corresponds to the one used in the FG model [28] to represent the directional degrees of freedom of LC molecules in 3D liquid crystal elastomers [32][33][34][35]. We also note that the FG model in this paper is identical to the one introduced in Ref. [36], where the surface tension and string tension of membranes are calculated on spherical and disk surfaces. In this paper, we use a cylindrical surface for calculating the diagram; thus, both the boundary conditions and the results in this paper differ from those reported in Ref. [36]. Let us assume that an external force f is applied to a square surface, the size of which is supposed to be L × H ( Fig. 1(a)). Let τ be the surface tension; then, we have f = τ L, and therefore, the accumulated surface tension energy is given by is the surface area. This surface area A p is the projected area of the frame, and therefore, A p is not always identical to the real surface area A if the surface is fluctuating. Thus, the surface tension τ is called frame tension if the real surface area deviates from the projected area A p due to the surface fluctuations [30,31]. This frame tension τ is the one that we would like to calculate in this paper. In this paper, not only macroscopic membranes such as human skin but also microscopic membranes are assumed as the research targets, where the thermal fluctuations are not always negligible.
A cylindrical surface is used for calculating the frame tension τ ( Fig. 1(b)). We use this cylindrical surface because the cylinder has no boundary except the one to which an external force is applied. A triangulated cylinder of size N = 297 is shown in Fig. 2(a), where the height H and the diameter D are assumed to be identical: H = D. Let N 1 (N 2 ) be the total number of vertices in the height (circumferential) direction; then, we have N = N 1 ×N 2 , H = √ 3(N 1 −1)a/2, and D = N 2 a/π, where a is the triangle edge length. Thus, the ratio , and all cylinders that we use in the simulations are characterized by this ratio.

B. Finsler geometry model
In this subsection, we introduce a FG model, which is identical to the one introduced in Ref. [36]. The outlines of the discrete model and the corresponding continuous model are shown in this subsection and in Appendix A, respectively, in a self-contained manner.
First, we introduce the variable σ i (∈ S 2 : unit sphere) to represent the directional degrees of freedom of liquid crystal molecules (or collagen molecules). The Hamiltonian of the FG surface model is simply obtained by replacing the surface metric g ab with a Finsler metric (see Appendix A). To describe the interaction between the variables σ themselves, we include the sigma model energy λS 0 in the Hamiltonian S with the interaction coefficient λ such that S(σ, r) = λS 0 + S 1 + κS 2 + U B , where σ || i (see Fig. 2(b)) is defined as In S 0 of Eq. (1), we assume the factor 3/2 for the non-polar interaction, because Lebwohl-Lasher potential for LCs includes this factor, although LCs are not always included in collagen fibers [37]. In our FG model, the variable σ represents the direction of collagen molecule as mentioned above. The collagen fiber is made of collagen fibrils, which is made of collagen molecules (a hierarchical structure), and therefore the fiber becomes relatively stiff [11]. In addition, the fibers are loosely connected by cross-linkers. For these reasons, the collagen fiber networks are always locally ordered. This is in sharp contrast to the case of polymers, which have not only crystalline but also randomly disordered structure.
Therefore, the collagen fiber networks change from locally ordered to globally ordered states when they are expanded by external tensile forces. This coarse-grained picture of locally ordered structure of fiber network is expressed by the energy term λS 0 with finite λ for σ in our FG model. The reason why the polar interaction is also assumed for σ is simply to compare the results with those of non-polar interaction. The coefficient λ of S 0 is fixed to λ = 1 in both polar and non-polar interactions in the simulations. The fact that λ is fixed to λ = 1 is the cause of locally ordered configuration of σ, although λ = 1 corresponds to the isotropic phase at least for small strain region.
The stiffness of the fibers can be measured by EI, where E is the Young modulus and I the second moment of area. This EI is called bending rigidity, which measures stiffness of macroscopic elastic materials. In contrast, the bending rigidity κ in Eq. (1) corresponds to stiffness of microscopic membranes such as red cells. Thus, κ in Eq. (1) should be simply considered as a microscopic parameter that can be controlled depending on the rigidity of fibers in consideration.
The vector N i is the unit normal vector of the surface at vertex i, and it is defined as where A j(i) and n j(i) denote the area and the unit normal vector of the triangle j(i) sharing the vertex i, respectively. Note that S 0 is implicitly dependent on r because σ || i depends on the surface shape. The expressions for the Gaussian bond potential S 1 and the bending energy S 2 are The derivation of these expressions from the continuous Hamiltonians is shown in Appendix A. The symbol ℓ ij is the length of bond ij, and v ij is given in Eq. (A2) (see also Fig. 11 in Here we should comment on the boundary condition for the cylindrical surface. Because of the definition of v 13 in Eq.(A2), the variable σ 1 at vertex 1 on the boundary cannot be vertical to bond 13 on the same boundary. In fact, if σ 1 · t 13 = 0 then we have v 13 = 0, and therefore γ 12 → ∞. This divergence of γ 12 implies that σ 1 never be vertical to the boundary.
Therefore, to remove such unphysical repulsive interaction with respect to the direction of tensile forces, we assume that the boundary vertices are able to move into the horizontal direction only slightly within small range δ B . This constraint for the boundary vertices is defined by the potential U B . The small value δ B is given by the mean bond length, and therefore we have This implies that the constraint potential U B is negligible in the limit of N → ∞: while σ can be vertical to the boundaries or parallel to the direction of tensile forces.
The discrete partition function Z is given by where Z(λ, κ; L) denotes that Z depends on the parameters λ, κ and the height L of the cylindrical surface. for those on the boundaries, are not constrained by the boundaries; therefore,

C. Formula for calculating stress-strain diagram
The surface position r is the variable that is integrated out in the partition function Z, and for this reason, Z becomes invariant under the change of the integration variable such that r → αr (α ∈ R). This property is called the scale invariance of Z, and it is used for calculating the stress-strain curve [29].
The scale invariance implies that Z should be independent of the scale parameter α [38]; The scaled partition function is given by where α −2 A p in Z(α; α −2 A p ) denotes the dependence of Z on α arising from the fact that the projected area A p is fixed. This dependence of Z on α implies that Z can be considered to be a two-component function. Thus, from Eq. (8), we obtain dZ/dα = ∂Z/∂α + The mean value of the Gaussian energy S 1 on the left-hand side is obtained by Monte Carlo (MC) simulations. However, the problem here is how to evaluate the final term (1/Z)∂Z/∂A p . To calculate this term, we assume that the free energy F is given by where τ is the frame tension as mentioned above. Because the corresponding partition function is given by Z = exp(−F ), we finally obtain τ = (2γ The problem now is how to obtain the projected area A p in this τ . Let D and H be the diameter and the height of the cylinder, respectively. Then, it is natural to define A p as A p = πDH if D is uniform in the sense that D is independent of the height position h of the cylinder. Note that D corresponds to L/π for the cylinder, such as the one in Fig. 1 However, the diameter D is expected to generally depend on h because the cylinder is not a three-dimensional one but rather a two-dimensional surface, and therefore, D at the height position h = H/2 may be different from D at h ≃ H or h ≃ 0, for example. For this reason, we use the diameter D 0 (= H 0 ) of the initial surface, which is used for the simulations of τ = 0. Thus, the formula for τ is given by

D. Monte Carlo technique
The multiple-dimensional integrations in Z are simulated using the standard Metropolis Monte Carlo technique [40,41]. The initial height H 0 (= D 0 ) for the frame tension τ in Eq. (12) is determined such that the equilibrium configurations satisfy τ = 0. From this definition of D 0 and that of A p , the frame tension τ in Eq. (12) is considered to be the nominal stress in the sense that τ is independent of π D the real length of the circumference of the cylinder. To be more precise, π D can be identified using the real length of the circumference of the cylinder only when the cylinder is sufficiently smooth and has no surface fluctuations. In the case of 3D LCE, the nominal stress is calculated with a constant sectional area, which is independent of the height of the cylinder [28]. In contrast, the projected area A p for τ in Eq. (12) is proportional to the height H of the cylinder. Note that the diameter D 0 (= H 0 ) is not always identical to D (the symbol is not used henceforth, for simplicity); however, the deviation between D 0 and D is expected to be small because the cylinder is constructed such that the diameter equals the height for τ = 0.

A. Snapshots
First, we show snapshots of surfaces of non-polar model in Fig. 3. From the snapshots, we confirm that the variable σ is locally ordered when H/H 0 = 1 for κ = 1 while it is globally

B. Canonical model
In this subsection, we present the results of the HP (or canonical) model and discuss why the canonical model is insufficient for explaining the existing experimental data of the J-shaped stress and strain diagram [1][2][3][4]. The canonical model is defined as where ℓ 2 ij = (r i − r j ) 2 is the bond length squares [23,26,29]. In Eq. (13), we use the symbols S −1 and S 3 for the canonical Gaussian energy and the bending energy, respectively, to distinguish them from S 1 and S 2 in Eq. (1) for the FG model. Note that these S −1 and S 3 are not assumed as Hamiltonians in the FG model; however, these quantities can be obtained (or calculated) from the surface configurations of the FG model.
As shown in Fig. 4(a), τ is linear with respect to H/H 0 for κ = 1.5 and κ = 1. This figure also shows that the shape of τ changes from linear to J-shaped when κ increases from κ = 1.5 which is defined as is also J-shaped if the corresponding τ is J-shaped ( Fig. 4(d) In this subsection, we evaluate the equivalence between the canonical model and the FG model for λ = 0. When λ is zero in the FG model, the variable σ becomes random, and hence, no anisotropy is expected on the surface [27,28]. Indeed, the FG model is an extension of the canonical HP model in the sense that the outputs of the FG model for λ = 0 are consistent with those of the canonical model. Note that the role of κ in the canonical model is not always the same as that in the FG model. Indeed, κκ ij plays the role of the bending rigidity in the FG model, whereas the constant κ is the bending rigidity in the canonical model. Moreover, the surface tension coefficient γ is fixed to γ = 1 in the canonical model, whereas in the FG model, γ ij plays the role of the surface tension coefficient. For these reasons, to clarify the relation between the coefficients in the canonical model and those in the FG model, we calculate the mean values of γ ij and κ ij such thatγ The reason for multiplying S 1 and S 2 by a factor of 3 in Eq. (15) is as follows. The sum of triangles ∆ in S 1 and S 2 in Eq. (4) can be replaced by the sum of bonds 2 ij because every term ℓ 2 ij or 1 − n i · n j is summed over twice in the sum ∆ . From this factor 2, the factor 1/6 in S 1 , and the expression of S 2 in Eq. (4), we include the factor 3 inγ andκ in Eq. (15).
We observe that the dependence ofγ on H/H 0 in Fig. 5(c) appears similar to that of τ in Fig. 5(a). Indeed, bothγ and τ have a plateau (linear behavior) for κ = 5, κ = 3 and κ = 1 (κ = 1.5). We also observe that the value ofγ is in the range 2 ≤γ ≤ 2.8, and it is larger than that of γ(= 1) of the canonical model. However,γ is included in S 1 (=γS −1 ), and therefore, we expect that this difference betweenγ and γ(= 1) does not make any difference between the stresses τ of the canonical and FG models. Indeed, τ in Fig Hamiltonian for the FG model such that S eff =γS −1 +κκS 3 , which can also be written as S eff =γ [S −1 +κ(κ/γ)S 3 ]. In these expressions of S eff , the multiplicative factor 3 is dropped for simplicity. Furthermore, because the factorγ in S eff can be dropped due to the scale invariance of Z, we finally haveS eff = S −1 +κ(κ/γ)S 3 .
From thisS eff , the differences between S 3 /N B of the canonical and FG models are understood. As shown in Fig. 5(d),κ is slightly larger than 1;κ > 1. However,κ is not included in S 3 , which is shown in Fig. 5(b) (κ is included only in S 2 (=κS 3 )) in contrast to the case ofγ, which is included in S 1 . For this reason, S 3 /N B of the FG model is expected to be smaller than that of the canonical model. However, the results are opposite; S 3 /N B of the FG model is slightly larger than that of the canonical model. This result can be understood from the effective HamiltonianS eff = S −1 +κ(κ/γ)S 3 described above. Indeed, the fact that κ/γ < 1 makes S 3 larger. In other words, the effective bending rigidity of the FG model quite analogous to the one reported in Ref. [42], where an XY model energy suppresses the crumpling transition on spherical surfaces, although the interaction between σ and the surface of the XY model in [42] is different from the one of the FG model in this paper.
The results of the nonpolar FG model are presented in Fig. 7, where λ is fixed to λ = 1.5.
These data are consistent with those of the polar model in the region of κ such as 1 ≤ κ ≤ 5.
For all values of κ assumed, τ has the J-shaped structure.
The mean triangle area A and the order parameter M for the variable σ defined by are plotted in Figs. 8(a)-(d). In Eq. (16), σ for the polar case is given by σ = | i σ i |/N.
A plateau can also be detected in A, like that in τ , in both the polar and nonpolar models, and the range H/H 0 of the plateau for A is almost identical to that for τ . The area A corresponds to the real surface area; hence, it is considerably different from the projected area A p , which is proportional to H/H 0 . This difference between A and A p implies that the radius of the cylinder shrinks in its plateau region. In fact, it is easy to understand that A has no plateau if the cylinder radius remains unreduced. The eigenvalues Σ of the tensor order parameter defined by are plotted in Fig. 9 for the non-polar model. The largest eigenvalue Σ 1 becomes Σ 1 → 1, and the other two eigenvalues Σ 2 and Σ 3 are expected to be Σ 2,3 → −0.5 if σ is completely ordered. We confirm also from Figs. 9 (a),(b) that the variable σ becomes ordered if H is enlarged. The behavior of ordering of σ is exactly consistent to that of M in Fig. 8 (d).  obtained under the assumption that the surface remains cylindrical in shape. However, the surface shape deviates from cylindrical and becomes very thin for the small κ region, such as κ = 0.5, and collapses into string-like configurations if κ is reduced to κ = 0.4. In such a very thin surface, the surface area becomes far different from the projected area. Moreover, for these highly wrinkled surfaces, τ appears to always be positive even for small H 0 . This is actually expected because the surface shrinks to a small sphere for sufficiently small κ.
For these reasons, we assume a relatively large bending rigidity (κ ≥ 0.75) such that the cylindrical surface shape is maintained in the range 1 ≤ H/H 0 ≤ 2. In contrast, for the region of large κ, which is denoted by κ → ∞, S 2 is expected to be zero, and therefore, the surface shape can be changed only in its tangential direction. In this case, only S −1 changes as H/H 0 increases, and there is no reason for τ to behave non-linearly with respect to H/H 0 .
However, in the case that S 2 is not always exactly zero, where the surface is expected to undergo buckling, we also expect a non-linear behavior in τ . However, in these large regions of κ, the model surface will be far from biological membranes. The nominal stress vs. strain ( ) of (a) blood vessel [5], (b) rat muscle [6], (c) collagen fibers [7] and (d) collagen hydrogels [8]. The solid triangle denotes the simulation data τ sim of Eq. (18) for the non-polar interaction.
In this subsection, we show that the simulation data can be compared to experimental stress-strain data of biological materials such as blood vessel, rat muscle, collagen fibers, and collagen hydrogels [5][6][7][8] (see Figs. 10(a)-(d)). First, we should comment on the unit of τ in Eq. (12) assumed for the simulations in more detail. In the simulations, the inverse temperature β(= 1/k B T ) is fixed to β = 1(⇔ k B T = 1), and under this unit the triangle edge length a is fixed to a = 1. This a corresponds to the lattice spacing in the lattice field theory language [39] and is suitably fixed such that the simulation data can be compared to the experimental data. However, the physical unit of τ is given by [N/m], which is different from the experimental one [P a] for stresses of macroscopic objects. For this reason, we obtain τ sim dividing the simulation data τ by a to compare τ with the experimental stresses. By including β, a and k B T in the calculation formula of τ , we have where the room temperature is assumed for T . Note that the simulation data τ in Eq. Young moduluses E exp and E sim can be obtained from their linear part of the corresponding experimental nominal stress and τ such that the following condition is satisfied: From this relation, the coefficient (4 × 10 −21 /a 3 ) is obtained and used to plot τ sim in Figs. 10(a)-(d). We find in the experimental data that the linear behavior terminates for large strain region (see Fig. 10(d)). In contrast, the simulation data τ sim behave only linearly for large strain region because no failure mechanism is implemented in the model.  Table I. The parameters a are approximately 10 times (or more) greater than the Van der Waals radius of atoms, and therefore these a are meaningful as the lattice spacing for the calculations of τ sim .

IV. SUMMARY AND CONCLUSION
We have studied the origin of the J-shaped stress-strain diagram using Monte Carlo simulations on triangulated surfaces. For such a non-linear behavior of the J-shaped diagram, it has been widely accepted that the collagen structure plays an essential role [12-17, 19, 20].
However, for the J-shaped diagram, no concrete result has yet been obtained in theoretical or computational evaluations of the curve from the perspective of two-dimensional surface models because the interaction between the collagen and the bulk material (including collagen itself) is too complex.
To understand the mechanism of the J-shaped diagram, we first calculate the frame tension τ of cylindrical surfaces using the canonical surface model of Helfrich and Polyakov (HP) [38,47,48]. From the Monte Carlo data of the HP model, we find that τ is J-shaped.
However, the J-shaped curve can be obtained only for some limiting cases, such as small and large bending rigidity κ regions. In fact, for the region of κ ≃ 1. The important point to note is that a structural change is essential for the J-shaped curve.
This structural change is associated with the directional degrees of freedom of σ, which has two different phases, such as ordered and disordered. A phase transition between these two phases, from the disordered phase to the ordered phase, is activated by an external force that expands the surface. In this expansion process, the external force changes the internal structure represented by σ in the small H/H 0 region. As a result of this structural change, the surface fluctuation property is altered such that a long-range correlation, expected for a certain region of κ in the canonical model, is suppressed due to the one-dimensional correlation of σ. Thus, the internal structural change during the process of surface expansion is the origin of the J-shaped stress-strain diagram of membranes, and this intuitive picture for the interaction between σ and the bulk polymer can be implemented in the FG surface model. We should note that the detailed information on the transition property of this internal structure and the dependence of the J-shaped curve on the internal phase transition remain to be studied. expected to share the same origin with the soft elasticity in 3D liquid crystal elastomers [32][33][34][35].
Let σ i (∈ S 2 : unit sphere) be the variable corresponding to the directional degrees of freedom of a polymer or molecule such as liquid crystals at the vertex i of the triangulated surface. Let t ij be the unit tangential vector of the triangle edge (or bond) ij, which connects the vertices i and j, such that t ij = (r j − r i )/|r j − r i |. Using this t ij , we define the tangential component of σ i along the bond ij by v ij = |σ i · t ij | .
We note that v ij = v ji in general.
Let 123 denote a triangle on M, and let the vertex 1 be the local coordinate origin of the triangle 123; then, the Finsler metric g ab on the triangle 123 is defined by Thus, by the replacements ∂ 1 r → r 2 − r 1 , ∂ 2 r → r 3 − r 1 , ∂ 1 n → n 0 − n 2 , ∂ 2 n → n 0 − n 3 , √ gd 2 x → 1 2 ∆ det g ab (A4) in Eq. (A1), we have Because there are three possible local coordinate origins on the triangle 123, all possible terms in S 1 and S 2 should be summed over with the coefficient 1/3. The sum over triangles Multiplying γ(= 1) by γ ij and κ by κ ij , we have γγ ij and κκ ij , which can be considered to be effective surface tension and effective bending rigidity. These quantities γγ ij and κκ ij are dependent on the position and the direction of the bond ij, although γ ij and κ ij are a part of energies S 1 and S 2 , respectively. This dependence of γ ij and κ ij on the position and the direction of the bond is the most interesting output of the FG model. Anisotropic coefficients are expected to play an important role for the anisotropy in LCE [49][50][51][52]. Note that both S 1 and S 2 are explicitly dependent on σ because γ ij and κ ij are determined by σ via Eq. (A2).