Mathematical modeling and simulations for large-strain J-shaped diagrams of soft biological tissues

Herein, we study stress-strain diagrams of soft biological tissues such as animal skin, muscles and arteries by Finsler geometry (FG) modeling. The stress-strain diagram of these biological materials is always J-shaped and is composed of toe, heel, linear and failure regions. In the toe region, the stress is zero, and the length of this zero-stress region becomes very large (≃ 150%) in, for example, certain arteries. In this paper, we study long-toe diagrams using two-dimensional (2D) and 3D FG modeling techniques and Monte Carlo (MC) simulations. We find that except for the failure region, large-strain J-shaped diagrams are successfully reproduced by the FG models. This implies that the complex J-shaped curves originate from the interaction between the directional and positional degrees of freedom of polymeric molecules, as implemented in the FG model.


I. INTRODUCTION
Biological materials such as muscles, tendons and skin are known to be very flexible and strong, and for this reason, these materials have attracted considerable interest with regard to the design of artificial materials or meta-materials [1,2]. The mechanical properties of these materials are of fundamental importance in their applications. The stress-strain diagram is a typical approach for characterizing the mechanical strength of these materials, and numerous experimental studies on this topic have been conducted [3][4][5][6][7][8][9][10][11][12][13][14][15][16][17].
It has been experimentally observed that the stressstrain diagram of soft biological tissues is J-shaped and that the curve is composed of toe, heel, linear and rapture (or failure) regions [3][4][5][6][7][8][9][10][11][12][13][14][15][16][17]. The failure region is beyond the scope of this paper and will not be taken into consideration. In the toe region, the stress is almost zero, implying that the materials freely extend without external forces, similar to the behavior of the soft-elasticity region of liquid-crystal elastomers [18][19][20][21][22][23][24]. The existence of this zero-stress region is the reason why we call the curve Jshaped. The origin of the shape of curve is qualitatively understood such that a variety of materials, e.g., myosin, titin, and collagen, and structures, e.g., the cross-linked, entangled or network structures of these materials, contribute to this phenomenon, and hence, the interaction range extends to different length scales similar to a multicomponent material [25]. The continuum mechanics approach based on the strain energy functional also successfully describes the J-shaped diagrams [5]. In the context of continuum mechanics, the non-linearity in the J-shaped curve is understood as hyperelasticity [26,27]. However, the mechanism of the interaction underlying the J-shaped curve is still unclear on a quantitative and a microscopic level.
In Refs. [28,29], we studied J-shaped curves by the Finsler geometry (FG) modeling technique and obtained Monte Carlo (MC) data consistent with previously reported experimental results, in which the toe length is up to 40% ∼ 50% on the strain axis. For these experimental J-shaped curves of small toe length, FG modeling successfully describes the diagrams.
However, the toe length reaches 150% in some biological materials [7]. The length of the toe region is generally, albeit not always, limited to less than 50% [3][4][5][6][7]. One of the reasons for such a large difference in the toe length comes from the difference in the component molecules. The main component of materials with a small (large) toe length is the collagen (elastin, etc.) molecule [1]. Moreover, the hierarchical structure, such as that in the collagen molecule, is also expected to be one of the reasons for the variety of toe lengths [2]. Due to this hierarchical structure of molecules, there must be a considerable difference in the interactions between the hierarchies. The mechanism for large toe-length materials is qualitatively different from that of small toe-length ones. As a result, we should check whether the J-shaped diagrams with a large toe length can be studied by the FG modeling technique or not.
In this paper, the existing experimental J-shape curves of biological tissues such as animal skin, muscles and arteries are compared with the simulation results. The experimental curves of these materials are grouped into two types: the group of diagrams with a small heel (Sheel) and the group of diagrams with a large heel (L-heel) (Figs. 1(a),(b)) [4][5][6][7]. The diagram with an S-heel ( Fig.  1(a)) is decomposed into two straight lines with different slopes, while the diagram with an L-heel is decomposed into two different linear lines and a curve for the heel between the two lines. Moreover, close to the failure region, some of the curves have a convex part where the collagen fibers start to break ( Fig. 1(b)).
We show that the curve with an S-heel is successfully reproduced by our two-dimensional (2D) FG model, while the curve with L-heel is not described by the 2D FG model but can be reproduced by a 3D FG model. The large strain curves with a convex part are found to be well-fitted by 2D FG model data. Our results in this A J-shaped diagram is composed of two different linear lines, and the region where two lines are smoothly connected is called the heel. The J-shaped diagrams are decomposed into two groups: the diagrams with (a) a small heel and (b) a large heel. The curve in (b) is convex upward in the large-strain region and hence is highly non-linear.
paper indicate that the FG modeling technique can be used to analyze a wide range of J-shaped stress-strain diagram of biological tissues such as tendon, skin, muscles and arteries. Since the main component of these materials is a polymer such as collagen fiber, the polymeric degrees of freedom can be coarse grained and are simply replaced by the variable σ(∈ S 2 : unit sphere) in the FG model. This simple coarse graining is key to understanding the mechanical properties of these biological materials mathematically.
We comment on the reason why we use the FG model instead of standard or canonical models, which are defined without the variable σ [30][31][32][33][34]. The problem is whether or not the results of the canonical 2D and 3D models are identical to those of the FG models. The answer to this question is that the FG model is more suitable than the canonical model. In fact, almost the same results are obtained for some limited range of κ; however, the stress-strain diagram obtained from the canonical 2D model for a certain finite value of κ is not J-shaped, as reported in Ref. [28]. Another reason is that only the FG model reproduces the diagrams that are convex upwards for the large-strain region. These are the reasons for why we use the FG modeling technique to analyze the experimental J-shaped diagrams. In addition to these features of FG modeling, it is important to note that the mechanical strength of real membranes, e.g., the surface tension, becomes dependent on the position and direction on the surface [3,5]. This property is the origin of why the FG model is considered more suitable for actual biological membranes than the standard surface models [35].

II. MODELS AND MONTE CARLO SIMULATIONS
We use cylindrical surfaces that are suitable for the calculation of the surface tension or the stress. The  Fig. 2(b) is constructed from tetrahedrons for the 3D FG model [36], and the surfaces inside and outside the 3D lattice are 2D cylinders that are exactly same as the one in Fig. 2(a) [37]. This 3D lattice is thin, and hence, all vertices are on the surface; there are no vertices inside the 3D structure. The lattice size is given by (N, N B , N T , N tet ) = (5022, 24624, 34182, 14580), where the first three symbols are the same as those for 2D lattice and N tet is the total number of tetrahedrons. This 3D lattice is topologically identical to the torus and has zero Euler number χ = N−N B +N T −N tet = 0. The height H is fixed during the simulations for the calculation of the tensile force ( Fig. 3(a)). The diameter D of the upper and lower boundaries is also fixed to D 0 , and this boundary condition protects the cylinder from collapsing for small bending rigidity in the simulations. It should be emphasized that this constraint for D on the boundaries is close to the experimental setup for the measurement of tensile force.

A. 2D model
In this subsection, we describe the 2D FG model. Although this 2D model is the same as the one in Ref. [28], we summarize the Hamiltonian here in a self-contained manner. The Hamiltonian is given by a linear combina- tion of four terms such that The variable r(∈ R 3 ) is the vertex position and represents the position of a polymer, such as collagen fibers. The direction of the polymer is represented by σ(∈ S 2 ), which has a non-polar interaction of the Lebwohl-Lasher type [38] described by S 0 ; hence, σ is identified with −σ. We should note that the edges of the triangles do not always represent linear polymers or polymer networks [39]. The triangles are simply introduced for the discretization of 2D materials [40][41][42][43][44][45]. Indeed, the triangle edges play a role as local coordinate axes for the discretization of the Hamiltonian. The variable σ || i in S 0 is defined by which is a component of σ i parallel to the tangential plane at the vertex i ( Fig. 3(b)). This tangential plane is determined by its unit normal vector N i ( Fig. 3(b)), which is defined such that 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 ( Fig. 3(c)). The Gaussian bond potential S 1 and the bending en-ergy S 2 are given by In these expressions, ℓ ij is the length of bond ij connecting the vertices i and j, and n i is a unit normal vector of the triangle i. The coefficient γ ij in S 1 is defined by using v ij , which is given by (see Ref. [36] for more detailed information on the discretization of S 1 and S 2 ). The quantities γ ij (= γ ji ) and κ ij (= κ ji ) are considered the position-and directiondependent surface tension and bending rigidity, respectively.
The potential U B allows the boundary vertices to move vertically in the z-direction within a small range ±δ B , which is fixed to the mean bond length. This constraint does not influence the results in the limit of N → ∞ because δ B /H is negligible in this limit because the mean bond length is independent of N , while H is proportional to N . The reason for this constraint U B is assumed to be avoiding a strong and non-physical force, which is suspected to appear when σ i aligns with the zdirection on the boundary. If σ i on the boundary aligns with the z-direction without U B , the corresponding v ij becomes v ij → 0 because of the definition of v ij in Eq. (5). Therefore, the corresponding γ jk becomes γ jk → ∞, and hence, S 1 → ∞. In this situation, the variable σ i never aligns with the z-direction, and therefore, U B is necessary for the well-definedness of the model. The partition function is given by where H is the height of the cylinder and is fixed during the simulation ( Fig. 3(a)). ∫ ∏ 2N1 i=1 dr i denotes the 1D integrations of the vertices on the boundaries, where N 1 is the total number of vertices on the upper and lower boundaries. The 2N 1 vertices are allowed to move along the circles of radius D 0 , and hence, the corresponding integration effectively becomes one-dimensional. The total number of remaining vertices is N − 2N 1 , and the positions of these vertices are integrated out by the 3D integrations represented by

B. 3D model
The 3D model Hamiltonian is defined on the 3D lattice discretized by the tetrahedrons shown in Fig. 4(a) [37]. Although the Hamiltonian is almost the same as that in Ref. [36], we briefly describe it in the outline below. The model is defined without the self-avoiding potential for the "surfaces" (not for the inside of the structure), and this is the only difference between the models in this paper and in Ref. [36]. The surface self-avoiding potential is a non-local potential and is time consuming for simulations. We expect that the results are not strongly influenced by whether this self-avoiding interaction is included or not because the upper and lower boundaries are fixed and the surface always remains relatively smooth. This is also expected for the 2D model, which has no self-avoiding potential.
The Hamiltonian S(r, σ) is defined by a linear combination of five different terms The variable r(∈ R 3 ) is the vertex position of a tetrahedron, and σ(∈ S 2 ) denotes the directional degrees of freedom of polymers exactly the same as in the 2D model. Each term shares the same role with the corresponding term in the 2D model. The definition of the Lebwohl-Lasher potential S 0 is slightly different from that of S 0 in Eq. (1), but the role of this term in the 3D model is identical to that of S 0 in the 2D model. The definition of S 1 is also slightly different from the 2D case; however, the continuous description of S 1 is the same (see Ref. [36]). In the coefficient of S 1 ,N is defined bȳ where n ij is the total number of tetrahedrons sharing the bond ij and N B (= ∑ ij 1) is the total number of bonds. The γ ij (tet) in S 1 is given by where the numbers 1, 2, 3, 4 denote the vertices of the tetrahedron in Fig. 4(a). In these expressions for γ ij , the symbol v ij is defined by the same expression as in Eq. (5) using σ i and the unit tangential vector t ij along the tetrahedron edge ij ( Fig. 4(b)).
The term S 2 in Eq. (8) is different from that of the 2D model in Eq. (1); however, the role of S 2 in Eq. (8), i.e., to keep the tetrahedron shape almost regular for positive κ values, is the same as that in the 2D model. The symbol ϕ i is the internal angle of triangles ( Fig. 4(a)). The role of the potential U 3D is to protect the tetrahedron volume from being negative. This potential U 3D introduces a repulsive interaction between the vertices so that the tetrahedron is hardly collapsed, and hence, this U 3D shares the same role with S 2 in part. For this reason, the tetrahedron hardly deforms for positive κ values, and therefore, we assume small negative κ values in the simulations to make the J-shaped diagrams have a large strain. The potential U B is exactly same as in Eq. (1) for the 2D model, and for this reason, its definition is not written in Eq. (8).
The partition function Z 3D can also be defined for the 3D model; however, its description is exactly same as Z 2D in Eq.(7) except the actual number N 1 for the boundary vertices; to avoid duplication, it is not written here.

C. Formula for stress calculation
In both the 2D and 3D models, the stress in the stressstrain diagram is calculated from the principle of scale invariance of the partition function dZ/dα| α=1 = 0 (for simplicity, the subscript 2D in Z 2D is not written henceforth) [46]. This invariance simply originates from the fact that the integrations in Z are independent of its expression for r, the position of the material.
If we change r to αr with a positive number α, we have the scaled partition function such that where A p is the projected area of the surface and N 1 is the total number of boundary vertices, as mentioned in the previous subsection. The expression Z(α; A p (α)) implies that Z depends both explicitly and implicitly on α. In the Hamiltonian S(σ, αr), the only term that depends on α is S 1 : S 1 (α) = α 2 S 1 . We should note that the coefficient α 3N −4N1 in the right hand side of Eq. (11) comes from the 3D and 1D integrations such that From the abovementioned scale invariance of Z, we have d log Z/dα| α=1 = 0 and For the last term on the left hand side, we assume that A p (α) = α −2 A p for the dependence of A p (α) on α because the projected area A p is kept fixed under the scale change r → αr for the evaluation of the tensile stress τ (Fig.  5(a)). Moreover, to evaluate ∂Z/∂A p on the left hand side of Eq. (12), we naturally assume that the surface is sufficiently expanded. Under this condition, the free energy F of the surface is given by where A 0 is the area of the surface corresponding to the zero-tensile force [46]. Thus, we have the partition function Z = exp(−F ). Inserting this Z into Eq. (12), we have where D 0 is the diameter of the boundary. We should note that A 0 = πD 2 0 , where the initial height is given by H = D 0 . This surface tension τ in Eq. (14) is called the frame tension because τ depends only on the area of the frame on which the surface spans. We should note that the formula for τ of the 3D model is the same as Eq. (14) for the 2D model. This is because the thickness of the 3D cylinder the for 3D model shown in Fig. 2(b) is sufficiently thin that this 3D cylinder is regarded as a 2D surface.
The reason why D 0 = H is assumed in the configurations for τ = 0 is because the lattice is constructed under the condition D 0 = H with a regular triangle (see Fig. 2). The edge length is expected to be uniform and independent of the direction in the initial undeformed configuration if D 0 = H is satisfied, at least for λ → 0. Therefore, we have no reason to fix D 0 , for example, D 0 ̸ = H, although a non-zero λ is assumed in both the 2D and 3D simulations.
To compare the simulation result τ in Eq. (14) with the experimental data, we have to change simulation units to physical units. For this purpose, we explicitly use k B T and the lattice spacing a [28,47], which are suppressed by k B T = 1 and a = 1 in the expression for τ in Eq.
In this expression, a is varied to modify the simulation data τ , and the modified τ sim can be compared to the experimentally observed stress τ exp . The Young's modulus E can also be determined from the linear region of the simulation data τ by dividing τ by the strain. The obtained E is modified by the same factor in τ sim such that This E sim is directly compared to the experimental Young's modulus E exp , which is called the stiffness and determined from the linear region of τ exp [3,[6][7][8]. The value of a in this E sim is the same as that of a in τ sim . Therefore, E exp is not independent of τ exp . For this reason, we compare only τ exp with τ sim in this paper. We should note that among the parameters used in the simulation, not all of them are always comparable to physical quantities. The only quantities that can be compared to the experimental ones are τ sim and E sim . In fact, we assume that κ is negative in the 3D FG model. The reason for the negative κ is that the tetrahedrons hardly deform for positive κ values, i.e., where the obtained diagram has no toe region, as mentioned above.

D. Monte Carlo technique
The standard Metropolis MC technique is used to update the variables r and σ [48,51]. The variable σ is updated by using three different uniform random numbers, and the new variable σ ′ is defined independently of the old σ. The variable r is updated such that r → r ′ = r+δr with a small random vector δr. This vector δr is randomly generated in a sphere of radius d 0 , which is fixed for an approximately 50% acceptance rate.
On the upper and lower boundaries, the new position r ′ is constrained such that the diameter of the boundaries remains constant at D 0 , as mentioned above. For the 3D model, two different diameters, D in 0 and D out 0 , are assumed: one for the inner cylinder and the other for the outer one. The difference is given by D out 0 −D in 0 = √ 3⟨ℓ⟩, where ⟨ℓ⟩ is the mean bond length of the initial configuration for the simulation. In the discussions below, we use only D 0 for simplicity. The constraint on the diameter allows the vertices to move only along the circles of diameter D 0 . Due to this free movement of vertices along the circle, it is possible that the surface will become folded for a range of relatively small κ values when the height H is sufficiently small and close to D 0 (Fig. 5(b)). However, as will be seen in the snapshots of the surfaces, no folding is expected.
Another constraint is imposed on the boundary vertices by U B in Eq. (1). Under this U B , the vertex positions can move in the vertical direction (⇔ z-direction) within the small range δ B . This δ B is fixed to the mean bond length, as mentioned in the previous subsection.
The lattice size for 2D simulations is (N, N B , N P ) =  (10584, 31416, 20832), and the size for 3D simulations is (N, N B , N T , N tet ) = (9761, 48124, 66965, 28602). For this 3D lattice, which is constructed using the same technique as the lattice in Fig. 2(b), there are no vertices inside the structure, and all the vertices are on the surface.

A. Comparison with experimental data
We show in Fig. 6 the experimental stress-strain data of snake's skin reported in Ref. [3]. The skin of snakes, which is composed of collagen and elastin fibers, has a relatively large deformation, as expected from their typical body elongation and bending. The units of the stress τ exp are [MPa], and the snake skin is relatively strong. The toe region ranges from 50% to 125% depending on the body position from which the skin is sampled. The sampling positions of the data (Exp) in Fig. 6(a) and Fig.  6(b) are 40% and 60% distant from the snake's snout, where 100% corresponds to the length between the snout and the vent. The toe length of the plotted data in Figs. 6(a),(b) is relatively large and almost equal to 100% and 125% in units of strain. These curves are typical examples of diagrams with an S-heel, as mentioned in the Introduction, and they are composed of two different linear lines. The parameters used for the simulations are summarized below in Table I.
The stress τ sim in the 2D model in Fig. 6(a) is calculated by Eq. (15) using the simulation data τ , and τ sim is found to be almost identical to the experimental data τ exp . The parameter λ in Eq. (1) is fixed to λ = 1, and the bending rigidity κ is varied for the simulations of the 2D model in this paper. The assumed bending rigidities are κ = 0.6 and κ = 0.55 for the 2D simulations in Figs. 6(a) and 6(b), respectively. The lattice spacing a in Eq. (15) used for the fitting is a = 0.83 × 10 −9 in Fig. 6(a) and a = 0.85 × 10 −9 in Fig. 6(b), both of which are larger than the Van der Waals radius (∼ 1 × 10 −10 [m]), i.e., the typical size of atoms. This is the reason why we call the FG model a coarse-grained model.
The second and third sets of experimental data are of soft biological tissues such as diaphragm and arteries, to which the 2D model data are not always well fitted.  [6]. The strain of the toe region is approximately 50%, which is relatively short. For the large-strain region, the experimental data are slightly convex upwards, which is well fitted by the 2D simulation data except in the terminal failure region.
The data plotted in Fig. 7(a) are those obtained in the study of the mechanical diaphragm in a model of muscular dystrophy [4]. The data in Fig. 7(b) are the diagram measured along the circumferential axis of arteries [5]. In Ref. [5], Arroyave et al. investigated a methodology for mechanical characterization of soft biological tissues. In that work, biaxial tensile testing was performed, and methodologies for testing and data processing were proposed. Sections of three groups of arteries, such as the abdominal aorta, thoracic aorta and left subclavian artery, were selected from pigs in their study. They also analyzed the experimental data by using continuous mechanical models (Fung's Model and Holzapfel's Model). It was found that abdominal aortas can support higher stress before rupture due to the presence of collagen in the samples. It is also possible to understand that high percentages of elastin are the reason for the large strain in all groups. The curvature of the heel region in both experimental curves plotted in Figs. 7(a),(b) is relatively small compared to that in the experimental data shown in Figs. 6(a),(b). For this reason, the fitting of the simulation data of the 2D FG model is not good for these data, while the 3D simulation data are well fitted except in the failure region.
The next experimental data reported in Ref. [6] are of the periodontal ligament of the molars of rats at (a) 6 months and (b) 12 months of age (Figs. 8(a),(b)). In the case of this material, the strain is over 250% and is larger than that in the previous examples shown in Figs. 6 and 7. Moreover, for the large-strain region, the curve of the experimental data starts to bend and becomes convex upwards. The simulation data of 2D model are well fitted even for the convex part (except the failure region), although the 3D simulation data start to deviate from the experimental data in the convex region in Fig. 8(a). We should emphasize that only the 2D FG model produces results that are in good agreement with the experimental data. In fact, the results of the canonical model are always linear for the large-strain region, and they cannot be fit to the experimental curve with the convex part. This is a non-trivial difference between the FG and canonical models, as is the fact that the diagram of the canonical model becomes linear at a certain value of κ, while the diagram of the FG model is always J-shaped independent of κ [28].
The final examples of experimental data shown in Figs. 9(a),(b) are of the passive tension vs. titin strain of rat muscles, namely, skeletal and cardiac muscles [7]. Granzier et al. studied the mechanical properties of cardiac muscle by investigating passive tension and stiffness in a stretch-and-release process. They estimated the contribution of collagen, titin, microtubules, and intermediate filaments to the tension and the stiffness by gluing experiments after chemical treatments. Because irregular isoforms of certain proteins might cause heart diseases such as dilated cardiomyopathy, as reported earlier [49,50], it is important to understand why cytoskeletons play significant roles in passive tension and how they function. In their study, it was found that passive tension of cardiac cells increases slowly at the shorter end of the working range (1.9 to 2.2 µm) in the heart but increases more steeply at longer lengths up to 4.5 µm. They concluded that titin and collagen are the best two candidates for explaining the experimental results over a wide range of the lengths. On the other hand, microtubules and intermediate filaments do not develop passive tension significantly in the working range. In addition, as opposed to passive tension in skeletal muscles, cardiac passive tension does not show a plateau at long lengths around 4 µm. However, they argued that this plateau can be observed by excluding the effect of intermediate filaments.
These materials are very soft and flexible, and the titin strain is very large and ranges from 400% to 700%. The experimental curves are convex upward in the largestrain region, and the convex shape is more clear than that in the data shown in Figs. 8(a),(b). In addition to the previous data with the convex part, the results of the 2D FG model better fit the experimental data. The parameters assumed in the simulations are summarized in Table I. All values of a are sufficiently larger than the Van der Waals radius. These are the microscopic parameters and do not always correspond to actual physical quantities, as mentioned in Section II C.

B. Behavior of the variable σ and snapshots
To observe how the variable σ aligns, we calculate the order parameter M of σ by which represents the alignment of σ along the z axis [38]. We also calculate the eigenvalues of the tensor order parameter The largest eigenvalue Σ 1 of Q µν and M corresponding to several simulation results is plotted in Figs. 10(a),(b). For the small-strain region, Σ 1 and M slightly deviate from each other; however, they are exactly the same for the large-strain region. This implies that σ aligns in the z-direction to which the tensile force is applied. Snapshots of the surfaces are shown in Figs. 11(a)-(h), where the scales of the figures are different from each other because the height difference is very large. We find that the variable σ locally aligns at least in the configurations for H = D 0 in Figs. 11(a) and 11(e). This implies that σ changes in a similar manner to the directional degree of freedom of polymers, which undergo a transition from a locally ordered polydomain phase to a globally ordered monodomain phase. We should note that the locally ordered configuration in the small-strain region is not always random but almost ordered. In fact, a configuration of σ with a circular direction, for example, has no contribution to M because the surface is cylindrical. However, it is expected that this almost-ordered phase is consistent with the fact that linear objects such as collagen fibers in biological tissues are almost aligned even for the released and zero-strain configurations.

IV. SUMMARY AND CONCLUSION
We studied the large-strain J-shaped diagrams of biological membranes such as animal skin, muscles and arteries by 2D and 3D Finsler geometry (FG) models.
These materials are very soft, and the zero-stress region of the diagrams ranges from approximately 50% to 150%. Because of this zero-stress region, the diagram is called J-shaped. The J-shaped diagram is roughly composed of two linear lines except the failure region: one is the toe region, and the other is a linear region. The region where these two lines are smoothly connected is called the heel. Based on the word "heel", we can divide the experimental large-toe diagrams into two groups: diagrams with a small heel (S-heel) and diagrams with a large heel (Lheel). The diagrams with an S-heel are consistent with the results of the 2D model, while the diagrams with an L-heel are well fitted by the results of the 3D model. Moreover, the experimental diagrams with a convex part in the large-strain region can be fitted by the 2D FG simulation data. These observations show that the FG modeling technique is suitable for analyzing J-shaped diagrams of biological membranes.
The FG model is constructed by extending the linear chain to a 2D surface or a 3D body. In the 2D and 3D models, a Finsler metric is assumed in the Gaussian bond potential, and this Finsler metric plays a role in introducing the interaction between the variables σ and r, which correspond to the polymer direction and position, respectively. As a result of this modeling, the mechanical strength, such as the surface tension, effectively depends not only on the position but also on the direction inside the material. This is the main advantage of the FG model over its canonical counterpart. We should note that this FG modeling is only possible for 2D and 3D models because 1D space is a line with no non-trivial metric function.
The fact that the FG model was successfully applied to large-strain J-shaped diagrams indicates that highly nonlinear large-strain diagrams of polymers, include those with rubber elasticity, can also be targeted [39,52]. This will be a biology-inspired challenge in understanding new materials and their possible mechanical properties.

V. AUTHOR CONTRIBUTIONS
K.M. and S.G. studied and analyzed the existing experimental data and wrote the presentation section of the manuscript. H.K. performed the simulations and wrote the other part of the manuscript.