Fractional viscoelastic models of porcine skin and its gelatin-based surrogates

Mechanics of Materials 177 (2023) 104559 Available online 13 January 2023 0167-6636/© 2023 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). Fractional viscoelastic models of porcine skin and its gelatin-based surrogates R. Mouˇ cka a , b , M. Sedlaˇ cík a , c , * , Z. P´ atíkov´ a d a Centre of Polymer Systems, University Institute, Tomas Bata University in Zlín, Tˇ rída T. Bati 5678, 760 01, Zlín, Czech Republic b Polymer Centre, Faculty of Technology, Tomas Bata University in Zlín, Vavreˇ ckova 275, 760 01, Zlín, Czech Republic c Department of Production Engineering, Faculty of Technology, Tomas Bata University in Zlín, Vavreˇ ckova 275, 760 01, Zlín, Czech Republic d Department of Mathematics, Faculty of Applied Informatics, Tomas Bata University in Zlín, Nad Str´ anˇ emi 4511, 76005, Zlín, Czech Republic ARTICLE INFO Keywords: Viscoelasticity Creep Fractional models Gelatin Porcine skin ABSTRACT Viscoelasticity of porcine skin and its material substitute, modelled by variously concentrated bovine gelatin, was determined in static (creep test) and dynamic (oscillatory test) mode by the means of rotational rheometry to obtain creep compliance and complex shear modulus. Mechanical properties characterization was also supplemented with large deformation compression test in order to determine and correlate shear and compression moduli of gelatin with its concentration dependence. Obtained data was fitted with fractional viscoelastic models (Poynting-Thomson, Maxwell) in order to quantify in detail gelatin’s transition from viscous-like behavior towards solid-like state with increasing gelatin concentration and hence crosslinking density. Potential of gelatin as biomaterial for skin surrogate was identified as well as a concentration region in which gelatin exhibits closest viscoelastic behavior to native porcine skin used. 1. Introduction During its lifespan skin has to be able to adapt itself to and resist the load/impact of many external mechanical stimuli varying significantly in the mode of deformation (White et al., 2013) (tensile (Yang et al., 2015), compression, friction (Bhushan et al., 2010; Leyva-Mendivil et al., 2017)), its intensity and time period of their impact. Additionally, skin has to maintain its properties in certain temperature and humidity window (Wu et al., 2006). All these requirements can be met to the large extent owing to unique viscoelastic properties stemming from skin inner structure/morphology (Depalle et al., 2015). Skin comprises of three layers epidermis, dermis and hypodermis (top to bottom) which are together about 1–3 mm thick. Detailed material knowledge of human skin in terms of its mechanical and stress-response properties has always been fundamental for many research fields ranging from cosmetics over dermatology and sensors applied directly onto skin (Zhou et al., 2019) to biomechanics and tissue engineering (Auger et al., 1998; Lee, 2000; Mansbridge, 2002; Vig et al., 2017). Another significant motivation fueling the research of viscoelastic properties of the skin is to be ultimately able to design material surrogate of the skin (Chanda, 2018). Although quite diverse materials can be used (Dabrowska et al., 2016; Morales-Hurtado et al., 2015), hydrogels in particular are a group of materials with high potential in this regard (Yi et al., 2021). Among them gelatin, a natural polymer extracted from collagen, is a promising candidate owing to its biocompatibility and mechanical properties tailorable through concentration and chemical crosslinking. Therefore, gelatin has been used in a number of biomedical applications (Alarcon-Segovia et al., 2021; Auger et al., 1995) including a component in additive technology of 3D printing of human skin (Jin et al., 2021). Its mechanical properties can be further improved through crosslinking, which is most frequently performed using glutaraldehyde (GTA) due to its easy availability and inexpensiveness (Bigi et al., 2001). In order to in detail determine skin’s viscoelastic properties variety of techniques (Pissarenko and Meyers, 2020) such as indentometry (Jachowicz et al., 2007), tension test (Yazdi and Baqersad, 2022), cutometry are employed, however rheological measurement remains a key method in this regard (Holt et al., 2008; Cheng et al., 2008; Pailler-Mattei et al., 2014; Verdier et al., 2009). Viscoelasticity of the material is typically determined through creep (Higgs and Ross-Murphy, 1990; Normand and Ravey, 1997) and relaxation tests and through dynamical measurement of its complex mechanical modulus, G*. All of these tests are designed to capture time response of the material to an impulse (aka forcing) which is static for creep (deformation response to * Corresponding author. Centre of Polymer Systems, University Institute, Tomas Bata University in Zlín, Tˇ rída T. Bati 5678, 760 01, Zlín, Czech Republic. E-mail address: [email protected] (M. Sedlaˇ cík). Contents lists available at ScienceDirect Mechanics of Materials journal homepage: www.elsevier.com/locate/mechmat https://doi.org/10.1016/j.mechmat.2023.104559 Received 25 August 2022; Received in revised form 13 November 2022; Accepted 12 January 2023 Mechanics of Materials 177 (2023) 104559 2 constant stress being applied to the sample) and relaxation (stress diminishing with time in the sample being deformed to a constant level) while time varying (typically harmonic deformation forcing within linear viscoelasticity range) in the case of oscillatory measurement. Standard approach to mathematical modelling of viscoelasticity typically involves the use of two basic elements – an elastic element represented by a spring and its viscous counterpart represented by a dashpot. These are linked either in series (Maxwell model), in parallel (Kelvin-Voight model) or more complicated network, such as in Zener or Poynting-Thomson models (Oyen, 2014). Although on the qualitative level these phenomenological models provide satisfactory results capturing the basic character of creep (Kelvin-Voight) or relaxation (Maxwell) test, they are far from being able to closely follow character of real material, whose steep initial increase and subsequent extremely slow levelling off are well beyond capabilities of a simple exponential function featured in the fundamental viscoelastic models (Bonfanti et al., 2020b; Macosko, 1996). Even though these shortcomings can be dealt with to a certain degree by including a number of basic elements (springs and dashpots) in the model, the resulting system still fails to fit the experimental data and also is overcomplicated due to large number of fitting parameters (two for each element) (Heymans and Bauwens, 1994). To overcome these disadvantages a new group of fractional viscoelastic models (Heymans and Bauwens, 1994; Long et al., 2018; Mainardi, 2010; Schiessel et al., 1995; Xu and Jiang, 2017) based on the idea of non-integer order of derivative have been proposed and proven to be an effective tool for description of experimental viscoelastic data (creep, stress relaxation) of the most diverse systems ranging from shape memory polymers (Fang et al., 2015) over gels (Faber et al., 2017; Holder et al., 2018; Rosalina and Bhattacharya, 2002; Zhang et al., 2018) and food (Mahiuddin et al., 2020) to biological (Bonfanti et al., 2020a; Carmichael et al., 2015; Craiem et al., 2006; Li and Tian, 2021; Mahiuddin et al., 2020) and geological systems (Di Giuseppe et al., 2009; Chen and Ai, 2020; Wang, 2021). The key concept is an introduction of another element called a spring-pot which exhibits behavior between the spring and the dashpot and by design captures power-law materials (Bonfanti et al., 2020b). The current work investigates a potential and feasibility of variously concentrated gelatin to model time dependent mechanical properties of skin. The study is based on detailed mechanical characterization obtained from viscoelastic rheological measurements as well as a compression test. Collected rheological data is fitted with suitable fractional viscoelastic models in order to quantify its time dependent behavior. 2. Experimental 2.1. Material 1. Porcine skin preparation Porcine ears were cleared and stored at −20 ◦C. Prior to sample preparation ears were defrosted naturally in room temperature and thoroughly cleaned under running water. Subsequently circular cartilage-free samples 25 mm in diameter were cut out. 2. Gelatin Gelatin from bovine skin, type B suitable for cell culture purchased from Sigma Aldrich was used in this study. 2.2. Gelatin samples preparation Gelatin samples were prepared in a concentration series spanning from 10 to 50 wt %. Corresponding amount of gelatin for preparation of x [%] concentrated gelatin was calculated (Eq. (1)) and dissolved in 10 mL of demineralized water at 60 ◦C for at least 15 min while mixing. Subsequently the solution was degassed (to remove air bubbles) in a vacuum temperature chamber and injected (due to high viscosity pouring was typically impossible) into a cylindric molds sized according to used test (rheometry diameter =25 mm, thickness =2 mm; compression test diameter =20 mm, thickness =8 mm). Gelatin was always let to cool down to lab temperature for 20 min prior to samples removal and other testing. x=mgelatine mwater (1) 2.3. Methods 1. Rotational rheometry All rheological tests were performed using an advanced modular rotational rheometer Physica MCR502 (Anton Paar, Austria) interconnected with water-cooled Peltier system P-PTD 200. A parallel-plate measuring system with a diameter of 25 mm (PP25) was used, while the plate gap was set to the sample thickness, i.e. 2.0 mm. All of the measurements were performed at 0.2 ±0.02 N normal force and at 25 ◦C. In frequency sweep tests angular frequency tested ranged from 0.1 to 100 rad s −1 with logarithmic sampling of 10 pt/decade. Strain applied to all investigated samples during this measurement was set to 0.5%, i.e. within the viscoelastic region determined from amplitude sweep beforehand. Creep test comprised measurement of strain for 300 s after exerting constant shear stress of τ 0 =10 Pa with sampling adjusted to the dynamics of the test. 2. Compression test Compression tests were carried out using a tensile test machine Testometric M350 5CT (UK) with a 10 kgf and a 1 kgf load cell employed for gelatin and porcine skin, respectively. Cylindric samples 20 mm across and 8 mm/1.5 mm (gelatin/porcine skin) thick were loaded at the center of two metal plates and compressed at a constant rate of 5 mm min −1 /1 mm min −1 (gelatin/porcine skin) while corresponding force was recorded. Measurements were performed at room temperature. 3. Scanning electron microscopy (SEM) Morphology of investigated systems, i.e. gelatin samples and porcine skin, was imaged using scanning electron microscope Nova NanoSEM 450 (FEI, Japan), prior to which the samples were first lyophilized overnight in order to remove water content and then sputter coated with a thin layer of gold. 4. Fractional viscoelastic models Fractional viscoelastic models work with a spring-pot (Fig. 1), which Fig. 1. A fundamental fractional viscoelastic element – a spring-pot; the value of α determines whether it behaves more like an elastic (spring) or a viscous (dashpot) element. R. Mouˇ cka et al. Mechanics of Materials 177 (2023) 104559 3 is a basic fractional element representing continuous transition between a spring and a dashpot and has a following governing equation (Bonfanti et al., 2020b): σ (t) = c α d α γ(t) dt α ,(2) where σ [Pa] denotes time dependent stress, γ [−] is time dependent strain, α is the order of Riemann-Louville fractional derivative, 0 < α <1 (also sometimes referred to as the fractional exponent) and c α [Pa s α ] is a “property” of the spring-pot. Substituting spring-pot in standard viscoelastic models, generalized fractional models are obtained. In our work generalized fractional Poynting-Thomson (FPT) model (Fig. 2) was used (Bonfanti et al., 2020b; Xu and Jiang, 2017). The creep compliancy of FPT model as a function of time is given as: J(t) = tλ η 3Γ(1+λ)+t α η 1 E α −β,1+ α (−t α −β η 1/ η 2),(3) where Γ is gamma function and E p,q (z) is Mittag-Leffler function Ep,q(z) = ∑ ∞ n=0 zn Γ(pn +q)(4) Complex modulus of the FPT model relevant for dynamic measurement of samples in oscillations is defined as follows: G∗( ω ) = η 3(i ω )λ[ η 1(i ω ) α + η 2(i ω )β] η 3(i ω )λ+ η 1(i ω ) α + η 2(i ω )β(5) 3. Results and discussion 3.1. Morphology Gelatin samples have rather uniform porous structure with pore size decreasing with gelatin concentration (Fig. 3 a-c). This in turn leads to denser system which is consequently also reflected in better mechanical properties in terms of higher shear G and compression E moduli. On the contrary to simple porous structure of gelatin, skin sample’s inner structure is much richer with complex hierarchical structure, from which only individual fibers are distinguishable in the SEM (Fig. 4 b), which is ultimately responsible for more graduate onset of deformation in the creep test (Fig. 5) as well as steeper increase of elastic part of complex shear modulus (Fig. 10) then observed in gelatin samples (Fig. 11). 3.2. Rheology 1. Creep test Typical creep test for porcine skin sample is shown in Fig. 5. Instantaneous elastic response of skin’s compliance to step shear stress is followed by gradually slowing increase of J with time, which within the investigated time frame showed no signs of reaching equilibrium compliance conditions. Although at the very beginning of the test a brief interval (0–0.01 s) marked with oscillations, where the skin behaves as an underdamped oscillator (so called creep-ringing (Ewoldt and McKinley, 2007; Goudoulas and Germann, 2016)), can be observed (Fig. 5 right), these manifestations are significantly smaller compared to gelatin. Porcine skin, owing to its extremely complex hierarchical structure comprising triple helical collagen molecules (1.5 nm in diameter) on the lowest level assembled in parallel into fibrils (50–500 nm thick), multiple of which form fibers and whole tissues (Gautieri et al., 2012), thus exhibits gradual creep of compliance with time when subjected to constant stress. On the contrary, gelatin samples (Fig. 6) exhibit pronounced creepringing in the first phase of the creep as gelatin together with rheometer forms an underdamped mechanical oscillator (Goudoulas and Germann, 2016). The oscillations tend to diminish sooner with gelatin’s concentration (Fig. 6 right) as the damping of the system increases due to gelatin’s changes in structure, i.e. higher crosslinking density and smaller pore size (Fig. 3). Apart from that increasing gelatin concentration also considerably affects value of immediate compliance, which decreases in non-linear fashion, as well as the shape of the creep curves, which flatten out. All of these manifestations clearly mark transition of gelatin from viscous-like to solid-like material with increasing concentration, which from material point of view stems mostly from higher crosslinking density and partially from decreasing pores’ size. In order to be able to quantify the character of the measured creep data in detail and correspondingly relate it to changes in the structure of the investigated system, data was approximated by the generalized FPT model (Fig. 2). Even though exceptional well agreement between FPT model and data was obtained, very symmetrical results in terms of the exponents of the parallel elements of FPT model and zero value for the third one in series (Table 1) hinted at the fact that possibly even simpler fractional model could suffice. Trying the very basic spring-pot model did not yield satisfactory fit, however two spring-pots in series i.e. generalized fractional Maxwell model (FMM - Fig. 7) did. Creep compliance for FMM has a following form: J(t) = t α η 1Γ(1+ α )+tβ η 2Γ(1+β).(6) Also, with one exponent in the Maxwell model being practically equal to zero (β =0.002) means that one of the spring-pots basically functions as a standard spring. Indeed spring-pot in series with a spring which accounts for immediate part of the compliance is sufficient for capturing the character of the porcine skin sample in creep, which is consistent with results from fitting using more complex FPT model where also the element in series functioned as a simple spring (λ =0). Therefore, the special case of fractional Maxwell model containing a spring-pot and a spring in series (Fig. 7) was employed for all further experimental creep data approximation. Although without rigorous physical meaning, for α <0.5 elasticity plays a leading role and the closer α approaches zero the more Hooke’s elasticity prevails (Di Paola and Zingales, 2012). Another interpretation relates the order of fractional derivative α to the Deborah number (Metzner et al., 1966) as follows: Fig. 2. Generalized fractional Poynting-Thomson viscoelastic model. R. Mouˇ cka et al. Mechanics of Materials 177 (2023) 104559 4 Fig. 3. SEM of variously concentrated gelatin (a–c: 10%, 30%, 50%). Fig. 4. SEM of porcine skin (a: sideview of the whole sample, b: detail of individual fibers). Fig. 5. Creep test of porcine skin in linear (left) and logarithmic (right) time scale to highlight creep-ringing at the initial phase of the test; experimental data (symbols connected with dashed line) fitted by the generalized fractional Maxwell viscoelastic model (full line). R. Mouˇ cka et al. Mechanics of Materials 177 (2023) 104559 5 De =1 α (7) with large De for solid-like materials while small De mark liquid. Thus, porcine skin with α =0.33 signifies rather elastic-type of behavior despite any clear signs of the creep curve reaching plateau in the investigated time frame. Indeed, especially in the case of natural materials, the final extremely long phase of the creep can often be misinterpreted as viscous flow due to limited (always finite) duration of the test. For gelatin samples FMM identifies instantaneous elastic part of compliance as well as in the case of skin sample yielding β =0 for all investigated concentrations (Table 2). Spring constant k (= η 2 ) in thus simplified FMM model is seen to increase non-linearly with gelatin concentration (Fig. 9) marking gradual increase of the samples!’ elasticity. Although transition of gelatin from more viscous-like system of at low concentrations to more elastic-like material as the loading approaches 50% is apparent from the flattening of the measured creep dependence, the usage of the FMM enables one to capture this transformation even quantitatively with α parameter, which gradually decreases from α 10% =0.83 to α 50% =0.27 (Table 2, Fig. 9) turning the spring-pot in the Maxwell model from dash-pot-like more into a spring (Fig. 1). This finding agrees with physical-mechanical expectations of the system rising from larger number of intermolecular interactions for more concentrated gelatin leading ultimately to higher crosslinking density (Fig. 8). Mechanical spectra, i.e. G*( ω ), were determined from subjecting samples to a sinusoidal shear deformation γ(t) = γ0ei ω t of amplitude γ 0 = 0.5% (within linear viscoelasticity region) and obtaining corresponding complex stress response. σ (t) = G∗γ0ei ω t= |G|eiθγ0ei ω t= (G′+iG′′)γ0ei ω t,(8) where G*( ω ) is frequency dependent complex dynamic shear modulus, real part of which (G′) defines storage (elastic) modulus in phase with Fig. 6. Creep test of gelatin in linear (left) and logarithmic (right) time scale to highlight creep-ringing at the initial phase of the test; experimental data (connected symbols) fitted by the generalized fractional Poynting-Thomson viscoelastic model (line); legend: gelatin concentration. Table 1 Fitting parameters of the generalized fractional Poynting-Thomson model. sample rheomode α β λ η 1 η 2 η 3 porcine skin creep 0.333 0.333 0.000 90.4 36.1 253.6 porcine skin oscillatory 1.000 0.115 0.354 16.2 6580 36 207 gelatine 10% oscillatory 0.004 0.270 1.000 1444 70 495 790 gelatine 20% oscillatory 0.000 0.261 0.698 3552 243 459 447 gelatine 30% oscillatory 0.000 0.324 0.282 5257 267 73 674 gelatine 40% oscillatory 0.000 0.361 0.256 8397 276 85 709 gelatine 50% oscillatory 0.001 0.421 0.243 12 501 303 89 703 Fig. 7. Generalized fractional Maxwell model (left) and its special case when β =0 (right). Table 2 Fitting parameters of the fractional Maxwell model employed for the creep test. Sample α β η 1 η 2 (=k) porcine skin 0.333 0.002 1272 2514 gelatine 10% 0.825 0.000 6574 1149 gelatine 15% 0.594 0.000 5670 1255 gelatine 25% 0.498 0.000 8399 2151 gelatine 30% 0.405 0.000 9687 3158 gelatine 50% 0.265 0.000 40 807 18 042 Fig. 8. Changes in interactions with gelatine concentration. R. Mouˇ cka et al. Mechanics of Materials 177 (2023) 104559 6 excitation deformation and imaginary part (G′′) gives loss modulus or viscous part responsible for energy dissipation. Although the creep test data is well fitted with it, the fractional Maxwell model failed to provide satisfactory results for frequency measurement of complex shear modulus and therefore FPT model (Eq. (5)) was employed with much better results even though both components (G′and G′′) of the complex shear modulus were approximated with the same set of parameters. Extracted parameters are given in Table 1. Even though the very same set of parameters for a concrete model determined from creep data fitting should work even for oscillatory data, in reality this is unfortunately not the case and even though the model is capable of fitting dynamic rheological data the parameters are quite different. This is at least partially due to numerous effects connected with sample geometry imperfections, loading and normal force application prior to measurement. Moreover, the measurement region is always restricted either in time (for creep) or forcing angular frequency (for oscillations) the result of which being mismatch between fitting parameters in each method. This transpires even in the case of complex modulus G* fitting, where basically two curves need to be fitted simultaneously. Nevertheless, at least qualitatively the fractional model captures its character. Spectra of both, porcine skin (Fig. 10) and gelatin (Fig. 11), systems exhibit higher values of storage modulus over loss modulus with tan δ= G′′/G′of about 0.3 and 0.05 for skin and gelatin, respectively. Particularly in case of gelatin, such low tan δ values indicate well-developed network, in which elasticity dominates over viscous flow, for which they are typical (Van den Bulcke et al., 2000). Storage and loss moduli of porcine skin both increase gradually with angular frequency which is characteristic for complex material comprising structural elements, e.g. fibers, fibrils and collagen Fig. 9. Fitting parameters of fractional Maxwell model for gelatin (symbols) and porcine skin (line). 2. Oscillatory test Fig. 10. Mechanical spectrum of a porcine skin; experimental data fitted with FPT model (left) and fitted FPT model extrapolated into wide frequency domain (right). Fig. 11. Mechanical spectra of gelatin fitted with FPT model; legend: gelatin concentration. R. Mouˇ cka et al. Mechanics of Materials 177 (2023) 104559 7 molecules, with different relaxation times. On the other hand, gelatin exhibits plateau or gently rising frequency dependence of G′for low (10% and 20%) or more concentrated (30% and above) samples, respectively. An exceptionally good fit of mechanical spectra with FPT model (Fig. 11) enables one to use the fitting parameters (Table 1) and predict the behavior of the material outside the measured frequency domain. Thus, master curves of variously concentrated gelatin samples were obtained (Fig. 12), from which it can be seen that the character of the gelatin’s spectrum considerably changes between concentrations 20 and 30%. A well-defined G′′ peak at low frequencies (10 −3 to 10 −2 rad s −1 ) observed for 10 and 20% is connected with relaxation of polymer chains’ disengagement (reptation), and marks cross-over frequency at which gelatin turns from liquid (G′′ >G′) into elastic state (G′>G′′). For more concentrated gelatin (from 30% on) G′>G′′ holds always as the density of crosslinks is too high for material to exhibit viscous flow even at very low frequencies. This is also reflected in the values of FPT model parameters with the spring-pot in series functioning at first as a dashpot (λ =1) with large viscosity ( η ≈5 ×10 5 Pa s) for 10% gelatin, while exhibiting more elasticlike behavior for 30, 40 and 50% gelatin (λ ≈0.25). Instantaneous elasticity present in the system due to physical crosslinks in gelatin’s chains is modelled by one of the spring-pots in parallel ( α =0). Its “modulus” ( η 1 ) was seen to increase from 1400 Pa to 12 500 Pa reflecting rise in the density of physical crosslinks in the system with gelatin concentration (Fig. 13). A fundamental topological parameter characterizing the polymer network introduced by de Gennes is the entanglement molecular weight, M e , which is defined as an average molecular weight between topological constraints (Degennes, 1971). It is typically inferred from plateau modulus G0 N often obtained as a value of storage modulus at frequency of loss modulus minimum (G0 N=G′( ω )G′′→ min) (Liu et al., 2006). Relationship between G0 N and M e is given as G0 N=4 5 ρ RT Me ,(9) where ρ [kg m −3 ] is density, R [J mol −1 K −1 ] molar gas constant and T [K] thermodynamic temperature. Availability of wide range mechanical spectra for gelatin, owing to FPT model predictions, allows one to identify G′′ minima and subsequently extract G0 N and calculate M e . Polymer network can be described in terms of density of its crosslinks ρ x [mol m −3 ] which can be estimated as: ρ x=|G∗| RT (10) and then used for the determination of mesh size: ξ= 6 πρ xNA 3 √,(11) where N A =6.022 ×10 23 [mol −1 ] is Avogadro’s constant. Entanglement molecular weight decreases with gelatin’s concentration (Fig. 13) as in the more concentrated gelatin there is higher incidence of interparticle interaction leading to creation of junction zones (Fig. 8). These are regions of partially reformed triple-helices (of gelatins’ polypeptide chains) stabilized by weak interactions (H-bonds), which act as physical crosslinks (Joly-Duhamel et al., 2002). Higher helix concentration consequently results in gelatin’s better mechanical properties manifested by higher shear and compression (shown) later on modulus. C. Compression test To expand the gelatin’s mechanical behavior to the macroscopic scale of large deformations a compression test of porcine skin and gelatin samples was carried out. Contrary to the tensile test for the compression test it is rather difficult to precisely identify at which point the breakdown of the sample occurs. Therefore, only initial phase, i.e. loading (compression) of the gelatin samples with pressure, which has been intentionally clipped at the strain of 60%, value attainable for all samples, is presented (Fig. 15). Porcine skin exhibits convexly shaped compression test curves (Fig. 14) with two linear regions. One is located at the beginning ( ε = 0–20%) of the test with compression modulus E 0 =1.75 kPa while the other appears at more substantial compression deformations (substantial scatter of the onset of the second phase between the four measured porcine skin samples is caused by imperfect nature of skin samples, whose top and bottom are not perfectly parallel) and has modulus of about 30 kPa (E 1 ). Both values approximately correspond to skin’s shear modulus |G*| determined at the ends of investigated shear rate range in rheological measurements. Compression tests of gelatin samples show region of linear elasticity for all investigated samples (10–50 wt % gelatin concentration) up to about 25% of strain above which the derivative d σ /d ε increases (Fig. 15). Compression modulus, E, of the individual samples was calculated from the linear region, namely for ε ∈(5%;10%), as a value of the derivative of the experimental engineer stress-strain curves: E=d σ d ε (12) Extracted modulus values as well as the compression test curves both show that in terms of compression deformation even lowest concentrated gelatin investigated (10%) has slightly better mechanical properties than porcine skin with more concentrated gelatin samples clearly exceeding compression modulus of porcine skin. Fig. 12. Extrapolation of fitted data with FPT model into wider frequency domain; legend: gelatin concentration. R. Mouˇ cka et al. Mechanics of Materials 177 (2023) 104559 8 Even though gelatin’s compression modulus, E, is about one order of magnitude smaller compared to the absolute value of complex shear modulus, |G*|, measured in oscillations at frequency of 10 rad s −1 , its concentration dependence has a quite similar slightly convex character (Fig. 16), which supports the idea of gradual reinforcement of the gelatin polymer network with crosslinks due to more frequent gelatin interactions as gelatin concentration is increased. Despite different deformation mode (shear versus compression) and magnitude (small versus large) the fundamental effect of gelatin concentration on its mechanical properties is maintained. To comprehensively assess viscoelastic properties of gelatin as a possible surrogate for natural skin it can be claimed that although much simpler in terms of its structural complexity, gelatin of appropriate concentration responds to both static (creep) and dynamic excitation in many regards in a way much similar to a real skin. This concentration for the used bovine gelatin was identified as 40%, at which complex dynamic modulus G* at relevant angular frequencies (5–10 rad s −1 corresponding to walking–running (Holt et al., 2008)) as well as parameters of FMM exhibited highest agreement. IV. Conclusions. Experimentally determined viscoelastic properties of porcine skin using rotational rheometry (creep and oscillatory test) were modelled with generalized fractional Maxwell model (two spring-pots in series) in order to improve upon accuracy and comprehensibility of the fit offered by traditional viscoelastic models used. Fitting revealed that even a special case of fractional Maxwell model comprising a spring-pot and a standard spring in series is capable of very precise approximation of the Fig. 13. Evolution of gelatin’s network characteristics (crosslinking density/black squares/, mesh size/blue circles/and entanglement molecular weight/black crosses/) with concentration. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.) Fig. 14. Compression test of porcine skin. Fig. 15. Compression test of gelatin; legend: gelatin concentration. Fig. 16. Compression and shear modulus versus gelatin concentration. R. Mouˇ cka et al. Mechanics of Materials 177 (2023) 104559 9 creep test data. Progressive decrease of α exponent captures gelatin transition from viscous-like to more elastic-like state with its concentration in the system. This finding is further supported by gradual increase of spring constant k in the model as well as both compression and shear moduli increasing with gelatin content. On structural level this transition of viscoelasticity more towards elastic part of the spectrum stems from higher density of intermolecular interactions between gelatin molecules. For fitting response to dynamic excitation (mechanical spectra) of material a more complex fractional model (Poynting–Thomson) was employed, which was able to fit gelatin exceptionally well with model’s parameters confirming an increase of its elasticity with gelatin concentration. This was also confirmed by gradual decrease of calculated entanglement molecular weight of gelatin with concentration. Improvement of mechanical properties, namely compression modulus, with concentration of gelatin was observed even at large scale deformations during compression test. Even though gelatin’s porous internal structure is much simpler compared to hierarchical structure to be found in biological system of porcine skin, it certainly has the potential to model the skin in terms of mechanical properties as has been shown for small deformations using rheometry. Performed large scale deformation compression test confirms this finding with even low concentrated gelatin exhibiting higher compression modulus than porcine skin. Thus, from the viewpoint of viscoelastic properties, gelatin and its modifications (e.g. through crosslinking) can be seen as suitable material for further testing (e.g. large deformation viscoelasticity) in the field of skin substitutes. CRediT author statement Robert Mouˇ cka: Conceptualization, Methodology, Formal analysis, Investigation, Writing – Original Draft, Visualization. Michal Sedlaˇ cík: Conceptualization, Investigation, Writing – Original Draft, Writing - Review & Editing, Supervision, Funding acquisition. Zuzana P´ atíkov´ a: Formal analysis, Writing - Review & Editing. Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Data availability Data will be made available on request. Acknowledgements The authors gratefully acknowledge project DKRVO [RP/CPS/2022/ 007] supported by the Ministry of Education, Youth and Sports of the Czech Republic. The authors also wish to thank the Czech Science Foundation [23-07244S] for their financial support. References Alarcon-Segovia, L.C., Daza-Agudelo, J.I., Rintoul, I., 2021. Multifactorial effects of gelling conditions on mechanical properties of skin-like gelatin membranes intended for in vitro experimentation and artificial skin models. Polymers 13. Auger, F.A., Valle, C.A.L., Guignard, R., Tremblay, N., Noel, B., Goulet, F., Germain, L., 1995. SKIN equivalent produced with human collagen. Vitro Anim. Cell Dev. Biol. 31, 432–439. Auger, F.A., Rouabhia, M., Goulet, F., Berthod, F., Moulin, V., Germain, L., 1998. Tissueengineered human skin substitutes developed from collagen-populated hydrated gels: clinical and fundamental applications. Med. Biol. Eng. Comput. 36, 801–812. Bhushan, B., Tang, W., Ge, S., 2010. Nanomechanical characterization of skin and skin cream. J. Microsc. 240, 135–144. Bigi, A., Cojazzi, G., Panzavolta, S., Rubini, K., Roveri, N., 2001. Mechanical and thermal properties of gelatin films at different degrees of glutaraldehyde crosslinking. Biomaterials 22, 763–768. Bonfanti, A., Fouchard, J., Khalilgharibi, N., Charras, G., Kabla, A., 2020a. A Unified Rheological Model for Cells and Cellularised Materials, vol. 7. Royal Society Open Science. Bonfanti, A., Kaplan, J.L., Charras, G., Kabla, A., 2020b. Fractional viscoelastic models for power-law materials. Soft Matter 16, 6002–6020. Carmichael, B., Babahosseini, H., Mahmoodi, S.N., Agah, M., 2015. The fractional viscoelastic response of human breast tissue cells. Phys. Biol. 12. Chanda, A., 2018. Biomechanical modeling of human skin tissue surrogates. Biomimetics 3. Chen, Y.F., Ai, Z.Y., 2020. Viscoelastic analysis of transversely isotropic multilayered porous rock foundation by fractional Poyting-Thomson model. Eng. Geol. 264. Cheng, S., Clarke, E.C., Bilston, L.E., 2008. Rheological properties of the tissues of the central nervous system: a review. Med. Eng. Phys. 30, 1318–1337. Craiem, D.O., Armentano, R.L., Ieee, 2006. Arterial Viscoelasticity: a Fractional Derivative Model, 28th Annual International Conference of the IEEE-Engineering-In. Medicine-and-Biology-Society, New York, NY, p. 4156. Dabrowska, A.K., Rotaru, G.M., Derler, S., Spano, F., Camenzind, M., Annaheim, S., Stampfli, R., Schmid, M., Rossi, R.M., 2016. Materials used to simulate physical properties of human skin. Skin Res. Technol. 22, 3–14. Degennes, P.G., 1971. Reptation of a polymer chain in presence of fixed obstacles. J. Chem. Phys. 55, 572. Depalle, B., Qin, Z., Shefelbine, S.J., Buehler, M.J., 2015. Influence of cross-link structure, density and mechanical properties in the mesoscale deformation mechanisms of collagen fibrils. J. Mech. Behav. Biomed. Mater. 52, 1–13. Di Giuseppe, E., Funiciello, F., Corbi, F., Ranalli, G., Mojoli, G., 2009. Gelatins as rock analogs: a systematic study of their rheological and physical properties. Tectonophysics 473, 391–403. Di Paola, M., Zingales, M., 2012. Exact mechanical models of fractional hereditary materials. J. Rheol. 56, 983–1004. Ewoldt, R.H., McKinley, G.H., 2007. Creep ringing in rheometry or how to deal with oftdiscarded data in step stress tests. Rheol. Bullet. 76. Faber, T.J., Jaishankar, A., McKinley, G.H., 2017. Describing the firmness, springiness and rubberiness of food gels using fractional calculus. Part I: theoretical framework. Food Hydrocolloids 62, 311–324. Fang, C.Q., Sun, H.Y., Gu, J.P., 2015. Application of fractional calculus methods to viscoelastic response of amorphous shape memory polymers. J. Mech. 31, 427–432. Gautieri, A., Vesentini, S., Redaelli, A., Buehler, M.J., 2012. Viscoelastic properties of model segments of collagen molecules. Matrix Biol. 31, 141–149. Goudoulas, T.B., Germann, N., 2016. Viscoelastic properties of polyacrylamide solutions from creep ringing data. J. Rheol. 60, 491–502. Heymans, N., Bauwens, J.C., 1994. Fractal rheological models and fractional differential equations for viscoelastic behavior. Rheol. Acta 33, 210–219. Higgs, P.G., Ross-Murphy, S.B., 1990. Creep measurements on gelatin gels. Int. J. Biol. Macromol. 12, 233–240. Holder, A.J., Badiei, N., Hawkins, K., Wright, C., Williams, P.R., Curtis, D.J., 2018. Control of collagen gel mechanical properties through manipulation of gelation conditions near the sol-gel transition. Soft Matter 14, 574–580. Holt, B., Tripathi, A., Morgan, J., 2008. Viscoelastic response of human skin to low magnitude physiologically relevant shear. J. Biomech. 41, 2689–2695. Jachowicz, J., McMullen, R., Prettypaul, D., 2007. Indentometric analysis of in vivo skin and comparison with artificial skin models. Skin Res. Technol. 13, 299–309. Jin, R.H., Cui, Y.C., Chen, H.J., Zhang, Z.Z., Weng, T.T., Xia, S.Z., Yu, M.R., Zhang, W., Shao, J.M., Yang, M., Han, C.M., Wang, X.G., 2021. Three-dimensional bioprinting of a full-thickness functional skin model using acellular dermal matrix and gelatin methacrylamide bioink. Acta Biomater. 131, 248–261. Joly-Duhamel, C., Hellio, D., Ajdari, A., Djabourov, M., 2002. All gelatin networks: 2. The master curve for elasticity. Langmuir 18, 7158–7166. Lee, K.H., 2000. Tissue-engineered human living skin substitutes: development and clinical application. Yonsei Med. J. 41, 774–779. Leyva-Mendivil, M.F., Lengiewicz, J., Page, A., Bressloff, N.W., Limbert, G., 2017. Skin microstructure is a key contributor to its friction behaviour. Tribol. Lett. 65. Li, X.Y., Tian, X.G., 2021. Fractional order thermo-viscoelastic theory of biological tissue with dual phase lag heat conduction model. Appl. Math. Model. 95, 612–622. Liu, C.Y., He, J.S., van Ruymbeke, E., Keunings, R., Bailly, C., 2006. Evaluation of different methods for the determination of the plateau modulus and the entanglement molecular weight. Polymer 47, 4461–4479. Long, J.M., Xiao, R., Chen, W., 2018. Fractional viscoelastic models with non-singular kernels. Mech. Mater. 127, 55–64. Macosko, C.W., 1996. Rheology: Principles, Measurements, and Applications. Wiley. Mahiuddin, M., Godhani, D., Feng, L.B., Liu, F.W., Langrish, T., Karim, M.A., 2020. Application of Caputo fractional rheological model to determine the viscoelastic and mechanical properties of fruit and vegetables. Postharvest Biol. Technol. 163. Mainardi, F., 2010. Fractional Calculus and Waves in Linear Viscoelasticity: an Introduction to Mathematical Models. Imperial College Press, London. Mansbridge, J., 2002. Tissue-engineered skin substitutes. Expet Opin. Biol. Ther. 2, 25–34. Metzner, A.B., White, J.L., Denn, M.M., 1966. Constitutive equations for viscoelastic fluids for short deformation periods and for rapidly changing flows - significance of deborah number. AIChE J. 12, 863–&. Morales-Hurtado, M., Zeng, X., Gonzalez-Rodriguez, P., Ten Elshof, J.E., van der Heide, E., 2015. A new water absorbable mechanical Epidermal skin equivalent: the combination of hydrophobic PDMS and hydrophilic PVA hydrogel. J. Mech. Behav. Biomed. Mater. 46, 305–317. Normand, V., Ravey, J.C., 1997. Dynamic study of gelatin gels by creep measurements. Rheol. Acta 36, 610–617. R. Mouˇ cka et al.