Percolation Breakdown in Binary and Ternary Monodisperse and Polydisperse Systems of Spherical Particles

Abstract

We perform computer simulations of an agglomeration process for monodisperse and polydisperse systems of spherical particles in a cylindrical container, using a simplified stochastic-hydrodynamic model. We consider a ternary system with three particle types A, B, and C, in which only connections of the type $A-B$ can be forged, while any other connections with particles of the same type or with C-particles are forbidden, and for comparison a binary system with two particle types A and C, in which only connections of the type $A-A$ can be formed. We study the breakdown of the percolation in the agglomeration at the bottom of the cylinder with an increasing fraction of C-particles.

This work has been partially financially supported by the European Horizon 2020 project ACDC – Artificial Cells with Distributed Cores to Decipher Protein Function under project number 824060.

1 Introduction

Fig. 1.

Left: Snapshot of an agglomeration of droplets recorded from an experiment. Right: Agglomeration of a polydisperse system of 2,000 spherical particles with types A and B, depicted as red and green, at the bottom of a cylindrical container, obtained in a computer simulation. (Color figure online)

We intend to develop a probabilistic compiler [3, 22] to aid the three-dimensional agglomeration of droplets filled with various chemicals (see Fig. 1) in a specific way in order to e.g. allow the creation of desired macromolecules via a successive reaction scheme [12, 13, 18, 19]. Neighboring droplets can form connections, either by forming bilayers [7] or by getting glued to each other by matching pairs of single-stranded DNA [4], as sketched in Fig. 2. Chemicals contained within the droplets can move to neighboring droplets either directly, as hydrophobic compounds can be exchanged between adjacent oil droplets at the contact face, or, if the oil droplets are contained in a hull comprised of amphiphilic molecules like phospholipids, through pores within these bilayers. An example for such a pore is shown in Fig. 3. Thus, a complex bilayer network is created [16], with the droplets being the nodes of this graph and the existing connections being the edges between the corresponding droplets. In such bilayer networks, a controlled successive reaction scheme can be effectuated to produce the intended macromolecules. As already demonstrated for a toy example, a gradual reaction network with three educts, two reaction steps, a desired product, and an undesired side product, can achieve a higher yield and a smaller amount of undesired side products in an agglomeration of droplets with restricted molecule transfer than in a scenario in which all educts would be put in one well-stirred pot only [17].

Fig. 2.

Sketch of a pair of oil-filled droplets in water, to which complementary strands of ssDNA oligonucleotides are attached: The surfaces of the droplets are composed by single-tail surfactant molecules like lipids with a hydrophilic head on the outside and a hydrophobic tail on the inside, thus forming a boundary for the oil-in-water droplet. By adding some single-strand DNA to the surface of a droplet, it can be ensured that only desired connections to specific other droplets with just the complementary single-strand DNA can be formed. Please note that the connection of the droplets in this picture is overenlarged in relation to the size of the droplets. In reality, the droplets have a radius of 1–50 $\upmu \text {m}$, whereas a base pair of a nucleic acid is roughly $0.34\,\textrm{nm}$ in length [1], such that the sticks of connecting DNA strands are roughly $5\,\textrm{nm}$ long.

Fig. 3.

Alpha hemolysin (aHL) pore: The left picture reveals how the aHL, which is comprised of seven macromolecules, sticks its trunk through the bilayer between two droplets. The right picture presents in detail a cut through the pore formed, revealing the channel through which molecules can move between the adjacent droplets.

For some applications, it is necessary to thin out the network, i.e., to reduce the number of edges connecting nodes in the network. This leads as we will show later to smaller numbers of nodes a node is attached to on average, to smaller cluster sizes and in turn to larger numbers of clusters which are isolated of each other, just as we need them for these applications, in which we either need to better govern the gradual chemical reaction process or in which we need to strongly reduce the maximum number of steps within such a gradual process. In order to achieve this thinning-out, particles which do not connect to any other particle can be added to the system. Within the scope of this paper, we study the effects of these auxiliary particles on the properties of the overall network for two basic scenarios:

• “Ternary” scenario:

  In the ternary scenario, systems with three particle types A, B, and C are considered. Only connections between neighboring A- and B-particles can be forged. Besides these $A-B$-connections, no other connections to particles of the same type or to C-particles are possible. Without the C-particles, this system would form a so-called bipartite network. In our simulations, we only consider the special case that the fraction $f_A$ of A-particles equals the fraction $f_B$ of B-particles. We will study the changes of the properties of the network for an increasing fraction $f_C$ of C-particles, with $f_C=1-f_A-f_B=1-2f_A$.

• “Binary” scenario:

  For comparison, we also study a binary scenario with two particle types A and C, in which only pairs of neighboring A-particles can form connections. Besides these connections of the type $A-A$, there are no other connections, such that also here the C-particles serve as auxiliary particles for thinning out the network. Again we want to study the effects of an increasing fraction $f_C$ of C-particles on the properties of the network. This scenario can also be considered as a site percolation problem [20], in which the locations of the A-particles represent the occupied sites and the locations of the C-particles represent the empty sites. The probability p for an occupied site is simply given by $p=f_A=1-f_C$.

In our computer simulations, we study both polydisperse systems, in which the radii of the particles differ from each other, and monodisperse systems, in which all particles share the same radius value, in order to mimick experiments of various kinds: The production of droplets using a microfluidic approach, in which an inner stream of fluid within an outer stream of another fluid is broken up in droplets in e.g. a t-junction under specific pressure conditions [7], leads to a rather monodisperse system of droplets. Contrarily, in other experiments, we are repeatedly rubbing a phial filled with water and one drop consisting of oil molecules and amphiphilic molecules over a rough surface, thus sending excitations into the system, which lead to a breakup of the large drop into many small droplets of varying sizes, resulting in a polydisperse system [4].

Within the scope of this paper, we present computational results for simulations based on a simplified stochastic-hydrodynamic model of an agglomeration process of a system of droplets, mimicking experiments. Here we want to focus on the influence of the fraction $f_C$ of auxiliary C-droplets on some specific properties of the networks created, which are of crucial importance for the gradual reaction scheme intended. We are especially interested in the question whether there is a percolation transition at some critical value of $f_C$: For an infinitely large system, one expects to get a sharp transition between two phases, with an infinitely large cluster for $f_C$ below the critical value and no such infinitely large cluster for $f_C$ above that value. For finite systems, one gets a smooth transition between a regime with a large cluster dominating the system for small $f_C$ and a regime with no such large cluster but many small clusters for large $f_C$. In order to focus on these questions and to exclude effects from other experimental properties, we simulate the droplets as hard spheres and ignore details of the surface structure of the particles, attractive forces as well as adhesion effects. As the extension of the bilayers is very small and as due to their small radii [2], the droplets keep their spherical shape during the experiments, as shown in Fig. 1, such that this simplified approach is justified.

This paper is organized as follows: In the next section, we sketch briefly how we simulate the agglomeration of droplets in a container. Then we give a short introduction to network analysis and percolation theory, focusing on those network properties for which we will present computational results. As we are mainly interested in the description of the percolation transition with an increasing fraction $f_C$ of auxiliary C-droplets, we present results depending on $f_C$ for the decreasing maximum and average number of nodes a node is attached to, for which we find power laws depending on $1-f_C$, for the increasing number of clusters, for which we get a linear behavior, and for the decreasing size of the largest cluster, which clearly exhibits a percolation transition. Finally, we provide a summary and give an outlook.

2 Simulation Details

At the beginning of the simulations, we place N spherical particles at randomly selected positions in a cylindrical container with radius $1\,\mathrm{{mm}}$ and height $4\,\mathrm{{mm}}$ in a way that they do not overlap with each other and that they do not overlap with the walls of the cylinder. For the polydisperse system, we randomly choose the particle radii $r_i$ uniformly from the interval [10–50] $\upmu \mathrm{{m}}$, whereas we set all radii $r_i \equiv 30\,\upmu \mathrm{{m}}$ for the monodisperse system.

After this initialization, we perform the main simulation which is comprised of $10^7$ time steps of a duration of $\delta t = 10^{-5}\,\textrm{s}$. In each time step, the particles are subjected to various forces:

• They sink in water due to gravity $\vec {F}_G$ reduced by the buoyant force $\vec {F}_b$:

  $$\begin{aligned} \vec {F}_G(i)-\vec {F}_b(i)=\frac{\displaystyle 4\pi }{\displaystyle 3} r_i^3 \left( \varrho _\textrm{oil} -\varrho _\textrm{water} \right) g \end{aligned}$$

  (1)

  For the oil density, we use the value $\varrho _\textrm{oil}=1.23\,\mathrm{kg/l}$, which is just the density of the oil used in some experiments.

• Secondly, the spatial components $v_{x,y,z}(i)$ of the velocity vectors $\vec {v}(i)$ are subjected to random velocity changes: They are randomly altered by up to ${\pm }5\%$ of their absolute values in order to take at least in this small random way into account that the containers are moved by the experimentalists during the agglomeration process.

• The particles are also subjected to the Stokes friction force $\vec {F}_S$:

  $$\begin{aligned} \vec {F}_S(i)=-6 \pi \eta r_i \vec {v}(i) \end{aligned}$$

  (2)

  The viscosity of water at 25 $^\circ$C is $\eta =0.891\,\textrm{mPas}$.

• As in classic hydrodynamics, the concept of added mass [21] is used. When applying Newton’s second law, we have to consider an effective mass of the particle, i.e., $\vec {F}(i) = m_\textrm{eff}(i) \vec {a}(i)$. This effective mass is composed of the mass m(i) of particle i and of the added mass $m_\textrm{added}(i)$. This added mass is caused by the inertia of the surrounding fluid, which needs to be deflected or attracted if the particle itself is accelerated or decelerated in the water, and can be determined to being half of the mass of the water displaced by oil particle i.

When working with such a set of second order differential equations governing the laws of motion for the particles, the question arises as to which integrator to use. Due to the stochastic nature of random velocity changes, only an Euler scheme with very small time intervals is suitable for the determination of new velocities and positions [6]. In the case of collisions between pairs of particles or between particles and walls, a mostly elastic collision dynamics with 90% elasticity and 10% plasticity is imposed. Overlaps occurring at the end of each time step are resolved as in [9, 14].

3 Network Analysis and Percolation Theory

For network analysis, we first of all have to define a network related to the problem we intend to study. As mentioned above, we are interested in generating gradual chemical reaction schemes performed in networks of droplets, with neighboring droplets being able to exchange molecules if pores within their bilayers exist or, more theoretically speaking, if a connection between the particle types of the two adjacent droplets is allowed. Then we can define an edge matrix $\eta$ with

$$\begin{aligned} \eta (i,j)= {\left\{ \begin{array}{ll} 1 &{} \text { if droplets } i \text { and } j \text { are neighbors of each other}\\ &{} \text { and a connection between them exists}\\ 0 &{} \text { otherwise} . \end{array}\right. } \end{aligned}$$

(3)

Two droplets i and j with their midpoints $(x_i,y_i,z_i)$ and $(x_j,y_j,z_j)$ and their radii $r_i$ and $r_j$ are neighboring each other if the condition

$$\begin{aligned} \sqrt{ (x_i-x_j)^2+ (y_i-y_j)^2+ (z_i-z_j)^2} \le r_i+r_j+0.1\,\upmu \textrm{m} \end{aligned}$$

(4)

is fulfilled, i.e., if the distance between their midpoints is smaller or equal to the sum of their radii plus some small offset which we need to introduce because of finite numerical precision. One usually sets $\eta (i,i) \equiv 0$ for all nodes i. Such a matrix $\eta$ contains all the information about the network. For this paper, we study both a binary scenario and a ternary scenario, for which we can generate two different edge matrices, considering the different conditions for the existence of a connection.

When analyzing a network, one mostly takes either an atomistic view, looking at the various nodes and determining their network related properties, or a global view, determining clusters of nodes. Clusters are defined as maximum subgroups of nodes in which each node within this cluster can be reached from any other node in the cluster by gradually traversing edges, thus walking along a path from this node perhaps via other nodes in the subgroup to the destined node. More seldomly, networks are considered at an intermediate level, e.g., for the detection of the maximum clique [8], or one asks for the importance of specific nodes for the overall network in a local-global view, see e.g. [11].

When looking at a network from a global point of view, one of the most important questions arising is whether the network is percolating. For an infinitely large network, this means that one has to ask whether there is an infinitely large cluster in the network [20]. In the finite networks resulting from computer simulations, one thus asks whether a dominating cluster exists in the network.

Mostly, the so-called site percolation is studied in which sites on a regular lattice are either occupied or empty and in which each site is connected to all neighboring sites. Theoretically, one finds for infinitely large systems that there is a critical probability $p_\textrm{crit}$ of occupied sites above which an infinitely large cluster exists in the system and below which there is no percolation anymore. Alternatively, also the so-called bond percolation is considered in which all lattice sites are occupied but only a fraction p of the edges exist. Also here one finds such a critical probability $p_\textrm{crit}$ dividing two such regimes. For some scenarios, this critical probability can be calculated exactly, but mostly, one has to make use of computer simulations with increasing system size, to determine the various clusters in the system and the size of the largest cluster, and finally to carefully determine $p_\textrm{crit}$ numerically [20].

4 Computational Results

The results presented in Figs. 4, 5, 6, 7 and 8 are averaged over the properties of the final configurations of 100 independently performed simulation runs.

The first observable we have a look at is the number e of edges, which can be derived from the edge matrix $\eta$ with

$$\begin{aligned} e=\sum _{i<j} \eta (i,j) . \end{aligned}$$

(5)

Figure 4 displays the results for e for the binary and the ternary scenario in simulations of monodisperse systems of 2,000 droplets and polydisperse systems of 2,000 droplets. We generally find that there are more edges in the binary scenario and that the number of edges decreases with an increasing fraction $f_C$ of C-particles, which do not connect with each other and with other particles. On average, there are slightly more edges in the monodisperse systems than in the polydisperse systems, both for the binary and for the ternary scenarios.

Fig. 4.

Decrease of the number e of edges with an increasing fraction $f_C$ of C-particles for the binary and ternary systems as described in the text: Results are presented for monodisperse systems consisting of 2,000 particles (left) and polydisperse systems consisting of 2,000 particles (right).

When taking a local perspective, one of the most important observables for a specific node i is its degree d(i), which can be calculated as

$$\begin{aligned} d(i)=\sum _{j=1}^N \eta (i,j) . \end{aligned}$$

(6)

Thus, the degree d(i) counts to how many other droplets the droplet i is connected. We are mainly interested in the maximum degree

$$\begin{aligned} d_{\textrm{max}} = \textrm{max}_i \, d(i) \end{aligned}$$

(7)

of all nodes. Of course, also $d_{\textrm{max}}$ decreases with increasing $f_C$, but in a specific way, such that we plot $d_{\textrm{max}}$ vs. $1-f_C$ in a double-logarithmic way in Fig. 5. The graphics reveal the existence of a power law for $d_{\textrm{max}}$. For the monodisperse systems comprised of 2,000 particles, we find a power law of the type

$$\begin{aligned} d_{\textrm{max}} = a (1-f_C)^{1/2}, \end{aligned}$$

(8)

both for the binary and the ternary scenarios, whereas we get a power law of the type

$$\begin{aligned} d_{\textrm{max}} = a (1-f_C)^{2/3} \end{aligned}$$

(9)

for the polydisperse systems comprised of 2,000 particles, both for the binary and for the ternary scenarios. The values for the various prefactors a are provided in Table 1.

Fig. 5.

Increase of the maximum degree $d_{\textrm{max}}$ vs. the remaining fraction $1-f_C$ of particles not being C-particles for the binary and ternary systems as described in the text: Results are presented for a monodisperse system consisting of 2,000 particles (left) and a polydisperse system consisting of 2,000 particles (right).

Table 1. Prefactors a found for the power laws as described in the text.

In the next step, we have a look at the mean value $\langle d \rangle$ of the degrees, which can be calculated as

$$\begin{aligned} \langle d \rangle = \frac{\displaystyle 1}{\displaystyle N} \sum _{i=1}^N d(i) . \end{aligned}$$

(10)

$\langle d \rangle$ is related to the overall number e of edges via

$$\begin{aligned} N \times \langle d \rangle = 2 e. \end{aligned}$$

(11)

Also for $\langle d \rangle$, we find a complex power law behavior depending on $1-f_C$, such that we plot $\langle d \rangle$ vs. $1-f_C$ in a double-logarithmic way in Fig. 6. We get both for the binary and for the ternary scenario, both for the monodisperse and for the polydisperse systems the power law

$$\begin{aligned} \langle d \rangle = a(1-f_C)^2 . \end{aligned}$$

(12)

The prefactors can again be found in Table 1.

Please note that both the maximum degrees and thus also the mean degrees are restricted in size. For the monodisperse system, $d_\textrm{max}$ cannot exceed the value of the so-called kissing number k [15]. The kissing number problem is stated as follows: How many spheres of equal size can be placed around a sphere in their midst touching it without any overlaps? This kissing number equals 12 in three dimensions, as already stated by Newton and proved in the 1950s. For the polydisperse system, there is a related restriction: Here $d_\textrm{max}$ cannot exceed a value of k, which depends on the ratio between the radii of the smallest and largest spheres. As the radii of the spherical particles are randomly chosen from the interval [10–50] $\upmu \mathrm{{m}}$, this ratio could be up to 1 : 5, for which we obtained a bidisperse kissing number of 120 [15].

Fig. 6.

Increase of the average degree $\langle d\rangle$ vs. the remaining fraction $1-f_C$ of particles not being C-particles for the binary and ternary systems as described in the text: Results are presented for monodisperse systems consisting of 2,000 particles (left) and polydisperse systems consisting of 2,000 particles (right).

Furthermore, we would like to compare our results for $\langle d \rangle$ for monodisperse systems with results obtained for other configurations of spheres. First of all, let us consider the densest packing of spheres of the same size: The densest packings can be achieved both in a face centered cubic (fcc) lattice and a hexagonal close packing (hcp). For these densest packings, one gets $\langle d \rangle = d_\textrm{max} =12$, in the case of infinitely extended lattices or lattices with periodic boundary conditions. Contrarily, if studying entirely random packings of spheres of the same size, one obtains $\langle d \rangle =6$ [10], a value which is almost in agreement with our values for the binary scenario at $f_C=0$ for the monodisperse system. The deviation is due to the finite extension of our agglomerations.

Now we turn to the global view in network analysis and have a look at the number n of clusters in the system, which is plotted vs. the fraction $f_C$ in Fig. 7. For $f_C=0$, we trivially have only a very small number of clusters in the binary scenario, whereas there is already a significant number of clusters of roughly $n \approx 0.05N - 0.07 N$ in the ternary scenario. For small $f_C$, n increases linearly with $f_C$, until it approaches sigmoidally the value $n=N$ in the limit $f_C \rightarrow 1$. For the ternary systems, there seems to be a little bending in the curves at $f_C \approx 0.5$.

Finally, we end up at the most important point of our investigation. We consider the size $s_\textrm{max}$ of the largest cluster in the system, which is plotted vs. $f_C$ in Fig. 8. Generally, we get a linear decrease of $s_\textrm{max}$ with increasing $f_C$ for small $f_C$, before a transition takes place, in which the percolation breaks down: For the binary scenario, we find a critical value of $f_C$ of roughly $0.55\pm 0.05$ both for the monodisperse and the polydisperse system. For the ternary scenario, we get $0.3\pm 0.05$ both for the monodisperse and for the polydisperse system.

Here we again would like to compare these results with other results obtained for spheres. As already mentioned, the densest packings of spheres of the same size can be achieved in a fcc and a hcp lattice. The critical probability for an infinitely large fcc lattice with periodic boundary conditions has been determined to be $p_{\textrm{crit}}=0.198$ [20]. But also hcp lattices on a slab with open boundary conditions and infinite extensions in two dimensions have been studied. The threshold depends on the thickness h (i.e., the number of layers on which the midpoints of the spherical particles are located) of the slab, one gets $p_{\textrm{crit}}=0.2828$ for $h=2$ and $p_{\textrm{crit}}=0.2086$ for $h=16$ in the limit of infinite extension in the other two dimensions [5]. For randomly packed spheres, a threshold of $p_{\textrm{crit}}=0.31$ [10] was obtained. Thus, we see that the critical values strongly depend on the systems under study. In our case, the spherical particles are neither located on the sites of a regular grid nor placed entirely randomly.

Fig. 7.

Increase of the number n of clusters with an increasing fraction $f_C$ of C-particles for the binary and ternary systems as described in the text: Results are presented for monodisperse systems consisting of 2,000 particles (left) and polydisperse systems consisting of 2,000 particles (right).

5 Summary and Outlook

In this paper, we presented results of simulations for the agglomeration of polydisperse and monodisperse systems of droplets. We were mainly interested in the effects the addition of auxiliary particles, which do not connect to any other particles, has on the networks and their properties. We found a power law behavior for the maximum degrees and mean degrees of the particles depending on the fraction of the auxiliary particles in the system. Furthermore, we detected a percolation breakdown if this fraction exceeds some critical value.

We will continue this study also with other connection scenarios, in which e.g. A-particles can connect to other A-particles and to B- and C-particles, while B-particles and C-particles cannot form connections. This scenario can be easily realized in experiments by only placing the constituents for pore macromolecules exclusively in the A-particles. Furthermore, we will extend our study to further system sizes in order to get better estimates for the critical values and also to find out in which way the prefactors a in the power laws found for the maximum and the mean degree depend on the system size N.

Fig. 8.

Decrease of the size $s_{\textrm{max}}$ of the largest cluster with an increasing fraction $f_C$ of C-particles for the binary and ternary systems as described in the text: Results are presented for monodisperse systems consisting of 2,000 particles (left) and polydisperse systems consisting of 2,000 particles (right).