Numerical simulation of the pattern formation of the springtail cuticle nanostructures

1. Introduction

Springtails (Collembola) are soil-dwelling arthropods [1] with a complex hierarchically structured cuticle surface [2] with strong repelling properties against water, low-surface-tension liquids [3,4] and sticky secretions of predatory insects [5]. These properties favour the collembolan cuticle adaptation to living in soil habitats [2,6]. Nanoscopic, comb-like structures are primarily responsible for the robust surface-repellent properties, as previously demonstrated by polymer replication methods [7,8].

In general, different layers of arthropod cuticle have different chemical composition [9]. Both exocuticle and endocuticle consist of chitin fibres and a proteinous matrix, which could be more or less sclerotized. The difference between exocuticle and endocuticle is mainly in the fibre orientation and in the degree of sclerotization (endocuticle is less sclerotized). In epicuticle, four layers with corresponding composition could be distinguished: internal layer (polymerized lipoprotein stabilized by quinones [10]), cuticulin layer (quinone-tanned lipoproteins), wax layer (mainly long-chain saturated alcohols esterified with acids [11]), and cement layer (proteins and lipids stabilized by various polyphenolic substances [12]). However, the chemical nature of the collembolan cuticle comb-like structure was not known until the publication by Nickerl et al. [13], where the chemical composition and architecture of the cuticle of Tetrodontophora bielanensis was studied. The authors applied a stepwise removal method of the different cuticle layers and by this method separately analysed the chemical composition of different layers. It was shown that the endocuticle and the exocuticle in representatives of Collembola consist mostly of chitin lamella. Within the epicuticle, the comb-like nanostructures consist of glycine-rich structural proteins, and the wax layer consists of fatty acids (palmitic and stearic acid), fatty and steryl esters, terpenes (lycopene), steroids (cholesterol and desmosterol) and triglycerides. A regular distribution of pore channels within the chitin-containing deeper cuticle layers was also shown, which may be used for different functions, such as fluid transport from epidermal cells to the surface and respiration. The presence of pores with possible respiratory function underlines the importance of the water-repellent properties of the collembolan cuticle surface. The main result of their study was that the non-wetting property of the collembolan cuticle is mainly based on cuticle topography rather than on surface chemistry. In spite of having knowledge about the structure, chemistry and properties of a biological surface, it is difficult to judge the development of such elaborate hierarchical nanostructures as those of springtails. According to the dimensions of the structures, they are produced at the level of single cells, but even at the lower resolution of the scanning electron microscope (SEM), the boundaries between single cells are difficult to reveal, which may support the hypothesis that the structures appear due to the self-organization process, based on several densities of materials synthesized by underlying epidermal cells.

In material science, surface coatings with regularly arranged particles have been produced by colloidal lithography, a technique that uses the self-assembly of nanoparticles on a substrate surface [14,15]. However, the range of producible patterns and properties of these coatings are restricted, as well as their durability. Therefore, the study of biological self-organizing nanopatterns seems to be promising for generating new solutions [16,17] for industrial applications.

In the present paper, we studied the cuticle surface of the springtail Orchesella cincta using a SEM. The images were obtained on dry specimens mounted on SEM aluminium stubs using double-sided carbon-containing sticky tape, sputter-coated by gold-palladium (thickness 10 nm) and observed in a Hitachi S4800 SEM at 3 kV accelerating voltage (figure 1). The cuticle of O. cincta is covered by nanoscale structures, which form complex hexagonal comb-like patterns. The sizes of the comb elements were slightly different, so structures with five or seven vertices occurred. A comb element possesses a double wall circular/oval circumference. At sites where boundaries of three comb elements meet together, elevating triangle-like structures appear. That is why five to seven triangles are located on each boundary of a comb element (some just circle). Each triangle contains a little pore in its centre, figure 1b. Sometimes, this pore is plugged by the nanoscopical droplet of apparently amorphous solidified fluid. Additionally, within the circles, other slightly larger pores occasionally occur (figure 1).

• Download figure
• Open in new tab
• Download PowerPoint

The mechanisms of nanostructure assembly are still not clearly understood, but some models already exist. In particular, self-assembly processes for ommatidia gratings of insect eyes, snake skin nanodimples and arachnid cerotegument [16–18], which are rather simple nanostructures in comparison to the collembolan cuticle patterns (figure 1). Several different approaches were previously used in modelling. As the molecular/chemical background of the nanoscale pattern formation in springtails remains unknown, in order to gain a better understanding of the process of self-assembly and self-arrangement of surface nanostructures, we apply a theoretical approach in the present paper. We introduce a numerical model that allows us to study the effect of different interaction properties between different substances on the morphology of the eventual structure. This approach was already used for the mathematical problems related to the so-called Tammes problem, which is looking for the optimal packing of a given number of pores or particles on a sphere with maximized distance between them [19–21]. The aim was to alter the model parameters in a way that 3D structures similar to those observed in springtails could be numerically reproduced based on the principle of frozen kinetics, previously discussed in another paper [16]. Different sizes (polydispersity) disturb the optimal hexagonal distribution and generate other locally optimal distributions described in the framework of random packing [22]. Here, the optimal distribution corresponds to the minimum amount of energy. We wanted to test which minimal alterations are necessary to transform the structures within a possible evolutionary scenario, and likewise how colloids can be tuned to tailor 3D functionalized surfaces by colloidal lithography.

2. Numerical model

From the physical point of view, one can expect that the observed surface nanostructures appear in some kinetic process during which spatially distributed substances are produced from some nucleation centres on the surface, mutually interact, redistribute and solidify. Such a process can be observed in many different chemical and physical systems including such underlying processes in biological systems. That is why closely related patterns can appear convergently in very different systems. For example, similar patterns with hexagonal or honeycomb lattices can form during ordering of the superconducting vortices [23] or at magnetic ordering [24] in quasi two-dimensional systems. This is also possible during self-organization at growth of surface structures in biological systems [16,25].

In all the cases, the process involves an interaction of many spatially distributed densities and its direct simulation can cause extremely time-consuming calculations. Furthermore, the simplest numerical model describing the hexagonal comb-like pattern formation, according to our considerations, is suggested. Taking into account the complexity of the real pattern under consideration in our particular case, one can expect that to reproduce it for the self-organization process, we will need several (at least three) different kinds of interacting substances with densities ρ_k. They should correspond to the dark grey fields inside the honeycomb-like formations (ρ₁), to the middle grey boundaries between them (ρ₂) and additional density (ρ₃) which reproduces the light grey ‘triangles’ situated on the top of the intersections of the two previous ones (figure 1).

The main ‘heuristic’ hypothesis, which we plan to use in the present study, is that all of the densities are generated during a common kinetic process and redistribute in the space due to mutual interactions. The second hypothesis, strongly related to the previous one, is that the material of the boundaries (ρ₂) is deposited due to an expansion of the dark grey fields of (ρ₁). In other words, it is synthesized and at the same time ‘pushed’ by the expanding regions of (ρ₁). The sources of densities (pores) are distributed by a particles repulsion frozen kinetic process described in [16].

Common expectation of the kinetic studies tells us that if many individual sources of some density (ρ₁) repulse one another, they redistribute on a uniform surface to the patterns visually close to slightly disturbed hexagonal lattice. This is close to the real patterns observed on the springtail cuticle surface. Such a lattice appears as a kind of ‘frozen’ (or stopped) kinetic process [16,17]. Besides honeycomb cell domains, it normally contains some frozen ‘dislocations’ with five- and sevenfold symmetry. This is also very close to what we see here in our SEM images (figure 1).

Taking all the above into account, one can start from the hypothesis about a set of the naturally distributed repulsing sources of the first density (ρ₁), which grows from some initial value to a stationary configuration (figure 2, ρ₁₁, ρ₁₂). During growth and expansion, this density gradually creates and at the same time dislodges the second density (figure 2, ρ₂). Such self-consistent process can naturally form observed honeycomb-like boundaries, which are clearly seen in the SEM images (figure 1).

• Download figure
• Open in new tab
• Download PowerPoint

However, from the mathematical point of view, the following important problem appears when we try to realize such a process in a numerical experiment. If we use common (unique) spatial distribution of the density ρ₁ for all the growing nucleation centres, the domains easily collide without competition. As a result, the boundaries between them and corresponding walls of the second density ρ₂, even being formed at intermediate stages of kinetics, will gradually smear with time.

The only way to allow such a ‘mathematical catastrophe’, which we are able to propose, is to numerate directly all the sub-densities: ρ₁ as formally independent particular density ρ_1,j interacting with all other densities of the problem.

For example, if the number of sources N is equal to N = 100 (this is a typical number which we normally use for the calculations below), the number of the interacting densities of the first type (only!) also extends to N = 100, j = 1 : N. Total (and experimentally measured physical density) ρ₁ of the substance is given by the sum of all the local ones:

2.1

Numerical experiments confirm that, because of the mutual repulsion, the local densities ρ_1,k are strongly localized inside their own domains. As a result, the difference ρ₁ − ρ_1,k with fixed number k reproduces the total density from which a particular subsystem ρ_1,k is removed. The corresponding term g₁₁(ρ₁ − ρ_1k) being added below to the kinetic equations will describe an interaction between the given density ρ_1,k and all other densities with some intensity g₁₁.

According to our SEM images, the third density (ρ₃), corresponding to the light grey triangles in the images, seems to be the last one which enters into the process at its very late stage. Thus, let us forget about this density and write down the kinetic equations only for the densities ρ_1,k and ρ₂.

A general set of equations describing an ordering of these densities can be written as follows:

2.2

k = 1 .. N. Here, terms with the Laplacian operator c_jΔρ_j = c_j(∂²ρ_j/∂x² + ∂²ρ_j/∂y²) simulate surface tension with the strength controlled by the coefficient c_j. Surface tensions are different for different densities, but common for all variables ρ_1,k of the same nature. That is why we have the terms c₁Δρ_1k and c₂Δρ₂. The bigger the coefficient c_j, the larger the impact on the density gradients (mathematically corresponding to the surface tension) into the energy of the system, and the smoother density distributions will be energetically favourable. The analogous physical process could be the displacement of solution on the cuticle surface by a solution excreted through the cuticle pores. Mutual repulsion of ρ_1,k may be provided by the formation of a layer with high surface energy between interchanging solutions/fluids.

Production of the second density is simulated by the term , which means that this density is created by the expanding boundaries of all the domains . Mainly near the boundaries, the absolute value of gradient essentially deviates from zero and causes strong production of ρ₂. It is exactly the case, if ρ₂ for some reasons has higher concentration near to the outer boundary of ρ_1,k.

Standard combination ρ_j(1 − ρ_j)(ρ_j − 1/2) corresponds to the two-valley potential free energy which favours two equivalent equilibrium states and the barrier between them. In the trivial case, when all other densities and gradients are absent, each solitary density ρ_j tends to one of the two values of the uniform equilibrium: ρ_j = 0 or ρ_j = 1. This term could be alternatively assigned to dynamical processes with much shorter characteristic time than in equation (2.2). Other terms, which mix different densities with one another in the equations, namely −g₁₁(ρ₁ − ρ_1k) − g₁₂ρ₂ and −g₂₁ρ₁, describe repulsion between them, due to which every density tends to push out another one from the same position in space.

A known difficulty of 2D simulations with a large number of spatially distributed densities is that they are relatively time consuming and need a complex presentation of many consequent time-depending configurations by the surfaces or colour density-maps. It makes it difficult to extract and analyse reasonable information about space–time evolution of the system. It especially relates to the interacting densities, which can overlap one another. So, before introducing a 2D model, we will reduce the description for the beginning to a much simpler 1D model.

From the mathematical point of view, it means that we treat all the variables depending on the one coordinate x only and rewrite Laplacian and gradient operators in their one-dimensional forms: c_jΔρ_j = c_j∂²ρ_j/∂x² and , respectively. This space reduction gives us an opportunity to record space–time evolution of the interacting densities and present time-depending results statically. It allows us to visualize how the densities are generated, suppressed or overlapped by one another.

One-dimensional reduction also allows us to reduce a number of variables in equation (2.2). In 2D space, where expanding domains of the densities ρ_1k can grow in different directions, one needs a relatively large number of neighbouring domains of ρ_1k (normally around 5 or 6 domains) to lock the density ρ₂ between them. Rigid confinement in 1D space allows this to be done by only two domains ρ₁₁ and ρ₁₂, which unavoidably meet one another expanding along sole line x. It strongly simplifies calculations and allows us to do preliminary simulations much easier.

The result of the simulations according to equation (2.2) for all three variables ρ₁₁, ρ₁₂ and ρ₂ in time–space coordinates is shown in figure 2. Here, we can see domains ρ₁₁ and ρ₁₂ growing and expanding from left and right. Each one produces its own walls of ρ₂. They move and combine with the time into static double wall ρ₂ between two mutually stabilized domains of ρ₁₁ and ρ₁₂.

An important feature of the real structure in figure 1 is the existence of light triangles formed by a third material which is associated with the density ρ₃. They certainly formed in the places, where double walls of ρ₂ meet one another on 2D surface. It causes an additional problem of finding these positions numerically. For now, we postpone this question and study it later and note that the simplicity of the 1D model allows the third density ρ₃ to be included in this preliminary study. The basic equation for the third density can be written as follows:

2.3

Its structure reflects the general hypothesis that the density ρ₃ is created during the final stages of the process on the walls of the secondary density ρ₂. In other words, it is expected to be produced by non-local term in the equation (gradient ). In the normal case, if it is already produced and passed the energy barrier, the density tends to the equilibrium in the two-valley free energy potential represented in the equation by the combination ρ₃(1 − ρ₃)(ρ₃ − 1/2). The density has surface tension c₃Δρ₃ and interacts locally with all other densities −g₃₂ρ₂ + g₃₁ρ₁ in the same manner as it was in equations (2.2), so it has high affinity to ρ₂ and low affinity to ρ₁.

Surface tension, local interactions and the influence of the gradient are regulated by the relations between all the coefficients of the problem c₃, g₃₂, g₃₁, c₃₂. In particular, direct observation of the natural pattern shows that the triangles of the third density ρ₃ demonstrate negative curvature and their vortices are turned in the direction of the first density ρ₁. It is highly probable that this combination of the properties appears due to the low surface tension, strong repulsion with the second density ρ₂, attraction to the first one ρ₁ and a strong enough rate of production regulated by the coefficients g₃₁, c₃₂.

It allows us to make some assumptions about the relationship between phenomenological constants c₃, g₃₂, g₃₁, c₂₁ and make preliminary adjustments using simple the 1D approximation. Typical distribution of all the densities in asymptotic (stationary) configuration at fixed set of the parameters is shown in figure 3.

• Download figure
• Open in new tab
• Download PowerPoint

Let us note a fine structure of the dashed curve of ρ₃. Its sharp central maximum appeared mainly due to the competition between the strong influence of the gradient term and repulsion from the maximums of ρ₂ by itself. At the same time, ρ₃ still conserves relatively wide wings penetrating both domains of the first density ρ₃ due to their mutual attraction. This fine structure reflects, in the restricted 1D manner, specific distribution and orientation of its triangle-like density ρ₃, which will be seen on the more realistic 2D surface.

Now, we can return to the study of the 2D case. Numerical complexity of the problem motivates us to do an additional reduction of the numerical model. The main idea of the simplification is to substitute an original continuous problem by its combination with some preliminary applied discrete approach. In particular, it can be used to get a more or less final or, at least, almost stationary configuration. A well-known example of such an approach is the study of topological phase transition in quasi 2D superconducting systems [26], where extremely complex evolution of the superconducting vortices was substituted, to some extent, by a motion of charged particles. This approach was found to be extremely fruitful and the researchers were awarded the Nobel Prize in 2016. It is also applicable for the Tammes problem [19–21], where the study of kinetic ordering on a sphere was substituted by a calculation of equilibrium ordering of N equal charges.

In all these cases and in the used approach, it is important that one can finalize the discrete distribution of nucleation centres for the future densities and only then simulate kinetic evolution of the interacting continuous fields. In the present study, we organized a discrete part of the procedure as follows. Initially, we randomly distributed some particles (further nucleation centres for densities) on the limited rectangular planar surface.

In such a simulation, the most important parameters are the relationships between the area of the surface L_x × L_y, the number of the particles N, and the characteristic distance of the repulsion interaction R₀ between them. To get desirable ordering of the particles, their interaction can be chosen simply as a pure short range (exponential) repulsion:

We used the limited area L_x × L_y for our numerical simulation. However, in reality, the particles can move on a very large surface (in comparison to the mutual distance between them). From the mathematical point of view, to study such a quasi-infinite system we must apply periodic (toroidal) boundary conditions in both directions. We have to take into account that the original process represents damped kinetics, so dynamic equations here should involve damping as well: , where dissipative constant γ defines the characteristic time scale and γ⁻¹ can serve as a time unit for the simulation.

During the calculation routine, the particles dynamically rearrange with a velocity that gradually decays. One can control this slow-down by calculating the mean velocity of the particle and stop the procedure when their relative motion becomes negligible. In fact, originally randomly distributed particles could not reach real ground state in the toroidal space. But, at a long time run, they demonstrate some quite realistic distribution. As we see from the images of the real system, the same distribution actually appears in Nature, where the patterns are not perfect. In both numerical and real cases, various possible realizations are produced as a result of frozen kinetics.

One particular realization of the almost periodically (but not perfectly) ordered particles is shown in figure 4 by the open circles. Below we will use them as nucleation centres for the individual domains j = 1, 2, …, N of the first substance ρ_1j. One has to pay attention to the couple of positions marked by dark circles, where nucleation centres are doubled. It is done on purpose, in order to reproduce real cases, where they are doubled (figure 1). Such configurations appear when two sources of the substance ρ₁ occasionally appear closer at the beginning of the kinetic stage. In this case, the expanding domains of ρ₁ in this area are not able to form a double wall of ρ₂ between them and remain connected to the very end of the process. Taking account of such observed possibilities, we added these configurations to the simulation.

• Download figure
• Open in new tab
• Download PowerPoint

Before starting with the continuous model, we still have one more question to answer concerning the discrete model. Let us note that when double walls of the second density ρ₂ meet one another at places of intersection of domain boundaries, the third density ρ₃ starts generating at these positions. To simplify the calculations, it is also worth defining the positions of nucleation centres for ρ₃ in a discrete style. To do this, we created a new array of ‘particles’, which repulse one another and the already fixed particles that define the positions of the growing domains of ρ₁. After a long enough run, they reach stationary positions as far as possible from all the previous ones. These positions for the already calculated configuration of open circles are marked in figure 4 by black dots.

Now we have everything ready for the solution of our ultimate problem on 2D surface. To do this, we distribute a set of density nucleation centres in the preliminary found positions {x_k, y_k}, k = 1, 2, …, N and solve numerically connected kinetic equation (2.2) for ρ_1k and ρ₂ with 2D Laplacians c_jΔρ_j = c_j(∂²ρ_j/∂x² + ∂²ρ_j/∂y²) and gradients . When distributions of the densities ρ_1k and ρ₂ are attracted to some state close to the equilibrium, we insert them in the two-dimensional version of equation (2.3) and find the third density distribution ρ₃ associated with ρ_1k and ρ₂. The result of the solution for a particular realization of the nucleation centres is shown in figure 4. For the relatively large and small scales, it is presented in figures 5 and 6, respectively.

• Download figure
• Open in new tab
• Download PowerPoint

• Download figure
• Open in new tab
• Download PowerPoint

The present study revealed that the cuticle surface of the springtail Orchesella cincta, showing the characteristic collembolan ornamentations, consists of a complex arrangement of nanostructures interacting during their formation. Three densities used in our model might correspond to three main chemically different layers recently revealed in collembolan cuticle [13]. The thickest and deepest layer of the cuticle, so-called procuticle consisting of exo- and endocuticle, is chitin-rich. The specifically arranged pore channels [27,28], observed on the springtail cuticle surface and within its procuticle [13,29] enable material transport, and are presumably involved in the complex pattern formation by depositing fluids interacting with the surroundings and with each other and, after solidifying, contributing to the formation of the surface pattern. Previous studies have demonstrated that epicuticular structures contain proteins with high amounts of glycine, tyrosine and serine, the amino acid composition resembling that of fibroin, collagen or resilin [13]. The outermost layer representing a lipid mixture of fatty acids, wax esters and terpenes terminates the epicuticle and forms the protective barrier for the animal.

3. Conclusion

We have shown that it is possible to explain complex 3D nanoscale morphology found in springtails, by the interaction of three different materials that define their arrangement on the cuticle surface. These interactions depend on the chemical and physical properties of the fluids and presence/distribution of pore channels. Thus, according to our model, the complex 3D hierarchical surface structure of the collembolan cuticle with strong anti-wetting ability might be assembled in a very simple self-organizing way. From the biological point of view, it will be interesting to reveal the underlying genes producing the building blocks, and their alteration throughout evolution. For taxonomists, it would be interesting to find some deeper relationships between different pattern morphotypes and derive some evolutionary hypotheses. Our results may also assist in the engineering of colloids used for the functionalization of surfaces in order to gain different 3D structures with tailored wetting, mechanical and optical properties.

Data accessibility

This article does not contain any additional data.

Authors' contributions

A.E.F. and S.N.G. conceived the study and drafted the manuscript; A.E.F. performed numerical experiments; all authors discussed the model, interpreted the results and contributed to writing the final version of the manuscript.

Competing interests

We declare we have no competing interests.

Funding

This work was partly supported by the Georg Forster Research Award (Alexander von Humboldt Foundation, Germany) to A.E.F. The preparation of this paper was partly supported by the Leverhulme Trust (project CARBTRIB ‘Nanophenomena and functionality of modern carbon-based tribo-coatings’).

Footnotes