Journal of The Royal Society Interface
You have accessResearch articles

Functional significance of graded properties of insect cuticle supported by an evolutionary analysis

M. Jafarpour

M. Jafarpour

Institute of Zoology, Functional Morphology and Biomechanics, Kiel University, Kiel, Germany

Faculty of Mechanical Engineering, University of Guilan, Rasht, Iran

Google Scholar

Find this author on PubMed

Sh. Eshghi

Sh. Eshghi

Institute of Zoology, Functional Morphology and Biomechanics, Kiel University, Kiel, Germany

Google Scholar

Find this author on PubMed

A. Darvizeh

A. Darvizeh

Faculty of Mechanical Engineering, University of Guilan, Rasht, Iran

Google Scholar

Find this author on PubMed

S. Gorb

S. Gorb

Institute of Zoology, Functional Morphology and Biomechanics, Kiel University, Kiel, Germany

Google Scholar

Find this author on PubMed

H. Rajabi

H. Rajabi

Institute of Zoology, Functional Morphology and Biomechanics, Kiel University, Kiel, Germany

[email protected]

Google Scholar

Find this author on PubMed


    The exoskeleton of nearly all insects consists of a flexible core and a stiff shell. The transition between these two is often characterized by a gradual change in the stiffness. However, the functional significance of this stiffness gradient is unknown. Here by combining finite-element analysis and multi-objective optimization, we simulated the mechanical response of about 3000 unique gradients of the elastic modulus to normal contacts. We showed that materials with exponential gradients of the elastic modulus could achieve an optimal balance between the load-bearing capacity and resilience. This is very similar to the elastic modulus gradient observed in insect cuticle and, therefore, suggests cuticle adaptations to applied mechanical stresses; this is likely to facilitate the function of insect cuticle as a protective barrier. Our results further indicate that the relative thickness of compositionally different regions in insect cuticle is similar to the optimal estimation. We expect our findings to inform the design of engineered materials with improved mechanical performance.

    1. Introduction

    Many biological materials have graded properties. Gradients are often produced through changes in material composition, density and microstructure. Several studies have highlighted the crucial role of graded properties in the functionality of biological systems. Examples are as diverse as improved adaptability (e.g. in adhesive systems of beetles [1]), controllability (e.g. in antennae of stick insects [2]), deformability (e.g. in wings of flying insects [3]), interfacial bonding (e.g. in soft tissue-to-bone interfaces [4]), damping capacity (e.g. in peel of pumelo fruits [5]), buckling strength (e.g. in plant stems [6]), impact tolerance (e.g. in pericarp of the green coconuts [7]) and resilience (e.g. in spider fangs [8]).

    A striking example of the effect of graded properties on the functionality of a biological system is found in insect cuticle. The cuticle is a biological composite that forms the exoskeleton of all insect species. It is the primary interface between insects and environment, and particularly plays a major role in protecting insects against imposed mechanical stresses [9]. Taking into account the key function of cuticle as a mechanical barrier, it should be able to carry applied loads (i.e. being load-bearing) and quickly recover from them (i.e. being resilient). The latter, i.e. resilience, is especially very important, because it represents the capacity of a material to withstand forces without being permanently deformed [10]. Although a variety of morphological and structural strategies may serve to improve the two properties of the cuticle, the use of material gradients is likely to be a key one [11,12].

    The presence of material gradients has already been reported for the cuticle of different insects. A few examples are: wings of dragonflies and damselflies [1315], attachment pads and claws of grasshoppers, flies and dragonfly larvae [1618], tibiae of stick insects [19], intromittent organs of cassidine beetles [2022], gula of fruit chafer beetles [23], ovipositor of parasitic wasps [24] and mouthparts of damselfly larvae [25]. Mechanical tests showed substantial changes in the material properties, in particular the elastic modulus, of some of these examples [22,23,26]. The results not only indicated variations in the elastic modulus over cuticle surface (i.e. within the same body part or different body parts), but also modulus variations over several orders of magnitude across cuticle thickness. Such variations are particularly important, because they suggest site-specific modifications that tune cuticle properties.

    Although graded materials evolved hundreds of millions of years ago, their use in ‘modern’ engineering does not exceed 40 years [27]. Inspired by biological materials, engineers developed the concept of functionally graded materials (FGMs) to develop materials with enhanced performance. Similar to their natural counterparts, FGMs are characterized by gradual changes in material properties, which are often caused by variations in the composition or structure. A few examples are graded metal–ceramic composites [28], graded glass–ceramic composites [29] and fibre-reinforced composites with graded fibre orientations [30]. Due to their improved mechanical resistance, such materials have found numerous applications in dental and medical implants, forming and cutting tools, and protective coatings to control damage by mechanical, or thermal, stresses [3133].

    Nowadays, modern manufacturing techniques provide satisfactory control over compositional and structural changes in FGMs. This, thereby, allows the development of materials with any arbitrary stiffness gradient. Despite this ability, the gradients in most existing FGMs can be described by only a few functions (i.e. the power-law, exponential and sigmoid) [3436]. This is partly because of our limited knowledge of the influence of material gradients on the mechanical behaviour of materials [37].

    The interest in understanding the influence of graded properties on two key properties of load-bearing capacity and resilience of materials prompted this study. Here, we use an optimization method and combine it with finite-element modelling to develop and analyse a rich set of graded models. The models are subjected to normal contacts to find those that can simultaneously withstand maximal loads with minimal permanent deformations. The reached optimal gradients are then compared with those observed in insect cuticle. Our results are expected to not only shed light on the evolutionary adaptations of cuticle as a protective coverage for insect body, but also inform the design of man-made FGMs with improved properties.

    2. Methods

    2.1. Development of a geometric model

    The finite-element software package ABAQUS/Standard v. 6.14 (Simulia, Providence, RI, USA) was used to develop a geometric model. The model was considered to be a cylinder, which is rotationally symmetric about the central axis (figure 1a). To decrease the computational time, the cylinder was reduced to a two-dimensional (2D) axisymmetric geometric model, as shown in figure 1b. The thickness of the model was set to be constant and equal to 100 µm, and its length was considered to be large enough, i.e. 650 µm, to prevent the influence of the boundary condition on the results. The model was meshed using the three-node axisymmetric elements (CAX3). These elements are general purpose and result in reasonably accurate solutions in a reduced computational time. Following a mesh sensitivity analysis, the size of the elements situated within a distance of 120 µm from the applied load was set to be 0.9 µm (table 1). The elements outside of this range had a length of 9 µm. Geometric nonlinearity was taken into account to capture the effect of large strains.

    Figure 1.

    Figure 1. Numerical simulation and optimization analysis. (a,b) Schematic diagram of the geometric model and indentation simulation. The geometric model (a) was reduced to an axisymmetric 2D model (b), which was fixed at the innermost layer and subjected to indentation with a spherical tip at the outermost surface. (c) Deformation and stress distribution in a representative model after loading (i) and after unloading (ii). (d) Distribution of optimal solutions along the Pareto front. (e) Distribution of the elastic modulus at 12 sites across the thickness of optimal models. The inset shows the enlarged view of the data for Sites 1–6. (f) NRMSE of six standard functions used to predict the optimal elastic modulus gradients. C-LIN, ‘continuous linear’ function; C-EXP, ‘continuous exponential’ function; C-LG, ‘continuous logistic’ function; DC-LIN, ‘discontinuous linear’ function; DC-EXP, ‘discontinuous exponential’ function; DC-CST, ‘discontinuous constant’ function. (g) Comparison between the mean optimal elastic modulus gradient (red line within the box in e) and six standard functions across the thickness of the geometric model. (h) Relative thickness of epi-, exo- and endocuticle in optimal models. EPI, epicuticle; ENDO, endocuticle; EXO, exocuticle. In all bar plots, the borders of the boxes indicate the 25th and 75th percentiles, a solid line (black/white) within them marks the median, a red line within them shows the mean and the whiskers (error bars) define the 10th and 90th percentiles.

    Table 1. Sensitivity of results to changes of modelling parameters.

    tested parameter quality/quantity penetration load (mN) relative difference (%) permanent deformation (µm) relative difference (%)
    boundary condition at the innermost layer of the models fixeda 38.11 2.48
    free slip 37.42 1.81 2.43 2.02
    friction between the indenter and models µ = 0 38.11 2.48
    µ = 1 39.42 3.44 2.34 5.65
    length of the models (µm) 550 38.05 2.48
    600 38.10 0.13 2.48 0
    650a 38.11 0.03 2.48 0
    700 38.12 0.03 2.48 0
    element size beneath the indentation site (µm), element size in the rest of the models (µm) 0.9, 9a 38.11 2.48
    0.8, 9 38.20 0.24 2.48 0
    0.7, 9 38.26 0.16 2.48 0
    0.9, 9a 38.11 2.48
    0.9, 8 38.10 0.03 2.48 0
    0.9, 7 38.10 0 2.48 0

    aConsidered in the current study.

    2.2. Assignment of material properties

    The elastic modulus of the geometric model was allowed to vary across the thickness over a large range from 2 MPa to 9 GPa. The two values correspond to the minimum and maximum elastic moduli of load-bearing insect cuticles, with protective functions, based on the most reliable experimental data [1,23,26,3840]. It is important to mention that some cuticles, such as maggot cuticle, may have much lower elastic moduli, in the range of a few kilopascals [41]. However, such pliant cuticles have only a minor load-bearing function that is the focus of this study; they mainly serve a locomotion function, according to [38]. To create a gradient of the elastic modulus, 12 sites that were always situated at an equal distance from each other were chosen across the thickness of the model (figure 1b). The elastic moduli at these sites were regarded as design variables and allowed to be varied, within the defined range of 2 MPa–9 GPa, by a genetic algorithm developed in this study (see §2.4). The elastic modulus at any site between the 12 sites was obtained by linear interpolation.

    To capture the deformation of our models, we used an elastic–perfectly plastic material model [42]. The onset of plasticity was predicted using the von Mises yield criterion, to account for plasticity that is likely to occur in the localized loading simulated here [42]. The yield strength was taken the same as that of insect cuticle, which is equal to 72 MPa [43]. The Poisson's ratio was also assumed to be constant and equal to 0.3, which is the same as that measured for many biological materials [44].

    2.3. Loading and boundary condition

    Models were fixed in all directions at their innermost layer (figure 1a). The outermost surface of the models was subjected to indentations with a spherical tip (r = 10 µm), to quantify the response of the models to normal contacts. A surface-to-surface contact pair was used to define the interaction between the models and the indenter at the contact site. A ‘hard contact’ relationship was used in the contact direction (i.e. normal direction) to prevent the penetration of the contact pairs into each other. The friction between the contact surfaces was assumed to be negligible. The indenter displacement was set to be equal to 7% of the total thickness of the models to avoid the effect of boundary conditions on the obtained results [42]. After reaching the maximum displacement, the models were unloaded to measure their permanent deformation. This was obtained by calculating the difference between the coordinates of a node located beneath the indenter tip on the model before and after loading. The reversible deformation of the models was measured by subtracting the permanent deformation of the models from the indenter displacement. All the simulations were considered to be semi-static and conducted using a general ‘static step’ [45].

    2.4. Optimization using genetic algorithm

    A genetic algorithm (GA), developed in MATLAB (R2012a, Mathworks, Natick, MA, USA), was used to perform a multi-objective optimization. The aim of the optimization was to obtain elastic modulus gradients that could enable our geometric model to simultaneously achieve (i) a maximal load carrying capacity and (ii) a maximum resilience. For this purpose, the load required to reach the defined indenter displacement (see §2.3) was used as a measure of the load-bearing capacity of the model and its reversible deformation was taken as a measure of the resilience. Taking into account that the resilience in this study was compared among models that were subjected to the same indentation displacement, a model with a higher resilience is the one that also underwent lesser permanent deformation. Therefore, the resilience in our study can also be used to estimate the resistance of the material to permanent deformations.

    The optimization started with 60 random elastic modulus gradients developed by the GA, as the initial population. No control was applied to these gradients, except that the values of the elastic modulus across their thickness should always be within the range defined in §2.2 (i.e. 2 MPa–9 GPa). For each elastic modulus gradient, the GA developed an ABAQUS input file, .inp, which contained all the information needed to develop and analyse a model including the number and location of nodes, connection of elements, material properties, type of the analysis, loading and boundary conditions, etc. A representative .inp file in provided in electronic supplementary material, data S1.

    Here the ‘whole arithmetic recombination’ and ‘uniform’ mutation were used as the crossover and mutation operators, respectively (for details, see [46]). Different rates of crossover and mutation were tested, and finally, they were set to be 0.6 and 0.2, respectively. These were the minimum reached rates that avoided local optima in a reasonable computational time. A ‘fitness proportionate selection’ method was used to choose individuals for the crossover and mutation. The stopping criterion was to reach a stable, well-distributed Pareto front. To test the accuracy of the results, the optimization process was repeated for multiple times. Electronic supplementary material, figure S1, presents the flowchart of the GA used in this study.

    To describe the profile of the optimal elastic modulus gradients, here we used two sets of mathematical functions, which have previously been used in finite-element models of insect cuticle [42]. The first set includes three continuous functions: (i) a ‘continuous linear’ (C-LIN), (ii) a ‘continuous exponential’ (C-EXP), and (iii) a ‘continuous logistic’ (C-LG) function. The second set includes three discontinuous functions: (i) a ‘discontinuous linear’ (DC-LIN), (ii) a ‘discontinuous exponential’ (DC-EXP), and (iii) a ‘discontinuous constant’ (DC-CST) function (table 2).

    Table 2. Mathematical functions used to describe the profile of the optimal elastic modulus gradients.

    function continuity/discontinuity equation determination of constants
    linear a continuous linear elastic modulus profile across the whole thickness (C-LIN) E(t) = At + B using the values of the elastic modulus at the bottom and top of each model
    a discontinuous linear elastic modulus profile across the thickness of individual layers (DC-LIN) using the values of the elastic modulus at the bottom and top of each individual layer
    exponential a continuous exponential elastic modulus profile across the whole thickness (C-EXP) E(t) = CeDt using the values of the elastic modulus at the bottom and top of each model
    a discontinuous exponential elastic modulus profile across the thickness of individual layers (DC-EXP) using the values of the elastic modulus at the bottom and top of each individual layer
    logistic a continuous logistic elastic modulus profile across the whole thickness (C-LG) E(t)=F+GF(1+HeI(tJ))1/k using the values of the elastic modulus at the bottom and top of each individual layer
    constant a discontinuous constant elastic modulus profile across the thickness of individual layers (DC-CST) E(t) = L using the mean values of the elastic modulus of individual layers

    3. Results

    In total, 2988 unique gradients of the elastic modulus were developed in this study. They were assigned to a geometric model and were subjected to indentations by a spherical tip (figure 1a). After each analysis, the load required to reach a constant indenter displacement (equal to 7% of the model thickness), as the measure of the load-bearing capacity, and the reversible deformation of the model, as the measure of the resilience, were recorded. Figure 1c shows the results of a representative simulation, which include the deformed shape of the model and the stress within it, at the end of the loading and unloading (electronic supplementary material, video S1). Our aim here was to find an elastic modulus gradient that could simultaneously maximize both load-bearing capacity and resilience of the model.

    Figure 1d represents the final spread of optimal solutions along the Pareto front (raw data are available in electronic supplementary material, data S2). The abscissa here is the normalized reversible deformation. The ordinate is the normalized load. Both the reversible deformation and the load were normalized using the min–max normalization method [47]. To visualize the optimal gradients, we plotted the distribution of the elastic moduli at 12 equally distributed sites across the thickness of the geometric model (figure 1e and inset). The results show an increasing trend in the elastic modulus from the innermost layer to the outermost surface of the model. The rate of changes of the elastic modulus is smaller at the deeper sites and notably increases at those closer to the surface.

    In order to examine the goodness of the fits, we calculated the normalized root mean square error (NRMSE) for two sets of continuous and discontinuous elastic modulus functions (see §2.4). As seen in figure 1f, the DC-LIN and DC-EXP functions provided best fits to the optimal gradients, compared with the others. The smallest NRMSE, however, was calculated for the DC-EXP function (electronic supplementary material, table S1).

    To better visualize the fit between the standard and optimal gradients, we compared them in figure 1g. Here, comparisons were made between the standard gradients and the mean values of the optimal gradients. Again, here the best fit with the optimal results was achieved by the DC-EXP gradient (electronic supplementary material, table S2). Although here the DC-LIN and C-EXP functions also exhibited acceptable fits, they slightly underestimated and overestimated the mean optimal gradient, respectively.

    Using insect cuticle as a ‘reference’, we estimated the relative thickness of epi-, exo- and endocuticle in models with optimal elastic modulus gradients (figure 1h). According to the most reliable data in the literature, the elastic modulus of the epicuticle is approximately 9 GPa [23] and those of the exo- and endocuticle range between 1 and 7.5 GPa [23,48] and 2 and 700 MPa [26,48,49], respectively. We used these data to estimate the thickness of the optimal models, through which the elastic modulus varies in the ranges mentioned above. We found that the elastic moduli of 8.0 ± 5.7%, 21.6 ± 6.2% and 70.4 ± 10.2% of the total thickness of the optimal models lie within the range of the moduli of the epi-, exo- and endocuticle, respectively. This finding, therefore, suggests that in an optimal cuticle, the three layers should comprise different fractions of the total thickness, as described above.

    4. Discussion

    4.1. Motivation for the use of an exponentially graded material

    Modern FGMs were initially invented, in the late 1980s, to minimize thermal stresses in components used in space planes [27]. However, due to their enhanced durability, they soon found use in a variety of low-temperature applications [31,50]. Our results suggest that FGMs with exponential elastic modulus gradients may have advantages over others; they can reach an optimal balance between load-bearing and resilience. This is an important finding, which can inform the design of impact-resistant materials. However, this also raises a question: why are materials with exponential elastic modulus gradients better than other graded materials?

    The answers to these questions lie in the fact that a compromise between the load-bearing and resilience in a material can be attained by combining two mutually exclusive properties of stiffness and compliance. Such a combination can be achieved in a material that consists of two layers: a stiff outer layer and a flexible inner layer. The stiff layer should be placed outside to reach a stiffer response at an applied displacement. On the one hand, the thickness of the stiff layer should be maximized to increase the material resistance to applied loads. On the other hand, the thickness of the flexible layer should be maximized to provide the material with a higher resilience. Based on our results, a compromise between the two properties is reached in a material that has a flexible layer notably thicker than stiff layer (approx. 9 times, figure 1h). However, it is also known that a discontinuous transition between two dissimilar layers can result in stress concentrations at their interface when subjected to external loads [51]. A graded transition in the composition can supress the stress at the critical interface of the two layers and, thereby, reduce the risk of plasticity [50]. All these together imply that an optimal balance between the load-bearing and resilience is reached in a graded material that consists of a comparatively thin stiff outer region and a thicker flexible region beneath that. The graded transition between the two regions should be smooth enough, to prevent stress concentrations, and sharp enough, to keep the elastic moduli of the stiff and flexible layers at their maximum and minimum values, respectively, while maintaining the optimal thickness ratio of the two layers. This condition is better satisfied by an exponential gradient of the elastic modulus than the others (figure 1g).

    4.2. Biological significance of graded properties in insect cuticle

    To protect insects against external mechanical stress, cuticle should be able to carry loads without being permanently deformed. To achieve these two purposes, a variety of strategies may be served. The simplest strategy is perhaps to increase cuticle thickness. A striking example is the thick cuticle of soldier ants of Carebara diversa, in which cuticle thickness exceeds 9 times that of workers [52]. Although a thick cuticle benefits soldiers for nest defence, it is costly to be produced and can restrict insect locomotion.

    A second strategy is to arrange chitin fibres along the direction in which the stress is maximum. This is because fibres are stronger in the axial direction. This strategy particularly works for compliant cuticles, in which the mechanical properties are governed by the properties of fibres, and not by those of the matrix. The arthrodial membrane cuticle of the African migratory locust, Locusta migratoria migratorioides, represents a good example. The membrane cuticle exhibits a strong anisotropy, so that in one orientation it is approximately 50 times more load-bearing than in the perpendicular orientation [38]. This anisotropy results from the preferred arrangement of the fibres in the membrane that are oriented around the abdomen [53]. This was suggested to enhance the deformability in the longitudinal axis of the insect, yet maintaining the load-bearing capacity orthogonal to the direction of extension. Orienting the majority of the fibres in a certain direction, however, can reduce the load-bearing capacity of the cuticle in the other directions.

    Sclerotization is the third strategy that enhances cuticle resistance against external loads. Although cuticles from different species show different degrees of sclerotization, the elytral cuticle of the weevil Pachyrhynchus sarcitis is an extreme example. Unlike the cuticle of all previously examined species, the weevil is equipped with a thoroughly sclerotized cuticle [12], which plays a significant role in protecting it from bites of lizard predators [54]. This strategy, however, comes with a disadvantage: a very low fracture resistance of the elytral cuticle [12]. The elytron of P. sarcitis catastrophically fails when the ultimate strength is reached.

    Here, we showed that the use of an exponential gradient of the elastic modulus is another strategy, which can provide an optimal compromise between load-bearing capacity and resilience. Surprisingly, this is similar to the elastic modulus gradient observed in the cuticle, which was shown to have an exponential profile [42]. Hence, we can conclude that the observed match between the optimal elastic modulus gradient and that of insect cuticle is a biomechanical adaptation to stresses imposed on the cuticle, which could permanently deform cuticle and impact insect survival. In comparison with the other mentioned strategies, the use of graded properties seems to be a less costly and more efficient solution to improve both load-bearing and resilience of insect cuticle. It is especially very important, because these two properties (i.e. load-bearing and resilience) cannot be easily combined in conventional engineering materials [55].

    To interpret our optimization results in the context of the classical definition of insect cuticle as a layered material, we subdivided the reached optimal exponential gradients into layers with different stiffness profiles. Our results suggested that an optimal gradient may consist of three distinct layers (figure 1h): a thick flexible inner layer, a medium-thick intermediate layer and a thin stiff outer layer. This is similar to our conventional understanding of the cuticle of most examined species [48]; insect cuticle is typically subdivided into three layers of endo-, exo- and epicuticle, from inside to outside, the thickness decreasing from the former to the latter. It should be noted that the ratio of the layers is not always constant and, in fact, may notably change during insect growth. In a developing adult locust, for instance, the cuticle thickness significantly changes over the first two to three weeks following the final moult [56,57]. During this period, the thickness of the exocuticle increases, in terms of stiffness, by the progressive sclerotization of the endocuticle that is rapidly deposited at an average rate of approximately 2 µm d−1 [57]. Therefore, our findings here may be applied to the relative thickness of the cuticular layers after the growth period.

    4.3. Generality of the results

    In this study, we employed an evolutionary algorithm to assess the influence of a wide variety of elastic modulus gradients on the mechanical response of materials to normal contacts. We used the most reliable literature data on the physical and mechanical properties of insect cuticle to provide our models with required inputs. This later enabled us to compare reached optimal gradients with those observed in insect cuticle. Considering the constraints imposed on the models, how much can we generalize the reached conclusions?

    To answer this question, we assessed the influence of the different inputs on the results. We altered some of the applied constraints, including the effects of the boundary condition at the innermost layer of the models, friction between the indenter and models, dimensions of the models and the element size. We found that the results were consistent and changing the mentioned parameters had only minor effects on the results (table 1). We also investigated the effects of a few other assumptions made in our study. These include the assumption of a constant Poisson's ratio and a constant yield strength. Although no direct experimental evidence is available on the Poisson's ratio of insect cuticle, our assumption of a constant Poisson's ratio of 0.3 fits that reported for a variety of biological materials including bone [44]. In addition, we do not expect a significant change in the results by the consideration of a varying Poisson's ratio. This is because a previous computational study showed that the response of materials to indentations is relatively less sensitive to variations in the Poisson's ratio than the variations in the elastic modulus [58].

    To test the sensitivity of the results to the gradients in the yield strength, we also analysed the variation of the stress across the depth beneath the indented surface in multiple randomly selected models in the Pareto front (see electronic supplementary material, figure S2, for a representative example). In all the examined models, the stress was the greatest right beneath the indenter, but suddenly dropped over a small distance from the surface and reached below the strength of the flexible cuticle. This suggests that the risk of plasticity is only high in regions right beneath the indentation site and, hence, the occurrence of plasticity is mainly influenced by the strength of the material in the uppermost parts of the model. Therefore, the results are unlikely to be very sensitive to variations in the yield strength. In-depth investigation of the strength gradients, however, would be an interesting area for future investigation. Taking into account the anisotropic properties of insect cuticle, it remains unclear how variations of the material properties in other directions, except the loading direction, can influence the trade-off between the load-bearing and resilience. We believe this to be an interesting area for future investigation.

    Finally, we asked whether the defined upper and lower limits of the elastic modulus can influence the reached conclusions. As seen in our results, the maximum and minimum values of the elastic modulus in our optimal models converged to the minimum and maximum allowed values. Hence, it is very likely that by changing the set upper and lower limits, the range of the elastic modulus in our graded optimal models will change too. However, we do not expect this to influence our conclusions regarding the role of the exponential elastic modulus gradient in obtaining an optimal balance between load-bearing and resilience. In fact, the exponential gradient is expected to be scaled to the new range of the elastic modulus.

    5. Conclusion

    Here using a combination of numerical simulations and multi-objective optimization analysis, we showed that:


    Materials with exponential elastic modulus gradients can achieve an optimal balance between the load-bearing capacity and resilience.


    The optimal elastic modulus gradient is very similar to the modulus gradient of insect cuticle.


    The similarity of the optimal elastic modulus gradient to that of insect cuticle suggests cuticle adaptation to imposed mechanical stresses.


    Not only the elastic modulus gradient, but also relative thickness of layers in insect cuticle is likely to be adapted to provide it with an improved mechanical performance.

    Data accessibility

    Supporting data are in the electronic supplementary material.

    Authors' contributions

    Conceptualization: M.J., S.E., S.G., A.D. and H.R.; data curation: M.J. and S.E.; formal analysis: M.J., S.E. and H.R.; investigation: M.J. and S.E.; methodology: M.J., S.E. and H.R.; project administration: A.D., S.G. and H.R.; resources: A.D.; software: M.J. and S.E.; supervision: H.R.; validation: M.J. and S.E.; visualization: M.J. and H.R.; writing—original draft preparation: H.R.; writing—review and editing: M.J., S.E., A.D., S.G. and HR.

    Competing interests

    We declare we have no competing interests.


    S.E. was supported by German Academic Exchange Service (DAAD). M.J. was supported by the Federal State Funding at Kiel University.


    We would like to thank Mr Vahid Nooraeifar (Ahrar Institute of Technology and Higher Education) for his helpful comments on our optimization analysis. We also appreciate the continuing support of Ahrar Institute of Technology & Higher Education and allowing access to their computing system.


    second email:

    Electronic supplementary material is available online at

    Published by the Royal Society. All rights reserved.