Team:Vilnius-Lithuania-OG/Merging Model

Collaborations

Mathematical description of droplet coalescence

Explore MODEL

Introducton

CAT-Seq system design incorporates a droplet merging step as a way to enrich the droplets containing the Product Nucleotides by combining them with droplets that contain DNA amplification mix and Reference Nucleotides in high-throughput.

Yet, as we began to grasp the droplet merging techniques, effective merging appeared to be an extremely difficult task. The merging was not efficient – only 1 to 5 droplets merged for every 20 of droplets that contained the Product Nucleotides. Each droplet that has not merged is a huge loss of information. As our goal is to create a biomolecule activity recording system that is not only fast, but also accurate, increasing the effectiveness of the droplet merging step was of the highest priority.

Videos displaying faulty and inefficient droplet merging due to poorly chosen parameters

The merging difficulty stemmed from the fact droplet merging techniques consist of large number of different parameters that depend on each other – droplet volumes, water and oil phase flow speeds, electrical field strength that helps the droplets to merge, microfluidic chip dimensions and other.

For this reason we have decided to develop an in-depth mathematical model to determine the appropriate merging system parameters that would provide effective droplet merging in CAT-Seq.

Here we present the full mathematical description of the active fusion of droplets in a microfluidic chip. The droplets from two populations are fused together one-by-one using specially designed droplet coalescence chip and an electric field. The full scheme of the chip is shown in Figure 1:

The electric field is created by plugging the voltage to the metal electrodes (Fig.1 A). One of the droplet population is being produced on a chip (Fig.1 B) and at the same time, the other is being reinjected in a form of emulsion (Fig.1 C). The droplets in emulsion are being “spaced” by using additional oil flow. Tuning the flow rates of the two droplet populations allow us to have the produced-reinjected droplet ratio 1:1. When the two droplets reach the electrodes in a chip (Fig.1 D), an electric field forces them to coalesce into a single one. The droplet coalescence is a result of the drain of spacing oil between them. This comes from the combined effects of viscous forces and the oil draining created by an electric field.

Our model presented here relies on a balance of energy. The energy terms are those for the kinetic energy and the energy dissipated as a result of viscous forces. In the next sections we explicitly present each of the terms, which are making an impact to the energy balance.

1. Determining the drainage rate

In this section we derive the approximate relation for drainage of oil separating two droplets before their coalescence. We use the relation for pressure difference \(\varDelta p\) the one as used by Barman et al.

(1)

Here \(\gamma_{WO}\) is the interfacial tension of oil-water interface. As we will find later, it is changed by adsorption of surfactant molecules. Notice that all voltages applied for the energy terms rely on voltage created at the electrode tip \(U_{tip}\) and not purely on applied voltage \(U_{app}\).

The action of electric field also changes the interfacial tension through a change in contact angle. Authors postulate the differential equation for oil-film drainage rate \(dT/dt\):

(2)

Integrating the differential equation (2)

We integrate this differential equation by using the mathematical methods of partial integration and guessing of coefficients.

The equation (2) is rewritten in the following form:

(A)

We mark the multiplicative factors of \(T\) as \(a, b\). Then split the right side into “possible” terms with unknown multiplicative factors:

(B)

We solve the system pf algebraic equations to find the coefficients A, B, C. These are found to be \(A=1/ab\), \(B=-1/b\), \(C=-1/a\). The explicit solution for the integral is the following:

(C)

Looking at the coefficient a, b, we find it to have dimensions \([a]=1m^3\), \([b]=1m^3\). This allows us to use only the ln(x)-part of the integral as the left side of equation (A) has the same dimension. Then we do an approximation:

(D)

We find this to have the same dimension as a left part of integration (A):

(E)

The oil used for droplet production has a surfactant dissolved in order to reduce the interfacial tension. In the same way the oil used for spacing the re-injected droplets must have the same fraction of surfactant to make sure their surface keeps saturated with surfactant molecules. We postulate the lowest thickness of the drainage film can be equal to the size of the surfactant monomer [2]\(a \approx 2 nm\). Placing back the thickness of the drainage-film \(x=T\), and rewriting the equation for \(T\), we get a function \(T(t)\):

(F)

Here \(T_0\) is the initial thickness of drainage film at time \(t=0\).

The initial thickness T0 is found from a continuity condition stating the total flow rate in a system is constant. The flow rates for two populations of droplets are expressed as following:

(G) [Re-injected emulsion of droplets V1 with spacer oil]

(H) [On-chip produced emulsion of droplets V2 with oil]

(I) [Emulsion of droplets after coalescence]

Here \(Q_{sp}\) represents the spacer oil flow rates for two droplet populations. In our experiments of droplet coalescence, the flow rates of liquids are tuned in a way to have droplets of two populations flowing into a fusion chamber at 1:1 ratio. We might assume the number of droplets per time unit is equal to one. This means the time in equations (E)-(F) has to be changed into coalescence time tc. The time for two droplet injection is the same unit time \(\varDelta t\), while the time unit for droplet fusion is \(\varDelta t + t\). In our system with idealized-conditions, the volume of spacer oil comes between two droplets through a channel of dimensions H×W1. This unit is then injected into a wider chamber H×Wch (See Fig.2). Using these arguments and the condition of continuity, we write the equilibrium condition:

(J)

The volume of oil between two droplets can be expressed as \(V_{sp} = (Q_1 + Q_2) \varDelta t = l_{init} H W_1 = l H W_{ch}\). Solving the eq.(J) for Δt and inserting it into the volume relation, we get the expression for l, which represents the initial thickness of oil in the coalescence chamber (with dimensions HWch):

(K)

Replacing T0 in eq. (F) with eq. (K), the full expression for drainage-film thickness is found:

(L)

The approximate solution of equation (2) is the following:

(3)

The drainage rate of the oil is expressed according to Barman et al. (See ESI of [1]):

(4)

When the expressions (1) and (3) are inserted, the equation (4) gives a function of drainage flow rate.

2. Interfacial tension and electrowetting-effect

The pressure difference in equation (1) and the expression for Qdrain accounts for the effects of surfactant adsorption and capillary pressure. However in our experiments we have the two populations of droplets with different interfacial tension γ. The difference of interfacial tension arises from different degrees of interface saturation of surfactant between re-injected and produced droplets. The re-injected droplets are incubated beforehand and already have saturated amount of surfactants on their surface. In comparison, the produced droplets only reach the saturation at a later point. All the droplets in our experiments are being produced by using nonionic surfactants.

The first population of droplets is being re-injected into the electrocoalescence chip after some incubation period. Liquid interface of these droplets gets saturated with surfactant monomers from oil solution. It is well known the interfacial tension of droplets reach the lowest (equilibrium) value γeq[3] after incubation. We use the value γ1=γeq for interfacial tension of re-injected droplets. This population of droplets is produced in a microfluidic chip of lower channel height, therefore the volume of droplets is smaller compared to that of the second population droplets. When the first population is re-injected into the coalescence chip with higher channel height compared to the channel height they were initially produced in, the droplets acquire spherical shape. The second population of droplets are produced in the coalescence chip. These droplets acquire a flat circle shape in the channels.

Scheme with a channel, wherein on droplet is flat circle shape (second population) and another is spherical shape (first population)

The interfacial tension of droplets produced with nonionic type of surfactant follows the nonlinear adsorption [4]:

(5)

Microfluidic droplet coalescence is based on application of alternating electric voltage. The value of interfacial tension created by surfactant adsorption is further modified by electrical field effect on a contact angle. This influences the interfacial tension γ for both produced and re-injected droplets.

The interfacial tension of the three-phase interface follows the Young-Laplace equation:

(6)

We also know the electric field impacts the contact angle θ:

(7)

The initial contact angle θ0 (before applying the electric field) is expressed from Young- Laplace equation (6) and inserted into the relation for electrowetting:

(8)

The applied electric field does not impact the interfacial tension between solid (PDMS) and oil γSO or solid and water γSw. We introduce the effective interfacial tension γ(eff) which will account the effect of electrowetting \(\gamma_{WO}^{(eff)}=(\gamma_{SO} - \gamma_{SW})/\cos \theta\). This relation is inserted into equation (8) to give a result:

(9)

In the last equation we find a dimensionless factor called the electrocapillary number ECN, showing the ratio of electric to capillary forces:

(10)

In our experiments we have ECN<<1, that shows the capillary forces have the highest impact in our system. Finally the interfacial tension is expressed in a following way:

(11)

Using the arguments from the electrowetting, the term for capillary pressure in equation (1) will be changed accordingly. The capillary pressure of two droplets is written as a sum of two pressure components coming from the different interfacial tension:

(12)

3. Energy supplied by electric voltage

The electric field energy is created by applying the voltage to the electrodes. To define the electric response of droplet coalescence system, we used the concepts of electrokinetic model described by Tan et al. This model determines the voltage at the tip of electrode:

(13)

Here i2= -1. CE and Cl are the electric capacitances in the chip which acts as an electric circuit. Next, we we have mathematically removed the imaginary unit from the equation (13).

The cross-sectional view of the droplet coalescence chip is shown in fig.3.

We can rewrite the eq. (13) as

(14)

Then we do a complex-conjugation to get rid of complex number

(15)

After this we can return to the previously-used ratio of voltages:

(16)

From this relation we might find the tip-voltage expressed without an imaginary unit.

From eq. (13) we find the voltage:

(17)

The electrical capacitances are expressed in terms of coalescence chip dimensions:

(18)

The energy density created by an electric field is released into the wider channel (chamber). Multiplying the energy density by a volume of chamber leads to the electric energy supplied by voltage:

(19)

4. Kinetic and dissipated energies

In this section we derive the kinetic energy for the two droplets reaching the coalescence site. We apply the kinetic energy expression, assuming the droplets move in channel with the same velocity as carrier oil. The energy is written as a sum of energies of two different-sized droplets. We write the kinetic energy terms for the energy-balance equation for droplets in a straight channel (HW1) and in the coalescence channel (HWch):

(20)

(21)

Droplets moving through a fluid are exposed to drag forces the same way as solid bodies. The drag force acts against the droplet movement direction. The drag force acting on droplets is expressed by using the Stokes law \(F_{drag} = -6 \pi \eta R v\). Here R is the droplet radius, η and v the viscosity and velocity of the carrier oil. With these arguments the drag-created energy loss is written for two sections analogically to the kinetic energy:

(22)

(23)

It was argued previously that the systems of multiple fluids are dissipative systems. This argument is strengthened for the systems of similar viscosities of two liquids. The coupling between the inner and outer flows around droplets arises and its impact can no longer be neglected [6-7]. The energy dissipated in a coalescence (merging) of droplets will be related to the volume of drained fluid between two droplets T3, \(E_{dissip} = \pi \eta D T^3 (\dot \gamma)^2 t_c\). Here D is the drag coefficient, related to the Reynolds number:

(25)

The strain rate \((\dot \gamma)\) represents the relative deformation of the fluid flow \(\dot \gamma = (Q_1 + Q_2) / H / (W - 2R_2)\). In the experiments we use high concentrations of surfactants in oil and that is why we have to account for the non-Newtonian behavior of the oil. This means the oil viscosity is not constant with the growing deformation rates. It was shown previously that the solutions of high concentration of ionic or nonionic surfactant have a shear thickening behavior [8-10], which means that viscosity of solution is increasing with the growing deformation. We use the Power-Law for viscosity function:

(26)

Here η0 is the viscosity value with no deformation applied. This is the value given in a manual description of oil and n=1.2 is the power law index. With these arguments the dissipated energy is written as follows:

(27)

It should be noted the power-law description of viscosity of oil is used in all energy terms.

5. Energy balance equation and its solutions

So far we derived all the components making influence for energy balance. We write the energy balance as equilibrium of the total energy. The left-side (LS) of equation represents the total energy of droplet pair moving in a straight channel of the cross-section H×W1. Right side (RS) of the balance equation represents the total energy in a coalescence channel of the size H×Wch. These two terms are written as follows

(28)

(29)

Here in equation (29) we include the energy loss associated with the volume of oil drained between two droplets in the coalescence. Taking together equations (28)-(29), the equilibrium equation

(30)

is being solved numerically for various experimental conditions (droplet volumes and voltage). To find the solutions of eq.(30) we used the package Mathematica 9.0.

Solutions of the derived model

We find the solutions of the full balance equation (30) for various volumes of on-chip produced droplets. These solutions are shown in a figure X. As smaller droplets are being reinjected, we keep their volume constant through the range of volumes of produced droplets.

Figure 1. Droplet coalescence time dependence tcoal on fused droplet volumes V1 and V2

As we can see the coalescence time for the single value of V1 has some variation with a maximum. Scanning values from the smallest volume of produced droplet to the largest, there is a range of V2 (the produced droplet) with a longest coalescence time. When we look at the results for the larger droplets V1, the coalescence time increases. Also, we notice that the interval of coalescence times is shorter and variation is smaller. This means that if we are increasing the volume V1, the coalescence time increases, but the interval of possible droplet coalescence becomes smaller. If there’s no solution found, it means no coalescence is occurring for given conditions (droplet volumes, voltage and other).

Figure 2. Droplet coalescence time tcoal dependence on fused droplet volumes at different applied voltage Uapp values.

Next we compared the coalescence time results for various applied voltages Uapp. We chose 3 different values of the applied voltage (0.3V, 0.4V and 0.5V). We have found that even a small increase (0.1V) of voltage strongly affects the coalescence time. With our calculations, the tcoal reduces approximately by factor of 3 as we increase the voltage from 0.3V to 0.5V.

What is even more interesting, is that the interval of possible droplet coalescence does not shorten when we increase the voltage. It becomes shorter only when the volume V1 is increased as we discussed previously.

Having the results for various conditions (V1, V2, Uapp) we now need to determine which of these conditions are suitable for our droplet merging experiments. It is our goal to find a range of droplet volumes which can provide an efficient droplet coalescence.

Firstly, we need to have flow rates tuned appropriately so that the droplets do not pass the whole coalescence chamber without merging with a droplet from another droplet population. For that reason we chose tflow as an upper limit of the time it requires for the droplets to pass through the coalescence chamber. The flow time is expressed as a ratio of channel volume to a unit-volume of fluid (total of the oil):

(31)

If the calculated time tcoal is longer than tflow, droplets will pass the coalescence chamber without fusion. That means if the droplets take a longer time to come together and merge then the time they spend in coalescence chamber, the merging will not occur. That is why we chose the tflow as the higher limit of allowed values.

Next, to find the lower limit must solve the different balance equation. We chose the lower limit as a time needed to remove the double volume of oil between two droplets. In other words, it is the situation where three droplets would merge instead of two. We express the volume of oil between three droplets as \(V_{oil}^{(limit)}=2 T_0 W_1 H\). This volume will be drained at the drainage flow rate. We then write the equation:

(32)

If the equation is True, three droplets will be fused together instead of two. Again, we need to limit our initially found tcoal values. Therefore, if the coalescence time of 3 droplets \(t_{coal}^{(3d)}\) is shorter than that of tcoal, we will find no pair-wise coalescence. The results are shown in Figure X below. The straight lines indicate the upper limit (time of the flow). The curves below are times for 3 droplet fusion. The points that are slightly below the curves for \(t_{coal}^{(3d)}\) represent the regions of stable 2 droplet coalescence with the 3 droplet coalescence time being longer \(t_{coal}^{(3d)} \gt t_{coal}\).

Figure 3. Droplet coalescence time dependence on fused droplet volume ratios at different applied voltage Uapp values. The black line displays tflow - the coalescence time limit, at which droplets leave the merging chamber without fusing. Dotted lines correspond to coalescence time limit tcoal(3d) at which coalescence of 3 droplets occurs. The points that are slightly below the curves for represent the regions of stable 2 droplet coalescence

In conclusion, we find there to be an optimal range of droplet volume ratios V2/V1 that allows stable effective 1:1 coalescence. Aproximately this is the interval from 3 to 5. We also see that this interval gets smaller as we increase the voltage. The interval also becomes smaller with increasing volumes of V1. We have also found that three-droplet merging should be present at the very low and the very high values of V2/V1. We have then experimentally tested these parameters and have successfully reached almost 100% effective droplet merging using volume ratio of 3 to 5. We have then proceeded to use the parameters in our CAT-Seq workflow.

Slowed down video of droplets merging

Physics alphabet
Symbol Description [units]/(sizes)
AH Hamaker constant (10-18J)
a Surfactant monomer size (~2nm)
C Electric capacitance [F] 
Cs Bulk concentration of surfactant [mol/m3]
D Drag coefficient [dimensionless]
d Spacing
E Energy [J]
f Frequency of AC-current [Hz]
H Height of the chip [m]
kads Surfactant adsorption constant (kads=18mol-1s-1)
kdes Surfactant desorption constant (kdes=6·10-3s-1)
l Length of the channel [m]
N Number of droplets  
n Power-law index for non-Newtonian viscosity (n=1.2)
Q Flow rate of fluid [m3/s]
R1,2 Radius of droplet [m]
R Universal gas constant (R=8.314 J/mol/K)
T Absolute temperature [K]
T Thickness of droplet spacing-oil [m]
t Time [s]
U Electric voltage [V]
V Volume of droplet [m3]
γ Interfacial tension [N/m]
ε0 Dielectric constant (ε0=8.85·10-12F/m)
ε Relative dielectric permittivity (εPDMS=2.5, εOil=5.8, εGlass=7.5)
η Dynamic viscosity [Pa·s]
θ Contact angle  
κ electric conductivity [S/m]
ρ Material density [kg/m3]
References
  1. Barman , Nagarajan AK, Khare K., RSC Adv.2015, 5.
  2. C. Baret, F. Kleinschmidt, A. El Harrak, and A. D. Griffiths, Langmuir 2009, 25(11),6088–6093.
  3. Brosseau , Vrignon J., Baret JC., Soft Matter, 2014, 10, 3066-3076.
  4. Riechers B., Maes F., Akoury E., Semin B., Gruner P., Baret JC., PNAS October 11, 2016 113 (41), 11465-11470.
  5. Tan SH, Semin B, Baret JC, Lab Chip, 2014, 14, 1099–1106.
  6. Baroud CN, Gallaire F., Dangla R., Lab Chip, 2010, 10, 2032–2045.
  7. Eggers J., Villermaux E. Prog. Phys. 71 (2008) 036601 (79pp).
  8. Shu-peng, Journal of hydrodynamics, 2012,24(2):202-206.
  9. Hellsten, Journal of Surfactants and Detergents, Vol. 5, No. 1 (January 2002).
  10. A. Khadom, A. A. Abdul-Hadi, Ain Shams Engineering Journal (2014) 5, 861–865.