Inﬂuence of the geometry on the agglomeration of a polydisperse binary system of spherical particles

Within the context of the European Horizon 2020 project ACDC


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 accumulated (Alberts et al., 2014).Even a new research field of synthetic biology has evolved (Gibson et al., 2017), in which natural objects like the DNA in cells are purposely altered or replaced in order to achieve some desired outcome, like producing a specific drug.
In our approach within the European Horizon 2020 project ACDC -Artificial Cells with Distributed Cores to Decipher Protein Function, we do not consider fully equipped cells but simplified cell-like structures.In the experiments in Cardiff, droplets comprised of some fluid, containing some chemicals, and surrounded by another fluid are considered (Li and Barrow, 2017), partially also being contained within some outer hulls, playing the role membranes have for cells.In Trento, mixtures of emulsion droplet populations are considered in which aggregates were found exclusively if the ssDNA oligonucleotides were of complementary sequence (Hadorn et al., 2012), as exemplarily shown in Fig. 1.Neighboring droplets, i.e., droplets whose midpoints have a distance close to the sum of their radii, can form connections, either by forming bilayers as in the experiments in Cardiff or by getting glued to each other by matching pairs of single-stranded DNA as in the experiments in Trento.Chemicals contained within the droplets can move to neighboring droplets through pores within these bilayers.Thus, a complex bilayer network is created, with the droplets being the nodes of this graph and the existing connections being represented by edges between the corresponding droplets.In such bilayer networks, a controlled successive biochemical reaction scheme can be accomplished to produce the intended macromolecules.The main task of the theory and simulation group in Winterthur is to simulate the various agglomeration processes occurring in the experiments (Schneider et al., 2020b) and to investigate the properties of the resulting bilayer networks (Schneider et al., 2020a) in order to finally create a chemical compiler designing an experiment for producing a desired outcome with a specific agglomeration of droplets (Flumini et al., 2020;Weyland et al., 2020).
Within the scope of this paper, we present computational results for basic simulations of a simplified agglomeration process of a polydisperse binary system of droplets, mimicking experiments performed in Trento.Here we want to focus on the question which influence the shape of the container and the initialization has on the agglomeration process of the particles.In order to focus on these questions and to exclude effects generated by other experimental properties, we simulate the droplets as hard spheres and ignore details of the surface structures of the particles, attractive forces as well as adhesion effects.As the extension of the bilayers is very small and as the droplets keep their spher- ical shape during the experiment, this simplified approach is justified.Systems of hard disks and spheres are widely used in computational physics as simple two-and threedimensional models for granular matter, colloidal systems, fiber-reinforced composites, and molecular crystals (Reiss et al., 1996;Russel et al., 1989;Metcalfe et al., 1995;Zallen, 1983).Usually, monodisperse and bidisperse sytems are considered, i.e., either all disks and spheres exhibit the same radius value or one of two different values.In recent years, the focus of research has shifted to densest packings of items with shapes for which the determination of overlappings between them is slightly more difficult, such as ellipsoids and spherocylinders.It has, e.g., been found that a random packing of ellipsoids with a specific aspect ratio (M&M candies) is denser than a random packing of spheres (Donev et al., 2004).Densest packings of multidisperse hard disks within some environment have been determined (Müller et al., 2009;Schneider et al., 2009).Usually, these packing problems are considered in either two or three dimensions.But also the packing of monodisperse hard spheres in higher dimensions is of interest as this problem can be mapped on finding efficient binary codes for digital communications (Conway and Sloane, 1984;Sloane, 1984).

Simulation details
Figure 2: Random initializations of 2,000 spherical particles with radii randomly chosen in the range 10-50 µm inside a cylinder of height 4mm and radius 1mm: In the "narrow" case on the left, all particles are initialized inside a virtual inner cylinder with radius 0.5mm, whereas the full width of the cylinder is used in the "wide" initialization on the right.
We consider a polydisperse binary system of N = 2000 oil droplets in water, which are modeled as hard spheres as mentioned above.There are two types of particles, Aparticles and B-particles, painted with the colors red and green in Figs.2-4.The binary character of the system is reflected in the way that if an A-particle touches a B-particle then a connection between them is created.Thus, there are A↔B-connections between pairs of A-and B-particles, but no A↔A-connections between pairs of A-particles and no B↔B-connections between pairs of B-particles.Please note that these A↔B-connections only exist virtually due to an A-particle touching a B-particle.But this connection has no further effect and is immediately destroyed if the involved particles do not touch each other anymore.To each particle, the color A or B is randomly assigned with equal probability.
The radii of the particles are randomly chosen from a uniform distribution in the range 10-50 µm.Initially, they are randomly placed in a cylinder of height 4 mm and radius 1 mm without overlaps among particles or between particles and walls.In a second scenario, the initialization is performed in the way that the particles are more closely randomly placed within a virtual inner cylinder of radius 0.5 mm, again without overlap.These two initialization scenarios, which we will refer to later on as the "wide" and the "narrow" scenario, are shown in Fig. 2. Furthermore, we consider two types of scenarios, in which we either use the already mentioned cylinder as container or in which we attach a halfsphere to the bottom of the cylinder in order to study the aggregation process both for flat and spherical bottoms, such that we have all in all four scenarios which we consider.The particles are initialized with zero velocity.
The simulation is divided in time steps of ∆t = 10 −5 s.During each time step, the particles are subjected to various forces: All particles are filled with oil of the density Figure 3: Intermediate configurations for a system with a "narrow" initialization in a cylinder after 1, 2, 3 (left), 5, 10, and 20 (right) seconds.The camera is directed at the bottom of the cylinder in order to observe the agglomeration process closely.= 1.23kg/l, such that they sink in water due to gravity reduced by the buoyant force.Furthermore, the three spatial components of the velocity vectors v i are subject to random velocity changes.Each component v x,y,z i is changed independently by a value equal to its current value times a random number uniformly chosen between −0.05 and 0.05.The particles are also subjected to the Stokes friction force with the radius r i of the particle, the current velocity v i , and the viscosity η = 0.891mPas of water at 25 • C. The concept of added mass is used (Stokes, 1851).This virtual mass is the inertia added to the mass of the particle because an acceleration or deceleration of a body in water must move or deflect some volume of surrounding fluid as it moves through it.It can be shown that the added mass for a spherical particle of radius r i is given by 2 3 πr 3 i fluid , i.e., it is half the mass of the fluid displaced by the particle.
After the new velocities of the particles are determined this way, their positions are updated according to (2) Due to the stochastic nature of the process as imposed by the random velocity changes, only this Euler scheme is suitable for the determination of velocities and positions (Kloeden and Platen, 2013).
After the determination of these new positions, some checks are imposed: First, it is checked whether a particle collides or even overlaps with a boundary of the cylinder or of the cylinder with the halfsphere attached at its new position and whether the overlap will increase if the velocity vector of the particle remains unchanged.If this is the case, the collision normal is determined and the velocity vector of the particle is updated according to the standard collision rules with an elasticity factor of 0.9.If there is an overlap, the position of the particle is updated in order to resolve the overlap.Analogously, then checks for collisions and overlaps between pairs of particles are performed and their velocity vectors and positions are updated accordingly.The overlaps have indeed to be resolved as otherwise they partially remain and can even increase over time, especially in the regime of slow velocities.
In total, 10 7 time steps are performed, which corresponds to a total time of 100 seconds.Please note that this time scale is sufficient as the smallest particles have a radius of at least r smallest = 10µm in our simulations.The time scale to finish the agglomeration process scales with ∼ r −2 smallest .Thus, more time is needed for even smaller particle radii, as in the experiments in Trento.The computing time for one simulation took roughly 12 hours on a standard laptop.Figure 3 shows some intermediate configurations for the scenario with "narrow" initialization and a cylinder with flat bottom as container.Figure 4 exhibits final configurations after 100 seconds for two scenarios with "narrow" initialization, both for a cylinder with flat bottom and for a cylinder with halfsphere attached as container.

Computational results
We are mainly interested in the temporal development of the network formed by the connections between the various droplets.We consider a connection between two neighboring droplets i and j to be formed if (3) Then we define an adjacency matrix η according to η(i, j) = 1 if a connection between i and j exists 0 otherwise.
This adjacency matrix contains all the information about the network formed by the droplets.(Please note that η is symmetric in our case, but we can extend the approach to directed connections and asymmetric η.) In order to compare unary systems (all particles have the same color) with binary systems (with the Aand Bparticles as explained above), we consider two different adjacency matrices, namely Figure 5: Time evolution of the average degree of the particles in the four simulation scenarios considered both for a unary system wih only one sort of particles ("1 color") and a binary system with A↔B connections only ("2 colors") • a matrix η 1 for the unary system, in which η 1 (i, j) = 1 if a connection between two nodes can be created if condition (3) is fulfilled, regardless the colors of i and j, and • a matrix η 2 for the binary system, in which η 2 (i, j) = 1 only if both condition (3) is fulfilled and droplets i and j have different colors.
In the graphics, we will refer to the results obtained for η 1 with the term "1 color" and for η 2 with the term "2 colors".There are various ways to investigate networks.One can focus on the local aspect and have a look at the neighborhoods of the various particles.Another approach is to look at the intermediate scenario, i.e., at groups of particles within a network, e.g., pairs or triples or other very small groups of particles.On the other hand, one can take the global perspective and consider the network as a whole and deal e.g. with questions as whether the network is percolating (Stauffer and Aharony, 1994;Stauffer, 1986;Naftaly et al., 1991).Last but not least, there is also the possibility to use a local-global scenario and to investigate the importance of specific nodes for the overall network or the role some cen-tral nodes play in the network (see e.g.Schneider and Kirkpatrick (2005)).In this paper, we exemplarily present results for the local and the global perspective.
We want to start with the time evolution of the average degree of nodes.The degree d i of a node i is defined as the number of nodes it is attached to via edges in the graph.Thus, the degree of a particle i in network η is given by and we can compute the mean value This mean value can also be considered as part of the global perspective on a network, as the overall number e of edges in the graph is given by e = i<j η(i, j), (7)  Figure 6: Time evolution of the numbers of clusters in the four simulation scenarios considered both for a unary system wih only one sort of particles ("1 color") and a binary system with A↔B connections only ("2 colors") such that we get the relation between the global parameter e and the local parameters d i .
Nevertheless, this part of the investigation of a network is usually seen as from a local perspective.
Figure 5 shows the mean value d of the degrees of the various nodes, both for the unary and the binary system, for the four simulation scenarios we consider.For all four scenarios, we find a sigmoidal increase of d over time.All four graphics can basically be divided into a short-time regime up to roughly one second, in which hardly any or rather few connections are created, and a long-time regime between one and 100 seconds, in which the average degree of the nodes increases sigmoidally till a final value.This final value is significantly largest for the scenario with a "narrow" initialization inside the cylinder to which the halfsphere is attached, both for the unary and for the binary system, which is not surprising, if we compare the final configurations depicted in Fig. 4. The half-spherical bottom allows a closer packing of the droplets in z-direction.However, this effect is more pronounced in the case of a "narrow" initialization.Thus, we already find here that both the temporal behavior of the increase of the average degree of nodes and the final values strongly depend both on the shape of the container and on the initialization process.
In the next step, we have a look at the number of clusters.Studying clusters in networks belongs to the class of global network analysis methods.An easy way to detect clusters is to assign a flag f i to each node i and to derive a neighborhood list for each node at the beginning.A node j belongs to the neighborhood of node i if η(i, j) = 1.Initially, all flags are initialized with f i = 0 and the number M of clusters is set to zero.Then one performs an iteration over the nodes.If the condition f i = 0 holds for node i, then this node is not yet part of a cluster.Then the number M of clusters is incremented by 1 and one sets f i = M in order to mark that node i belongs to the new cluster with cluster number M .Then one goes through all nodes j in the neighborhood list of node i and sets f j = M as well.This inner loop is recursively repeated, among neighbors, neighbors of neigh-bors, and so on, until no further neighbor with f j = 0 can be found.Finally, all flags f i are non-zero and M equals the number of clusters in the system.Please note that according to this definition, a cluster can consist of one node only.
Figure 6 shows the temporal decrease of the number of clusters for the various scenarios from a common initial value of M = N .Again we can distinguish a short-time phase till roughly one second and a long-time regime for longer time durations in which the number of clusters decreasses strongly in a sigmoidal way.Please note that the curves for the unary system and the binary system are rather close to each other for times t ≤ 10s in case of the wide initialization and for times 2s ≤ t ≤ 20s in case of the narrow initialization.Thus, a huge gap between the curves for the average degree as shown in Fig. 5 does not lead to a huge gap here for the number of clusters.Furthermore note that it is very likely that only one large cluster containing all nodes remains at the end for the unary system, especially in case of the cylinder with a halfsphere attached to it.However, some small isolated droplets can lie at the bottom of the cylinder, without being connected to the large cluster, if only a cylinder is used as container, such that one does not necessarily get M = 1 for the unary system at the end.For the binary system, one finds not only the temporal behaviors but also the final values for the numbers of clusters depend both on the shape of the container and on the initialization process.The largest number of clusters remains for the "narrow" initialization with a cylindrical container, whereas the smallest number of clusters is achieved for the "narrow" initialization with a container consisting of a cylinder and a halfsphere.
Finally, we wonder what the configurations look like if up to 200 different clusters remain till the end of the simulation.We are especially interested in the question whether there is a dominating cluster also in the case of the binary system.Thus, we protocol the sizes of the clusters during our cluster detection and have a look at the size of the largest cluster.Figure 7 shows the time evolution of the size of the largest cluster.Again we see that the time evolution can be separated into a short-time regime up to roughly 1 second and a long-time sigmoidal increase and that the curves for the unary systems and their binary counterparts are again rather close to each other.However, the main result of this investigation is that we finally achieve a dominating cluster comprised of more than 90% of the particles in all scenarios for the binary system.

Summary and outlook
In this paper, we presented results of simulaions for the agglomeration of droplets.As we are interested in the effects of the container shape and of the initialization process, we study a very simplified system, in which the droplets are represented as hard spheres, subjected to gravity reduced by buoyancy as well as the Stokes friction, the added mass effect, random velocity changes, and almost-elastic impacts.
Connections between these particles are virtually formed if they touch or overlap.The particles gradually agglomerate at the bottom of the container.The analysis of this agglomeration process from a local and a global point of view shows that the temporal results strongly depend on the shape of the container and on the initialization process, but generally, we find a dominating cluster for a binary system.
We intend to extend our investigations on ternary systems, in which only A↔B-connections can be formed between pairs of particles, while the C-particles are unable to make connections.Hereby we want to study the breakdown of the size of the largest cluster with increasing density of Cparticles.Furthermore, we want to add gluing forces between particles to find out how they change the results.Figure 7: Time evolution of the sizes of the largest cluster in the four simulation scenarios considered both for a unary system wih only one sort of particles ("1 color") and a binary system with A↔B connections only ("2 colors")

Figure 1 :
Figure 1: Exemplary agglomeration of droplets at a bottom of the well in an experiment (top) and clusters of droplets in closeup view via confocal microscopy (bottom)