## Abstract

Microstructures responsible for temporary arresting of contacting surfaces are widely distributed on surfaces in different organisms. Recent morphological studies show that these structures have different density of outgrowths and not ideal distribution pattern on both complementary parts of the contact. One can suggest that this difference is optimized by natural selection to get stronger mechanical arrest within the system. In this paper, we simulate such a system numerically, both in the frames of continuous contact and discrete dynamical models to prove this hypothesis and elucidate other aspects of optimization of such mechanical adhesive systems.

## 1. Introduction

Specialized pairs of microstructure-covered fields of arthropod cuticle often function for mechanical adhesion (attachment) between different body parts within an organism [1–5]. These systems are often called *probabilistic fasteners* due to the fact that they are composed of cuticular micro-outgrowths of different origin [6] and because the size, shape and density of the outgrowths on both surfaces do not correspond exactly to each other. Thus, mechanical interactions between the micro-outgrowths of both surfaces take place without precise positioning of both surfaces [1]. Probabilistic fasteners prevent sliding of surfaces, when in contact. Mechanical adhesion (attachment), in this case, must depend on the size, shape, distribution density and the material properties of outgrowths (called elements in this paper), and is fast, precise and reversible. Individual elements are not necessarily hook-shaped (figure 1). The mechanism of attachment in such systems is definitely different from the typical hook-and-loop Velcro principle, but has been only marginally studied for the case of probabilistic fasteners with parabolic contact elements [5].

Because this kind of attachment structure occurs most widely among insects: in the unguitractor plate [8], intersegmental fixators of leg joints [7], synchronizing mechanisms of contralateral legs in cicada [9], dragonfly head-arresting system [10], and especially in wing-attachment devices [2,3,11,12], it is of the importance for the functioning of numerous biomechanical systems of insects. The major similarity of these devices is that co-opted fields of cuticular outgrowths are present on two separate parts of the body. Sometimes, as in the case of the forewing fixators of beetles, more than just one pair (up to eight pairs) of co-opted structures are involved in the proper functioning of the mechanical system. In such cases, each pair favours movement of the closed wings in one preferred direction. The head arrester of adult dragonflies immobilizes the head in some critical moments when it may experience strong mechanical disturbances [10]. Different shapes of outgrowths in functionally corresponding fields have been previously described [3]: (i) hook-like or cone-shaped or parabolic elements (outgrowths) and mushroom-like elements, (ii) cone-shaped or parabolic elements on both surfaces, (iii) clavate elements on both surfaces and (iv) plate-like elements on both surfaces [5].

## 2. Some structural properties observed in biological systems

As mentioned in §1, surface outgrowths of probabilistic fasteners vary in different biological systems and even between corresponding fields within the same system. However, in the pair of corresponding surfaces, one could expect some constraints in size and distribution of co-opted areas of outgrowths. In the dragonfly head arrester, the size of the outgrowths is usually larger in large species, and their density is lower, but the density of outgrowths on co-opted fields of one species is surprisingly never the same (figure 1). Similar observation has been made for the co-opted microtrichia fields of the forewing fixators of the beetles [2]. In the example of the dragonfly arrester, the density is higher on the neck sclerite than on the rear surface of the head [7]. The difference in the microtrichia density on the corresponding structures is more distinct in the larger species [7]. An increase in size of single microtrichia in larger species may be explained by possible correlation with an increasing cell size in larger species. It is more difficult to explain the difference in the outgrowth density on the corresponding structures. In addition, the distribution of the outgrowths shows a roughly hexagonal pattern, which is generally accepted to represent the highest package density of structures, but this pattern is far from an ideal one. It remains unknown whether these features have any functional significance on the function of co-opted fields.

Understanding the role of dimensions and distribution of outgrowths of co-opted fields is crucial for understanding the mechanical behaviour of such systems in contact. To the best of our knowledge, there is only the model of frictional behaviour of two single parabolic pins in sliding contact [5], but there are no models predicting the behaviour of arrays of outgrowths in such systems. This study was undertaken to establish simple models that realistically describe behaviour of co-opted fields of probabilistic fasteners in contact. Especially, we were interested in understanding the role of distribution of outgrowths on corresponding fields. To the best of our knowledge, the main question raised in this paper has not previously been discussed in the literature: how does the difference in densities and in the type of distribution of single elements influence the contact forces generated by the system of co-opted fields?

## 3. Continuous contact model

First, we analysed contact properties of a continuous numerical model. The conceptual structure of the model is shown in figure 2. We constructed numerically two flexible contacting surfaces with regular, or partially randomized, arrays of the outgrowths having different lengths, periods of distribution, symmetry and widths of the asperities on both upper and lower surfaces. For a definiteness and numerical simplicity, we simulate the outgrowths by the arrays of Gaussians
3.1*k* = 1, 2, … , *N _{x}*;

*j*= 1, 2, … ,

*N*

_{y}with numerically regulated height *G _{k}*

_{,j}, width and individual positions of their centres {

*x*,

_{k}*y*}. By manipulating these parameters and by varying distances between outgrowths of upper and lower surfaces, we were able to simulate different combinations of forms of outgrowths on both contacting surfaces.

_{j}It should be noted that the majority of the structures observed in biological systems possess more or less randomized hexagonal ordering [5,10,12]. Such a type of structure is well expected intuitively, as a self-organized densely packed array of relatively rigid outgrowths repelling one another by their cores. However, in many specific cases in real biological systems, these structures have a lower degree of symmetry and rhombic or striped quasi-periodic ordering rather than the ideal hexagonal one. In any case, one can generally suggest that the real ordering as a result of natural selection is some sort of compromise between self-organized dense packing of outgrowths and the specific functional requirements of the particular system; for example, its arresting efficiency, generated friction forces in contact and probably other functional properties, which can demand a symmetry specific for the required function. This hypothesis has not been proven rigorously in an experiment before, in spite of some experimental studies showing that this type of structure is generally very promising for the purposes of biomimetics [5,13–16]. Unfortunately, it is rather difficult to approach this problem experimentally because of the difficulty of designing structures with a low grade of symmetry for comparison with structures with a high grade of the symmetry. This is why we decided to build and analyse a numerical model, which allows simulation of how particular ordering of the outgrowths satisfies different prescribed functional demands.

We started with analysis of the contact properties of two surfaces with an ideal hexagonal distribution of their outgrowths:3.2where arrays of the coordinates {*x*_{k}, *y*_{j}} belong to the regular hexagonal lattice {*x*_{0k}, *y*_{0j}} with all constant Gaussian parameters for each of the two (upper and lower) surfaces: *G*^{u} = *G*^{0u} = const., *G*^{d} = *G*^{0d}, Δ^{u} = Δ^{0u}, Δ* ^{d}* = Δ

^{0d}. Despite the perfect symmetry of each surface, we allow different periods of outgrowths on upper and lower surfaces

*λ*

^{0u},

*λ*

^{0d}. Other fixed parameters of the surfaces can also be different:

*λ*

^{0u}≠

*λ*

^{0d},

*G*

^{0u}≠

*G*

^{0d}, Δ

^{0u}≠ Δ

^{0d}. This corresponds to our biological observations [2,7].

Below, after the study of the perfect symmetry case, we shall modify the surfaces by randomizing some of the mentioned parameters, or in the general case, all of them3.3

Here, *ζ*_{x,k} and are *δ*-correlated random numbers, which are independent for every surface, position and coordinate {*x*,*y*}:3.4

The intensity of the randomization can easily be regulated by choosing amplitudes of the random noise *ζ*. It is obvious that in the limit corresponding surfaces degenerate to the perfectly ordered ones, as well as in the opposite limit the system becomes absolutely disordered. This freedom in choice of the constants allowed us to model different combinations of weakly and strongly disordered surfaces as well as to study various contacts between completely ordered and differently randomized surfaces. Typical realizations of ordered and intermediately randomized surfaces are shown in figure 2 ((*a*, less randomized) and (*b*, more randomized), respectively).

When the surfaces contact one another, they deform proportionally to their elasticity *k*^{u,d}. If both are made of the same material, their elasticities coincide and they deform equally. In a numerical simulation, it means that, when converging surfaces *Z*^{u}(*x*, *y*) and *Z*^{d}(*x*, *y*) formally intersect one another, in the regions of intersection they form some equilibrium contact surface, which coincides with the mean value: . The deformations of both surfaces cause local elastic forces, which in this case are equal to . These forces correspond to the load acting on both contacting surfaces. Even at an unknown coefficient *k*, the load can be characterized by the deviation of the surfaces *δZ*(*x*,*y*) from their equilibrium positions. Nevertheless, from the observations of biological structures, one can estimate typical deformations *δZ*(*x*, *y*) for dragonfly species depicted in figure 1 in the interval 10–20 µm and corresponding elastic forces *f*^{elastic}(*x*, *y*) = *k δZ*(

*x*,

*y*) of about 10–20 µN at the typical spring constants of insect micro-outgrowths of 1 N m

^{−1}[17]. This estimation will be useful to validate the proposed model in further experimental studies.

If both surfaces represent ideal hexagonal lattices, one can expect *δZ*(*x*, *y*) to depend on a relation between their periods *λ*^{u,d} and mutual orientation between the surfaces (angle *φ* between their axes of symmetry). It is clear that the minima of *δZ*(*x*,*y*) should appear at some special orientations *φ*, when asperities of upper and lower surfaces fall between one another. It is especially predictable when the periods *λ*^{u,d} are commensurate or simply equal (*λ*^{u} = *λ*^{d}). The opposite situation is expected for the *strongest incommensurate* relation *λ*^{u}/*λ*^{d} = *τ*, where is the so-called *golden ratio*. For example, in the dragonfly head arrester, this relationship ranges from 1.5 to 4.0. Thus, it is interesting to note that this relation *λ*^{u}/*λ*^{d} = *τ* between the periods of outgrowths on both surfaces is quite close to some previously found relations in biological systems [7]. It means that the hypothesis that some incommensurateness between the periods *λ*^{u,d} is biologically important seems to be very plausible. To prove it, we have performed calculation of angle-depending configurations for different *λ*^{u,d}. As expected, at some angles *φ*, regular patterns of *δZ*(*x*, *y*; *φ*) were observed for the commensurate cases (especially for *λ*^{u} = *λ*^{d}). However, at the majority of arbitrary angles *φ*, projections of both lattices do not coincide, and rather complex contact patterns appear.

From a biological point of view, such an integral characteristic as the mean absolute value of *δZ*(*x*, *y*; *φ*) averaged over spatial coordinates should be important. We calculated it for different values of *λ*^{u,d}. The typical angle dependence of normalized to the maximum values for the two most important and representative cases *λ*^{u} = *λ*^{d} and *λ*^{u}/*λ*^{d} = *τ* are shown in figure 3*a*. It is seen directly from figure 3 that despite naive expectations, the curves in both commensurate and incommensurate cases have their own fine structure of minima corresponding to some special orientations *φ*.

According to the general definition [18], to estimate the local impact on friction between two rough surfaces, one also needs an absolute value of the inclination (or mathematical gradient) of the contact surface . Being integrated over the whole system at different mutual orientations *φ* and normalized per unit area, this gives the angle-dependent mean gradient . The results of the integration for two representative cases *λ*^{u} = *λ*^{d} and *λ*^{u}/*λ*^{d} = *τ* are shown in figure 3*b*.

Local impact to friction is proportional to the product of both values: local gradient and deviation determining local load, which is proportional to the elastic force *f*^{elastic}(*x*, *y*) = *kδZ*(*x*, *y*). The mean friction coefficient *μ* is given by a sum of these products normalized to the total load3.5and at constant *k*, it does not depend on the elasticity: . As we found, for all relations *λ*^{u}/*λ*^{d}, this practically important combination is much less dependent on the orientation *φ* than on solitary values of both and , and what is even more important, there is a big difference between commensurate and incommensurate cases. This is seen directly from figure 3*c*, where the thin line, representing the case *λ*^{u} = *λ*^{d}, has a well pronounced deep minimum, but the bold curve corresponding to *λ*^{u}/*λ*^{d} = *τ* is practically flat (or has small maxima at some angles).

As mentioned in §2, the distribution of outgrowths on biological surfaces is not perfectly hexagonal, but partially randomized. The same calculations as for the cases of regular distributions were repeated for the randomized hexagonal lattices , with conserved equal amplitudes of the random noises , … {*x → y*; *u → d*}. It was found that even for quite intensive randomization , the curves for and still have well pronounced minima. However, physically meaningful friction coefficients for both commensurate *λ*^{u} = *λ*^{d} and incommensurate *λ*^{u}/*λ*^{d} = *τ* cases demonstrate relatively small fluctuations. These fluctuations are related only to the limited size of the numerically simulated systems, so each individual calculated line varies from one realization to another.

All the results obtained for the values , , and friction coefficients *μ* of randomized system, are summarized in figure 4*a–c*. It is seen directly from figure 4*c* that the friction coefficient in the incommensurate case *λ*^{u}/*λ*^{d} = *τ* is definitely higher than in the commensurate case *λ*^{u} = *λ*^{d}.

## 4. Discrete model and dynamic simulations

Unfortunately, our continuous model may only partially reflect the real situation. Besides, in many cases, it also becomes time consuming to calculate dependencies of the system properties on the varied parameters, especially in the case of the dynamics simulation. It would be useful to create an alternative model which will allow the same system to be studied by numerical simulation of Brownian dynamics of simple mutually interacting subsystems. Below, we have developed and studied such a discrete model.

The conceptual structure of the discrete model is presented in figure 5. We simulated two contacting surfaces of probabilistic fasteners by interacting arrays of outgrowths, which are attached by one of their ends to the upper or to the lower flat surfaces. The outgrowths are allowed to rotate around the attachment roots with some flexibility and interact one with another according to simple rules, defined by corresponding interaction potentials. These properties of the model are qualitatively represented in figure 5, where two arrays of the outgrowths are depicted by two different kinds of the lines with attachment roots shown by the black and open circles, respectively. Owing to the interactions one with another and their flexibility, the outgrowths appear instantly bent as is shown in figure 5. As in the above continuous model, we can distribute the roots of the outgrowths differently. In particular, they can be ordered into a perfect or partially randomized hexagonal ‘dense packed’ lattice, or be placed in some other order.

By default, the outgrowths are vertically oriented at equilibrium *α*_{0} = *π*/2. The simplest way to simulate the deviations of the outgrowths from perfectly vertical positions is to use an effective elastic interaction *f ^{φ}* =

*B*(

*π*/2 −

*α*), where

*α*is the instant inclination angle. When the deviations are relatively small, the difference between

*α*and

*π*/2 is proportional to the distance between the projections of instant positions of the outgrowth ends {

*r**} and their root coordinates . Here, the discrete index*

_{j}*j*= 1, 2, …

*N*runs over all

*N*outgrowths and is the two-dimensional vector . Thus, in the following, we simulated elasticity of the outgrowths by the simple elastic relation4.1with an effective elastic constant

*k*

_{eff}.

Real biological outgrowths consist of sufficiently stiff polymers and in the nearest approximation can be modelled by rigid outgrowths, which can rotate around their root positions, but practically do not bend in their middle segments. Owing to the interaction one with another, they are repulsed by their hard cores and cannot penetrate one into another closer than to a characteristic distance corresponding to their widths. The widths of the outgrowths *w* can vary from one system to another, and can also be different for upper *w ^{u}* and lower

*w*surfaces. In real biological systems, these widths

^{d}*w*

^{u}^{,d}normally correlate with the density of outgrowths, or in other words with the periodicity

*λ*

^{u,d}of the structure . In the following, we estimated this choice of parameters in the interaction between the outgrowths.

The simplest way to simulate numerically the repulsion between the outgrowths is to use Gaussian interaction potentials4.2where *k* = 1, 2, … , *N _{x}*,

*j*= 1, 2, … ,

*N*; and are intensities and decays of the interactions, and indices

_{y}*a*,

*b*can independently take the values

*a*=

*u*,

*a*=

*d*,

*b*=

*u*and

*b*=

*d*, corresponding to upper and lower surface lattices. To perform a numerical experiment, the centre of mass

**of the upper lattice was pulled by an external elastic force with velocity . If the system has hexagonal (or any other symmetry), the angle**

*R**β*is calculated starting from one of the coordinates

*y*coinciding with any of the symmetry axes.

This model has been used to perform a number of experimental runs at different pulling angles *β*, for different symmetries and mutual orientation *φ* of upper and lower lattices, amplitudes of the randomization of the outgrowth distributions , and their rotational stiffness *k ^{a,b}*. These simulations extend and generally confirm the results obtained in the study of our continuous model. In figure 6, we summarize the most important results for the representative general orientation

*φ*=

*π*/5 and pulling direction

*β*= 3

*π*/10 obtained for two strongly different stiffness constants equal for both surfaces,

*k*=

^{a}*k*= 10 and

^{b}*k*=

^{a}*k*= 1. The simulations for all mentioned sets of the parameters were done at varied relations

^{b}*λ*

^{u}/

*λ*

^{d}between the periods of the perfect or randomized lattices in the interval 1 <

*λ*

^{u}/

*λ*

^{d}< 2.5.

The numerical experiment was organized as follows. At different *λ*^{u}/*λ*^{d}, the external force pulls the upper plate from its static initial position, while after some transient process, it starts to move from stationary. We record a maximum friction force which, according to the definition, determines the static friction threshold, and show its dependence on the relation *λ*^{u}/*λ*^{d} in figure 6.

It is seen directly from figure 6 that for relatively stiff outgrowths the results for the frictional behaviour of the discrete dynamic model are in a very good agreement with those of the continuous model. The function *F*_{max}(*λ*^{u}/*λ*^{d}) has its maximum around *λ*^{u}/*λ*^{d} ≈ *τ*. The vertical dashed line in the plot indicates the golden ratio *λ*^{u}/*λ*^{d} = *τ*. With an increased degree of randomization, the curve *F*_{max}(*λ*^{u}/*λ*^{d}) generally shifts to higher values. In particular, the dependence *F*_{max}(*λ*^{u}/*λ*^{d}), calculated at , is shown in figure 6 (the curve is marked by the black diamonds).

Good agreement of this behaviour and the above results of our continuous model probably arises from proper selection of the outgrowth rigidity *k ^{a}*

^{,b}= 10, allowing the outgrowths to remain almost vertical in the course of motion. This should be expected, because the continuous model completely ignored horizontal motion of the outgrowths, which were represented by Gaussian maxima with fixed centres.

It is interesting that in the system with softer outgrowths, friction behaviour is not much affected. For *k*^{a,b} = 1, both curves *F*_{max}(*λ*^{u}/*λ*^{d}) obtained in cases of perfect and randomized distributions (figure 6) are slightly shifted to higher values of *F*_{max} at relatively small *λ*^{u}/*λ*^{d} ≤ 1.5 values and to lower values of *F*_{max} for larger *λ*^{u}/*λ*^{d} ≥ 2 values. Hence an optimal relationship between the periods of outgrowths of co-opted surfaces, corresponding to the maximum of *F*_{max}(*λ*^{u}/*λ*^{d}), moves towards smaller values of *λ*^{u}/*λ*^{d} ≈ 1.45. Perhaps it happens because at extremely low stiffness of outgrowths it is easier for the fibres from one surface to adapt from the very beginning to the positions between the outgrowths of another surface. The same adaptation reduces a difference between the curves for perfect and randomized cases, because almost freely rotating ends of the outgrowths can leave their ideal lattice positions more easily.

## 5. Biological and biomimetic significance of obtained results

A previous analytical model of a pair of sliding parabolic outgrowths in contact showed that the pull-off force which can be generated by probabilistic fasteners comprising individual outgrowths, is of the same order of magnitude as the loading force needed to bring such arrays of individual outgrowths into contact [5]. Furthermore, that analytical model is valid only for the small deformations at the initial phase of the contact formation between sliding outgrowths. A full analysis considering outgrowth deformation, even for a single pair of outgrowths, can only be made numerically. However, in reality, the forces will strongly depend on the collective behaviour of numerous outgrowths, and will therefore be strongly dependent on their spatial distribution on co-opted lattices. This hypothesis can only be tested numerically.

It has been previously shown that in biological systems the density of outgrowths correlates to the length, width and especially the distance between them, so that the longest and widest outgrowths are usually sparsely distributed within the field [2]. However, there is only a weak correlation between length and width of outgrowths: longer outgrowths are not always wider. Distances between tips of longer and wider outgrowths are larger. The interplay of these parameters may result in different behaviour of contact element assemblages. This paper studied two numerical models predicting the behaviour of the probabilistic fasteners with different density of single outgrowths and different regularity of their distributions in co-opted lattices. The main finding of our virtual experiments is that there is an optimal relationship of densities on co-opted lattices for generation of the strongest contact force. Another important result is that the optimal distribution of the outgrowths should be not absolutely perfect (grid-like), but slightly randomized. Both these features are observed in various biological probabilistic fastener systems capable of fixing body parts as a result of friction enhancement.

Our finding that co-opted lattices within the same organism have different distribution of individual outgrowths was initially rather surprising for us [7], because the organism can definitely keep similar density of outgrowths on both surfaces by controlling the size of the cells responsible for development of such outgrowths. Furthermore, by keeping the same cell distribution on both surfaces, a highly ordered type of outgrowth distribution can be quite easily generated as, for example, in the case of the tarsal adhesive hairs in the fly. Numerical simulations, presented here, provide a strong argument for the functional significance of both certain difference in outgrowth density on co-opted lattices and certain imperfections in outgrowth distribution on each lattice in biological probabilistic fasteners.

The Velcro principle was transferred to industrial applications from biological observations. Recent patent databases contain a large number of ideas dealing with applications of existing hook-and-loop fasteners. They are widely used in the textile industry, medicine, sports, transport, etc. However, the vast majority of these applications use the same types of hook-like tapes. Indeed, these systems generate very high attachment forces, if compared with the applied load. Their function is quite long-term attachment, whereas probabilistic systems, discussed in this paper, might be of importance for designing short-term fast attaching and detaching devices, as previously described for macroscopic cases [5] and for microscopic systems [13–16]. For the purpose of optimization of such technical devices, knowledge concerning the optimal relationship between distributions of single outgrowths on co-opted lattices in biological systems, and the numerical study presented here might be of great importance.

## Funding statement

This study was partly supported by the SPP 1420 priority programme of the German Science Foundation (DFG) ‘Biomimetic Materials Research: Functionality by Hierarchical Structuring of Materials' (project GO 995/9-2) to S.N.G.

## Acknowledgements

We thank V. Kastner (Tübingen, Germany) for linguistic corrections of the manuscript.

## Footnotes

One contribution of 19 to a discussion meeting issue ‘Cell adhesion century: culture breakthrough’.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.