Mathematical modeling and computer simulation of the three-dimensional pattern formation of honeycombs

We present a mathematical model, a numerical scheme, and computer simulations of the three-dimensional pattern formation of a honeycomb structure by using the immersed boundary method. In our model, we assume that initially the honeycomb cells have a hollow hemisphere mounted by a hollow circular cylinder shape at their birth and there is force acting upon the entire side of the cell. The net force from the individual cells is a key factor in their transformation from a hollow hemisphere mounted by a hollow circular cylinder shape to a rounded rhombohedral surfaces mounted by a hexagonal cylinder shape. Numerical simulations of the proposed mathematical model equation produce the rounded rhombohedral surfaces mounted by a hexagonal cylinder patterns observed in honeybee colonies.

We present a mathematical model, a numerical scheme, and computer simulations of the threedimensional pattern formation of a honeycomb structure by using the immersed boundary method. in our model, we assume that initially the honeycomb cells have a hollow hemisphere mounted by a hollow circular cylinder shape at their birth and there is force acting upon the entire side of the cell. the net force from the individual cells is a key factor in their transformation from a hollow hemisphere mounted by a hollow circular cylinder shape to a rounded rhombohedral surfaces mounted by a hexagonal cylinder shape. numerical simulations of the proposed mathematical model equation produce the rounded rhombohedral surfaces mounted by a hexagonal cylinder patterns observed in honeybee colonies.
Although it is well known that the shape of the honeycomb cell is hexagonal, the three-dimensional structure of it is less known. Actual honeycombs consist of two layers of congruent cells, each one with a hexagonal opening and is capped by three rhombic faces 1 . Figure 1a represents a schematic illustration of four unit honeycomb cells.
Let us consider a hexagonal cylinder with side a and height h, see Fig. 1b, then volume is = . V ah 1 5 3 2 . If we cut the three shaded tetrahedra and put them together on the top of the hexagonal cylinder, then we have a three-dimensional geometry of a honeycomb cell as shown in Fig. 1c. The area of the surface with one open side a h x a a x 3 (2 ) 3 3 0 25 2 2 , where x is the height of the shaded tetrahedron. With a fixed volume, i.e., a and h are fixed, S has the maximum value at = x a/(2 2). Karihaloo et al. 2 reported that the fresh cells in a natural honeybee comb have a circular shape however quickly transform into the rounded hexagonal structure, while the comb is being built. Bird nests in nature are typically hemispherical in shape. In this study, it is proposed that rounded hexagonal and rhombic pattern formation in the honeybee comb is the result of the net force coming from the forces exerted outward by the cell boundaries. This work extends a previous study of honeycomb formation that only considered two dimensions 3 ; we examine the structure from a three-dimensional perspective.

Simulation Results
The proposed method was implemented using C language and the visualization of results was performed using MATLAB software (The MathWorks, Natick, MA, USA). For simulation, we assume the initial array of cells as shown in Fig. 2a, where each cell contacts other cells and has thin soft layer. We consider the initial unit cell as a hollow hemisphere mounted by a hollow cylinder as shown in Fig. 2b. The radius of the hemisphere is = R 1 with a center = x y z ( , , ) (0, 0, 1) and the height of the cylinder is = H 2. We generate the surface triangulation of the unit cell using DistMesh 4 , which is an open source MATLAB mesh generator, see Fig. 2c.
We investigate the temporal evolution of 50 cells by the proposed model. We denote the gap distance between the cells by g, see Fig. 2d. Figure 3 shows the snapshots at time = ∆ t t 0, 300 , and ∆t 700 . Here, we use = N 73 h 0 2496, ∆ = . t h 0 1 2 , and = .
g 0 1. Figure 3a,b are two different views. To examine the evolution of the cells in detail, we isolate four cells including the highlighted cell from the other cells as shown in   www.nature.com/scientificreports www.nature.com/scientificreports/ Fig. 3c. We can observe the cells transform from a circular cylinder to a rounded hexagonal cylinder and from a hemisphere to a rounded rhombohedral surface.
Let us consider the morphology evolution of the highlighted cell in more detail. Figure 4a-c show the different views of the snapshots of the highlighted cell with a best fitted reference frame at = ∆ t t 0, 120 , and ∆t 600 , respectively. The reference frame is the mathematically optimized shape. The result indicates that the highlighted cell is getting fit to the optimized shape, which is observed in nature.
So far, we performed computational tests with uniform boundary force as shown in Fig. 5a. As the final test, we investigate the temporal evolution of 50 cells with random boundary force as shown in Fig. 5b. Numerical parameters are the same in the previous test except time step, ∆t. Here, we use time step size as ∆ = t h 2 . The boundary force is randomly given in space and time, i.e., we randomly apply boundary force to 10% of cells. The numerical results are represented in Fig. 5c-e. Figure 5c,d are two different views and Fig. 5e is the four cells including the highlighted cell. Here, we can observe the cells transform from a circular cylinder to a rounded hexagonal cylinder and from a hemisphere to a rounded rhombohedral surface.

conclusions
We proposed a mathematical model, a numerical scheme, and computer simulations of the 3D pattern formation of a honeycomb structure using the IBM. In the proposed model, the cells having a hollow hemisphere mounted by a hollow circular cylinder shape at their birth becomes an optimized honeycomb structure by the net force from the outward forces of the individual cells. Computational experiments demonstrated the rounded rhombohedral surfaces mounted by hexagonal cylinder patterns observed in nature. In this paper, we focused on the mathematical modelling for the three-dimensional honeycomb formation and its numerical method. As the future work, it would be interesting to compare the theoretical results with those of experiments.

Method
In this paper, we propose a mathematical model and perform computational simulations of the three-dimensional honeycomb formation using the immersed boundary method (IBM) [5][6][7] . The basic mechanism of the three-dimensional model is similar to that of the two-dimensional model 3 and is as follows: We set the two layers of congruent cells and each cell consists of a hollow hemisphere mounted by a hollow cylinder.
First, we calculate forces acting on the entire surface of the individual cells. Second, we compute the net force from those forces. Third, we move the cell boundary points to new positions according to the net force. We repeat these steps until a specified number of iterations is reached. Now, let us describe the proposed scheme more precisely. We discretize a set of honeycomb cells by a triangular surface mesh are the lists of vertices and triangles of a honeycomb cell p, respectively.
We use the following evolution equation proposed in 3 : www.nature.com/scientificreports www.nature.com/scientificreports/ where α is a proportional constant and possess the dimension [time]/[mass]; and F is the net force resulting from the forces exerted outwardly by the individual cell boundaries such as attaching wax, movement of larva, and storing honey. Let a computational domain embedding the discrete honeycomb cells be partitioned in Cartesian geometry into a uniform mesh with mesh spacing h. The center of each grid cell, Ω ijk , is located at Here, N x , N y , and N z are the numbers of cells in the x-, y-, and z-directions, respectively. On the Cartesian grid, we applied the zero Dirichlet boundary condition. Let I X ( ) be the set of the index of the triangles containing a vertex X. Then, the unit normal vector at the vertex X is given by , and G q is the centroid of T q 8 , see Fig. 6. We define the boundary force per unit surface area σ = f N X ( ), p n p n where σ is a magnitude constant. To compute the net force density, we spread the boundary force over the nearby lattice points:

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.