Keywords

1 Introduction

Over the last decades, huge progress has been made in biochemistry. A large amount of knowledge about the constituents and the processes within a cell has been gathered [1]. Even a new research field, that of “synthetic biology”, has evolved [2], in which natural objects like the DNA in cells are purposedly altered or replaced in order to achieve some desired outcome, like producing some drug. Still, some questions remain unanswered so far, like one of the basic questions for the origin of life: Which constituent of a cell came first, the RNA [3] or the cell membrane [4]?

Fig. 1.
figure 1

Sketch of the initial and final states of the spatial rearrangement of droplets.

Full size image

In our approach, which we intend to follow within the European Horizon 2020 project ACDC, we do not consider fully equipped cells but the most simplified cell-like structures, being droplets comprised of some fluid and surrounded by another fluid. As an additional feature, we also allow droplets being contained within some outer hulls, playing the role membranes have for cells. Droplet generation, especially in the field of microfluidics, has been extensively studied over the past years [5,6,7,8] and has become an easy-to-use technology after the introduction of 3D printing technologies [9, 10]. A stream of fluid is broken up into droplets within a T-junction or some other antechamber, as the form of spherical droplets is energetically favorable when compared to a continuous stream of fluid. Hereby the applied pressure should be neither too small nor too large, but in a range so that the system is in the so-called dripping regime, in which droplets of fluid are produced in equal time intervals [11]. The size of the droplets can be controlled by the flow rates of the two fluids. In the experiments of our collaborating group in Cardiff, the droplets leave the antechamber, then move lined up in an almost one-dimensional ordering, and enter an expansion chamber while several of them are surrounded by some newly generated hull. Within this capsule, the droplets rearrange themselves in a three-dimensional way [11], as shown schematically in Fig. 1.

This paper is organized at follows: in the next two sections, we describe the steps to be taken before we can start developing a chemical compiler. For this purpose, we first need to simulate the arrangement process and second to compare the resulting configurations. In the last section, we give an outlook to the development of the chemical compiler itself.

2 Simulating the Arrangement Process

In order to study this process and its outcome, we develop a computer simulation. After a short overview of existing and widely used simulation techniques, we present our plan for the generation of a new Monte Carlo Movement Simulation Technique.

2.1 Simulation Techniques for Microfluidic Systems

Over the past years, various approaches for simulating droplets moving in fluids have already been developed, from macroscale approaches, in which not single droplets but only droplet densities are considered, to microscale approaches, in which the state variables of the various droplets are changed gradually and individually. Often the methods had been originally developed for other systems but then adopted for the application to microfluidics.

One macroscale approach is the Lattice Boltzmann method [12, 13], with which the time evolution of the density and velocity field of a fluid is simulated on a two-dimensional or three-dimensional lattice. Alternately, collision steps and streaming steps are applied. For the collisions, often the simplified Bhatnagar-Gross-Krook relaxation term [14] is used. While part of the density remains at its current lattice site, other parts are then usually allowed to flow to all sites within a Chebyshev distance of 1, i.e., not only the directly neighboring sites but also the diagonally displaced neighbors are used for flow directions. The huge advantage of this model is that it is very fast and ideally suited for parallel enablement, such that only small amounts of computing time are needed for simulations. The disadvantages are that sometimes lattice artefacts occur and that one has to find out about appropriate rules for flows in various directions. Sometimes two lattices displaced by half a lattice unit in all spatial directions or even more lattices are used instead of one lattice only. A further disadvantage is that this method only considers densities of droplets but not the singular droplets themselves. Thus, this method is not applicable for our investigations, as we need to know about the exact locations and velocities of the various droplets.

On the other hand, the probably most microscopic but also most computer time consuming method is Molecular Dynamics [15, 16]. Hereby, the forces between the various particles are considered. The velocities and the locations of the particles are iteratively and simultaneously updated using specific problem-dependent time integrators, which e.g. preserve the total energy of the system. While this method is suited for considering our problem on an atomic or molecular level, we are unable to use it due to the large system size on the one hand and the lack of computing time on the other hand. Thus, we now turn our attention to two types of simulation techniques with an intermediate requirement for computing time, but also with the possibility to simulate the movements of the various droplets in a way that their exact locations and velocities can be determined exactly.

A wide variety of Monte Carlo simulation techniques has been applied to microscopic simulations of discs in two dimensions and spheres in three dimensions for decades, see e.g. [17]. A subclass of these techniques is called the Direct Simulation Monte Carlo technique [18]. It can be applied to study movements in systems for which the mean free path of a particle is of the same size or larger than its representative physical length scale. The method assumes that free movement phases and collision phases can be decoupled over time periods that are smaller than the mean collision time. There are various ways to model collisions [19], some of them seeming to be rather artificial. A widely used modeling of collisions assumes the particles to be point particles. Then a small box is created around a randomly chosen particle. One of the other particles within this box is then randomly chosen for the simulation of a collision process. The velocities of the two particles are taken into account and the rules for a direct collision of these particles if they were point particles are determined. However, some randomness is added to the direction of the relative velocity vector before updating the velocity vectors of the two particles in order to mimick also a non-direct collision of extended particles in a random way. There are also other Monte Carlo approaches like the Griesbauer method [20]. It introduces springs between the various particles, such that we also consider this method not to be applicable for our problem, as the various droplets move entirely independent of each other at the beginning in the experiments, as can be observed in movies generated by Jin Li [21]. Only at a later stage when they are already surrounded by some hull, the droplets gradually settle down, reducing their individual behaviors, and start to move coherently.

A further widely used approach called Dissipative Particle Dynamics [24,25,26] attempts to relate macroscopic non-Newtonian flow properties of a fluid to its microscopic structure. For the determination of the velocity of the particles, three types of forces are considered which act on a particle. A particle interacts with all other particles within some predefined cut-off distance. There are conservative forces with which the particles interact, then there is a dissipative force, and finally also a random force with zero mean is added. The dissipative and random forces can be chosen in a way that they form a thermostat keeping the mean temperature of the system constant. By choosing the random force between each pair of particles in a way that it acts antisymmetric on both particles as required by Newton’s third law, the local momentum of the particles is conserved. Also this technique is well suited for parallel enablement if using spatial domain decomposition. The diameters of the various domains of course need to be much larger than the cut-off distance. However, artefacts due to the spatial decomposition can occur and the geometry of the experiments to be simulated can become rather complex, such that we consider also this method not to be applicable for our problem.

2.2 A New Monte Carlo Movement Simulation Approach

Fig. 2.
figure 2

From left to right: Two droplets in a hull can either stay standalone or touch each other, forming a bilayer to some smaller or larger extent, and either stay almost spherical or lose their spherical shapes.

Full size image

Summarizing, we intend to create our own simulation technique, with which we want to simulate the experimentally found transition from an originally one-dimensional lineup of droplets into a three-dimensional arrangement. During this rearrangement process, some droplets touching each other will form bilayers [27]. These bilayers can be broken up and reformed, depending on the stability of the bilayers [28]. When bilayers are created, the droplets can lose their spherical shape, as shown in Fig. 2. We will test various ways to simulate the formation, change, and destruction of bilayers and the change of the shape of the cores in a computationally not too expensive way. While the specific spatial setup of an experiment with proposed values for widths and lengths of various parts of the junction can be easily employed also in the Monte Carlo simulation, it is a harder task to find appropriate values for the probabilities for deceleration and acceleration of droplets as well as for bilayer formation and destruction and also for some introduction of random movement. These values depend on various experimental parameters, like pressure and viscosity, and also on the various radius values of the droplets. We intend to adjust the parameters for the Monte Carlo simulation in a way that the resulting configurations reflect the three-dimensional arrangements of droplets as found in experiments.

3 Comparing Resulting Configurations

Fig. 3.
figure 3

Top: Two resulting three-dimensional arrangements of droplets filled with various chemicals. Bottom: Corresponding bilayer networks between the droplets. We made these bilayer networks visible by reducing the size of the droplets and printing a connecting edge between a pair of neighboring droplets if they have formed a bilayer.

Full size image

After this first part of our objective has been achieved, we have a closer look at the resulting arrangements. As the experiments performed by Jin Li have already shown [21], there is not the one and only resulting three-dimensional packing of droplets. Instead, various arrangements are possible, as depicted in Fig. 3. However, when looking closely at the resulting configurations, one finds that they are not entirely random but often rather similar to each other and that maybe even some configurations can be considered as part of a group of configurations having several properties in common. Whether such groups of configurations have an entirely identical backbone [29,30,31] in common or whether they share some properties with some larger probabilities, as found for dense packings of multidisperse systems of hard discs [32] will be seen. We also need to question the influence of the excess of polydispersity of the radius values.

Fig. 4.
figure 4

Left: Overlap matrix with a structure dominated by iterated replica symmetry breaking. Right: corresponding ultrametric tree.

Full size image

Complex systems often exhibit the property of ultrametricity in configuration space [22, 23]. A standard metric d has to obey to the triangle inequality

$$\begin{aligned} d(i,j) \le d(i,k) + d(k,j), \end{aligned}$$
(1)

with d(i, j) denoting the distance from node i to node j, i.e., a direct connection cannot be longer than a detour via a third node k. For an ultrametric, this inequality is replaced by the ultrametricity condition

$$\begin{aligned} d(i,j) \le \textrm{max} \{ d(i,k), d(k,j) \} . \end{aligned}$$
(2)

If permuting the nodes i, j,  and k, one finds that this condition is fulfilled if the nodes are placed on the edges of equilateral triangles or isosceles triangles with short base. A distance between two configurations can be defined using an overlap measure between the configurations. The larger the overlap is, the smaller is the distance. After the application of an appropriate permutation of these configurations, which can e.g. be found with an optimization technique leading to a clustered ordering of configurations [33], the overlap matrix can exhibit a structure as schematically depicted in the left graphics of Fig. 4. In this example, one finds that the overall set of configurations is split in two large groups of configurations. The overlap values of the configurations to other configurations in the same group are larger than those to configurations in the other group. Each group can then be split in four subgroups in this example, which in turn exhibit even larger overlap values within each subgroup. The subgroups are then split again. In statistical physics, one speaks of iterated replica symmetry breaking if such a behavior is observed. Ideally, this replica symmetry breaking property corresponds to ultrametricity, i.e., one can derive also another representation by generating an ultrametric tree. The right graphic in Fig.  4 shows such a tree. At the root, the tree splits into two branches, which in turn split into four subbranches each. These subbranches then split again into three subbranches each. The various configurations then form the leaves of the tree on the right side.

As the property of ultrametricity was also found for a related hard disc packing problem [34], we expect it will also turn up for this problem. As already mentioned, ultrametricity is related to iterated replica symmetry breaking and the possibility to generate ultrametric trees. For their generation, we will use the neighbor-joining method, which is a standard tool to reconstruct phylogenetic trees [35, 36], as well as finding a clustered ordering of configurations [33].

4 Final Steps Towards a Probabilistic Biochemical Compiler

If we have achieved this second part of our objective of understanding and predicting the outcome of an experimental setup, i.e., when the various groups of three-dimensional arrangements of droplets have been generated, we will be able to create a probabilistic chemical compiler in the final stage of this project. We aim at creating plans for e.g. a step-wise generation of some desired macromolecules, which are gradually constructed from smaller units, contained in the various droplets, with the successive chemical reactions enabled via the bilayers formed between neighboring droplets. Thus, the compiler has

  • to determine bilayer networks with which the desired reaction chains leading e.g. to the macromolecules we want to produce can be performed and

  • to design and to govern the experiment leading to such a bilayer network.

Such a compiler has been exemplarily already developed for one specific molecule [37]. In this project, this compiler has to be generalized and also made probabilistic because of the various possible outcomes in the rearrangement process.