Keywords

1 Introduction

Over the last decades, huge progress has been made in discovering the properties and functionalities of the various constituents of cells [2]. In the new field of synthetic biology [7], this knowledge is used to produce some desired macromolecules by altering or replacing the DNA. Within the scope of the European Horizon 2020 project ACDC – Artificial Cells with Distributed Cores to Decipher Protein Function, we follow a different approach: We use droplets containing some fluid and being surrounded by another fluid as simplified models of cell-like structures. Droplets touching each other can form bilayers with pores through which chemicals can be exchanged. Thus, an agglomeration of droplets filled with various chemicals can form a bilayer network, allowing a gradual chemical reaction process. Our aim is to create a specific agglomeration of droplets in order to generate a desired gradual chemical reaction process for some purpose, like the generation of some desired macromolecule [17, 18]. In order to reach this goal, we need to study the dynamics and outcome of agglomeration processes [20], but also to develop a chemical compiler [6, 23], creating a plan of how to design an experiment leading to an appropriate agglomeration for a bilayer network for the desired gradual reaction processes.

Besides, we are also interested in the origin of life [15]. One of the open questions in this field is the huge difference between the time needed for a totally random evolutionary process, starting from the primordial soup and finally leading to life of larger complexity, up to human beings, and the age of the earth, which is too small to allow such an entirely random evolutionary process performed at a constant rate, such that faster rates in the past have to be assumed [10]. But if we replace the primordial soup at least partially with an agglomeration of droplets [27], the production of the first complex molecules of life might become likelier and thus faster, as the compartimentalization of semi-stable molecules within small droplets instead of within the overall primordial soup might stabilize these educts [11], thus increase the probability for the production of such complex molecules [15], and in turn shorten the time for some initial steps of the evolutionary process.

For both of these problems we want to deal with, the basic question arises which molecules could be created by an agglomeration of droplets and the resulting bilayer network. And vice versa, one could ask which agglomeration was needed to produce a specific desired molecule. But here an additional problem arises: Even if a bilayer network for the creation of that molecule can be determined theoretically, still the question remains whether this bilayer network can indeed be realized with an agglomeration of droplets filled with various chemicals. Geometric restrictions may form the most prominent obstacles towards the realizations of such bilayer networks.

In this paper, we restrict ourselves to the case that the droplets can be modeled as hard spheres, i.e., we consider the case of large surface tensions, such that the droplets have in general a spherical shape, and of small bilayers, such that the droplets only touch slightly each other and stay spherical if a bilayer is formed. Furthermore, we mostly restrict ourselves to the case of all spheres exhibiting the same radius value R. In the following sections, we present two examples of geometric restrictions of great importance for our underlying problems, show how one of these problems can be solved analytically, develop a heuristic optimization algorithm for finding (quasi) optimum solutions for the other problem, and present results obtained with that algorithm.

2 The Buttercup Problem

2.1 Description

Fig. 1.
figure 1

Solutions for the buttercup problem for \(N=2\), \(N=3\), \(N=4\), and \(N=5\), (top row from left to right) unit spheres on a ring (printed in red), with all of them touching one further unit sphere (printed in blue), with a glossy sphere encapsulating them and touching all of them. For \(N=6\), such a configuration is already impossible for a finite radius of the surrounding hull, such that we print here the configuration with the smallest hull sphere encapsulating a ring of 6 unit spheres (bottom). (Color figure online)

Full size image

The first problem we want to consider here has been named the buttercup problem (see Fig. 1): N unit spheres (We use the term unit sphere throughout this paper for a sphere with radius 1. But all results shown for these unit spheres can be easily transferred to other systems with spheres having an other common radius value R.) shall be placed around an additional \(N+1\)th center unit sphere in their midst in the way

  • that each of the N unit spheres touch the center unit sphere in their midst,

  • that the N unit spheres form a ring, i.e., each of them touches only its left and its right neighbor besides touching the center sphere and no further connections are formed, and

  • that the whole configuration is placed within a surrounding hull sphere of radius \(\mathcal{R}\), with \(\mathcal{R}\) being chosen in a way that the overall configuration is stabilized without further ado, i.e., the surrounding hull touches each of the N spheres on the ring and the center sphere.

From these conditions, it follows that the midpoints of the N unit spheres on a ring lie in one plane. This structure is reminiscent to a buttercup flower with its petals and center. It is important for our problem of generating specific agglomerations for some desired purposes, as the center sphere can serve as a control sphere which governs e.g. the opening and closing of pores in the bilayers between the petal spheres.

2.2 Analytic Solution

Now we have to determine the configurations for various N. Without restriction, we place the midpoints of the N petal spheres in the xy-plane on a circle with radius \(\varrho \), such that the midpoint of the petal sphere No. k has the coordinates \((\varrho \cos (2\pi k/N),\varrho \sin (2\pi k/N),0)\). For symmetry reasons, the \(N+1\)th center sphere has the coordinates \((0,0,z_{N+1})\) and the midpoint of the surrounding hull, which has a radius value of \(\mathcal{R}\), lies at (0, 0, h). First we determine the radius \(\varrho \) with the condition

$$\begin{aligned} (\varrho \cos (2\pi /N) - \varrho )^2 + (\varrho \sin (2\pi /N))^2 = 2^2, \end{aligned}$$
(1)

which resolves to

$$\begin{aligned} \varrho =\sqrt{\frac{\displaystyle 2}{\displaystyle 1-\cos (2\pi /N)}}. \end{aligned}$$
(2)

There are of course also some more elegant ways to determine \(\varrho \) for some specific values of N. For example, for the case \(N=5\), one can make use of the properties of the golden cut and finds \(\varrho =\frac{1}{5}\sqrt{50+10\sqrt{5}}\). In the second step, we determine \(z_{N+1}\) with the condition

$$\begin{aligned} z_{N+1}^2=4-\varrho ^2. \end{aligned}$$
(3)

So far, the locations of the \(N+1\) spheres have already been determined. Now we have to choose h and \(\mathcal{R}\) in a way that the outer hull touches all \(N+1\) spheres, which is only the case if there is a common distance value \(\mathcal{D}\) between the midpoint of the surrounding hull and the midpoints of the \(N+1\) spheres. For determining the two parameters h and \(\mathcal{R}\), we need two conditions:

$$\begin{aligned} h^2+\varrho ^2 & = & \mathcal{D}^2 \\z_{N+1}-h & = & \mathcal{D} \\\end{aligned}$$
(4)

Please note that there are two cases: If the midpoints of the \(N+1\)th sphere and of the surrounding hull lie on opposite sides of the xy-plane, then h is negative, as in the case \(N=5\). If they lie on the same side, then h is positive, as for \(N\le 3\). By eliminating \(\mathcal{D}\) from Eq. (4), we get

$$\begin{aligned} h=\frac{\displaystyle z_{N+1}^2 - \varrho ^2}{\displaystyle 2 z_{N+1}}. \end{aligned}$$
(5)

The radius \(\mathcal{R}\) can then easily be determined as

$$\begin{aligned} \mathcal{R} = z_{N+1}-h+1. \end{aligned}$$
(6)
Table 1. Approximate numeric values of the parameters (see text) of the solutions to the buttercup problem for various numbers N of unit spheres in a ring. If the spheres have a common radius value R different from 1, each of the numbers provided has to be multiplied by R.
Full size table

Table 1 provides the parameters of the buttercup problem for \(N=2\) to \(N=6\). Please note that the configuration for \(N=6\) in Fig. 1 violates the condition that the surrounding hull must touch all spheres, including the center sphere. This is only possible in the limit \(\mathcal{R}\rightarrow \infty \). Instead the picture for \(N=6\) depicts the case of the smallest surrounding sphere. But this configuration is not stable, as the center sphere can move freely in the z-direction. The configurations for \(N=2, 3,\) and 4 are identical with the densest packings of 3, 4, and 5 spheres in a sphere [26]. In the configuration for \(N=4\), even a further unit sphere can be added symmetric to the \(N+1\)th sphere at the other side of the xy-plane without overlaps and without the need to enlarge \(\mathcal{R}\). Furthermore note that no buttercup configurations can be created for \(N>6\) (see also the configurations for the related problem of packing circles in a circle in [24]), if not at least one of the conditions mentioned above is abandoned or the radius of the center sphere enlarged.

3 The Kissing Number Problem

3.1 Description

Fig. 2.
figure 2

The kissing number problem in three dimensions

Full size image

The second problem we want to study in this paper is the kissing number problem [3, 4]. It can be stated as follows: find the maximum number N of unit spheres touching a further unit sphere in their midst without overlaps in D dimensions [16]. This problem, which has been studied for centuries, can be trivially solved in one and two dimensions. In two dimensions, \(N=6\) unit disks can be placed around a center disk. In one dimension, the kissing number is given by \(N=2\). But in three dimensions, the exact value remained unknown for a long time. Newton and Gregory debated about this problem, with Newton claiming that the kissing number in three dimensions was \(N=12\), whereas Gregory believed it was 13. Only in the 1950s, it was proved that Newton had been right [22]. A configuration of an arrangement of spheres for this kissing number problem in three dimensions is shown in Fig. 2, with the center sphere printed in blue and the surrounding spheres printed in red. A glassy orb is drawn around the configuration to better show that the surrounding spheres indeed touch the center sphere in their midst. The kissing number problem is also considered in higher dimensions, but only for some dimensions like \(D=4\) (\(N=24\)), \(D=8\) (\(N=240\)), and \(D=24\) (\(N=196560\)), exact values of N have been determined. For most dimensions, only lower and upper bounds to the kissing number are known. A list of currently known values can be found in [25]. They indicate that the kissing number N grows exponentially with the dimension D. In the next step, we want to develop an optimization algorithm for the kissing number problem.

3.2 A Heuristic Optimization Algorithm for the Kissing Number Problem

Optimization algorithms can be divided in two classes, namely exact mathematical algorithms, which return one solution and guarantee that this solution is indeed optimal, and heuristic methods, which create one or several configurations of good quality, hopefuly being even (quasi) optimum. These heuristic methods can be subdivided in two subclasses, one of them being construction heuristics, which start at a tabula rasa and gradually create a solution for the overall problem, and the other one being iterative improvement heuristics, which usually start at a randomly chosen configuration and iteratively apply changes to the configuration in order to gradually increase the quality [21]. Mostly, one follows the local search [1] path and applies only small changes which do not alter a configuration very much. The simplest improvement heuristic is the greedy algorithm: starting from a randomly chosen initial configuration \(\sigma _0\), it performs a series of moves \(\sigma _i\rightarrow \sigma _{i+1}\) by successively applying small randomly chosen changes to the configuration. A move is accepted with probability

$$\begin{aligned} p(\sigma _i\rightarrow \sigma _{i+1}) = \left\{ \begin{array}{ll} 1 &{} \quad \text {if}\, \varDelta \mathcal{H} \le 0 \\ 0 &{} \quad \text {otherwise} \end{array}\right. \end{aligned}$$
(7)

with the difference \(\varDelta \mathcal{H}=\mathcal{H}(\sigma _{i+1})-\mathcal{H}(\sigma _i)\) between the cost function values of the current configuration \(\sigma _i\) and the tentative new configuration \(\sigma _{i+1}\). In case of rejection, one sets \(\sigma _{i+1} := \sigma _i\).

As the greedy algorithm does not accept any deteriorations, it often gets soon stuck at local minima of mediocre quality. In contrast, simulated annealing [9], which uses the Metropolis criterion [12]

$$\begin{aligned} p(\sigma _i\rightarrow \sigma _{i+1}) = \left\{ \begin{array}{ll} 1 &{} \quad \text {if}\, \varDelta \mathcal{H} \le 0 \\ \exp (-\varDelta \mathcal{H}/T) &{} \quad \text {otherwise} \end{array}\right. \end{aligned}$$
(8)

with the temperature T as acceptance criterion, also accepts deteriorations with some probability, which is the smaller the larger the deterioration \(\varDelta \mathcal{H}\) and the smaller the temperature T is. Its deterministic variant [13], which is also called threshold accepting [5], uses the transition probability

$$\begin{aligned} p(\sigma _i\rightarrow \sigma _{i+1}) = \left\{ \begin{array}{ll} 1 &{}\quad \text {if}\, \varDelta \mathcal{H} \le T \\ 0 &{} \quad \text {otherwise} \end{array}\right. \end{aligned}$$
(9)

and thus accepts all improvements and all deteriorations up to a proposed threshold value T. During the optimization run, T is gradually decreased from a large initial value, at which the system to be optimized is in a high-energetic unordered regime, towards a small vanishing value, at which the system gradually freezes in a low-energetic ordered solution.

When applying threshold accepting to a specific optimization problem like the kissing number problem, one thus has to define

  • a routine for creating a random initial configuration,

  • one or several move routines, which alter a configuration slightly, and

  • a function for measuring the cost function value of a configuration.

In our optimization apprach, we use threshold accepting, as its acceptance criterion is better at avoiding small overlaps at the end of a simulation than the Metropolis criterion. We work with a fixed number N of spheres, which shall be placed in a way that all of them touch the center sphere and that they do not overlap with each other. If a feasible configuration meeting these constraints is found for a specific value of N, we increase N by 1 and perform a new optimization run with this incremented value of N. We iterate this approach until no feasible configuration can be found anymore for a specific \(N_{\textrm{max}}\), such that the heuristic algorithm returns that the kissing number has the value \(N_{\textrm{max}}-1\).

For the initialization process, we could randomly place the N spheres anywhere. However, we want to restrict the search space to configurations in which the midpoints of the N spheres lie on a virtual sphere with radius 2, as this condition must hold for any solution of the problem. For all initial coordinates \(x_i(j), j=1,\dots ,D\) of sphere i, we first randomly select a value from the interval \([-2,2]\), such that the midpoints of the spheres are already placed in a D-dimensional cube of side length 4, which is centered around the origin, at the beginning. Then we calculate the distance

$$\begin{aligned} r_i=\sqrt{\sum _{j=1}^D x_i(j)^2} \end{aligned}$$
(10)

of the midpoint of sphere i to the origin and renormalize the coordinates of the midpoint of sphere i according to

$$\begin{aligned} x_i(j) := 2 \times x_i(j) / r_i, \end{aligned}$$
(11)

thus placing the midpoint of sphere i on the virtual sphere of radius 2. Therefore, the condition that each of the N spheres has to touch the sphere at the center is already fulfilled.

Furthermore, we implemented three move routines in which one of the N spheres is randomly chosen and transferred to a tentative new location:

  • In the first of these three move routines, the new location of the chosen sphere is determined as in the initialization process, i.e., first each of its coordinates are randomly chosen in a way that its midpoint lies in a cube of side length 4 centered around the origin, then they are renormalized using Eq. 11, i.e., in the way that the midpoint lies again on the virtual sphere of radius 2. Thus, the midpoint sphere can lie at an entirely new location on this virtual sphere.

  • In the second move routine, the tentative new location is created by adding randomly chosen values from the interval \([-1,1]\) to the current coordinates of sphere i, such that the midpoint of sphere i in the tentative new configuration first lies in a cube of side length 2 centered around the current position. Then the coordinates are renormalized as above, such that the midpoint lies on the virtual sphere of radius 2 again. We chose this maximum value of 1 for shifting a coordinate, as this maximum value proved to provide superior results for a formerly studied multidisperse disk packing problem [14, 19].

  • The third move routine is almost identical to the second move routine, except that the interval from which the values to be added to the coordinates are randomly chosen is not \([-1,1]\) but \([-T,T]\). With decreasing threshold value T, the amount by which the coordinates can be changed at a maximum decreases during the optimization run, thus increasing the likelihood for randomly selecting small changes with larger acceptance probability.

Each of these three move routines ensures that the N spheres touch the center sphere as required, however, configurations can contain overlaps between pairs of these surrounding spheres.

This brings us to the last point of our considerations of how to apply threshold accepting to the kissing number problem. For the application of an optimization algorithm like threshold accepting, a cost function \(\mathcal{H}\) is needed which is to be minimized. However, there is no such cost function for the kissing number problem in our approach with a fixed number N of surrounding spheres a priori, one can only state whether a configuration is feasible because all constraints are fulfilled or whether it is not feasible because at least one pair of surrounding spheres overlaps. However, during the optimization run in the search of a feasible configuration, we can allow to accept non-feasible configurations as well, but introduce penalties for the overlaps between pairs of spheres, which we sum up and thus get a cost function value. Therefore, we can define the cost function \(\mathcal{H}\) as

$$\begin{aligned} \mathcal{H}(\sigma )= \sum _{i,j=1 \atop i<j}^N \varDelta (i,j) \varTheta (\varDelta (i,j)) \end{aligned}$$
(12)

with

$$\begin{aligned} \varDelta (i,j) = 2-\sqrt{\sum _{k=1}^D (x_i(k)-x_j(k))^2} \end{aligned}$$
(13)

and the Heaviside step function

$$\begin{aligned} \varTheta (a)=\left\{ \begin{array}{ll} 1 &{} \quad \text {if}\, a\ge 0 \\ 0 &{}\quad \text {otherwise}. \end{array}\right. \end{aligned}$$
(14)

When a tentative new configuration is put to the decision whether to accept or reject it, not its overall cost function value has to be calculated. When determining \(\varDelta \mathcal{H}\) for an algorithm using local search moves as in our approach, it is much faster only to sum up the differences between those addends which have changed between the current and the tentative new configuration, i.e., one has to check only for the sphere which is to be displaced whether there are overlaps at its current position and whether there are at the tentative new position.

The threshold T is gradually decreased from an initial value of 10 by a cooling factor of 0.999. In each of the 10000 temperature steps, one million MCS are performed. An optimization run can be prematurely stopped when a feasible configuration without overlaps is reached.

3.3 Computational Results for the Kissing Number Problem

Fig. 3.
figure 3

Two solutions each of our optimization program for the kissing number problem in \(D=4\) dimensions with \(N=24\) surrounding spheres (left) and in \(D=6\) dimensions with \(N=72\) surrounding spheres (right): Only the midpoints of the N surrounding spheres are displayed at \((x_i(1),x_i(2))\). If surrounding spheres touch each other, their midpoints are connected with edges. The central sphere is omitted for visibility reasons.

Full size image

We applied the threshold accepting algorithm as described above and got the solution shown in Fig. 2 for the kissing number problem in three dimensions. But the kissing number problem does not provide much of a challenge in three dimensions, as the available space on the surface of the central sphere is almost sufficient for a thirteenth sphere. In order to show the effectiveness of our heuristic algorithm, we thus applied the algorithm to the kissing number problem also in higher dimensions. In several optimization runs, we achieved feasible solutions for the optimum number of \(N=24\) surrounding spheres in \(D=4\) dimensions and for \(N=72\) surrounding spheres in \(D=6\) dimensions, four of which are shown in Fig. 3. We also tried to find a feasible solution with \(N=73\) surrounding spheres in \(D=6\) dimensions, but were unsuccessful. The best solutions we achieved were configurations with a cost function value of 2, i.e., in these solutions, the 73rd surrounding sphere is placed exactly where the 72nd surrounding sphere lies. Although having no exact proof by using threshold accepting, we thus assume that the value of \(N=72\) is not only the lower bound to the true kissing number in six dimensions but that it truely is the kissing number in six dimensions.

3.4 Computational Results for a Bidisperse Kissing Number Problem

Fig. 4.
figure 4

Solutions for a bidisperse kissing number problem with a center sphere with larger radius R: \(N=28\) unit spheres can be placed around this larger center sphere without overlaps for \(R=2\), \(N=52\) for \(R=3\), \(N=83\) for \(R=4\), and \(N=120\) for \(R=5\) (from left to right).

Full size image

Finally, we apply our optimization algorithm to a bidisperse variant of the kissing number problem in which the center sphere has no longer a radius of \(R=1\) but a larger radius value, while the surrounding spheres are still unit spheres, i.e., still have a radius of 1. For this problem, we have to alter our optimization algorithm simply in the way that we have to ensure that the midpoints of the surrounding spheres lie on a virtual sphere of radius \(R+1\), therefore Eq. (11) has to be replaced by

$$\begin{aligned} x_i(j) := (R+1) \times x_i(j) / r_i. \end{aligned}$$
(15)

Otherwise, the optimization algorithm remains unchanged. We consider this altered problem only in three dimensions. Figure 4 displays the best solutions found for various radius values. The configurations are drawn as in Fig. 2. The kissing number N increases with increasing radius R of the center sphere. Due to the small problem sizes studied, we cannot generally state according to which laws N will furthermore increase if we move on to even larger values of R. However, if considering that the surface of the center sphere is given by \(4 \pi R^2\), we can construct an upper bound function to the kissing number N, which increases \(\propto R^2\), such that the kissing number itself cannot increase stronger than with \(\propto R^2\).

4 Conclusion and Outlook

In this paper, we considered two problems which impose restrictions to our aim of generating bilayer networks in agglomerations of droplets for the purpose of producing specific macromolecules. In the buttercup problem, spheres have to be placed like the petals of a buttercup around a center sphere while fulfilling some constraints. In the kissing number problem, the largest number of spheres being able to touch a further sphere in their midst without overlaps has to be determined. We provided analytic results for the first problem and obtained computational results from a heuristic optimization algorithm we had developed for the second problem, achieving feasible configurations with optimum kissing numbers in 3 and 4 dimensions and with the lower bound in 6 dimensions. Furthermore, we provided results for a newly introduced bidisperse variant of the kissing number problem. We will continue studying these problems and their variants with additional constraints. We will also consider other such problems related to our work, like the problem of sticky spheres [8], in which agglomerations of spheres have to be found in a way that the number of bilayers between them is maximum.