Open Access Article

This Open Access Article is licensed under a Creative Commons Attribution-Non Commercial 3.0 Unported Licence

DOI: 10.1039/D0SM01113D
(Paper)
Soft Matter, 2020, Advance Article

François Villemot,
Antoine Calmettes and
Marc Durand*

Université de Paris, CNRS, UMR 7057, Matière et Systèmes Complexes (MSC), F-75006 Paris, France. E-mail: marc.durand@univ-paris-diderot.fr

Received
16th June 2020
, Accepted 8th September 2020

First published on 10th September 2020

Analysis of thermal capillary waves on the surface of a liquid usually assumes incompressibility of the bulk fluid. However, for droplets or bubbles with submicronic size, or for epithelial cells whose out-of-plane elongation can be modeled by an effective 2D bulk modulus, compressibility of the internal fluid must be taken into account for the characterization of their shape fluctuations. We present a theoretical analysis of the fluctuations of a two-dimensional compressible droplet. Analytical expressions for area, perimeter and energy fluctuations are derived and compared with Cellular Potts Model (CPM) simulations. This comparison shows a very good agreement between theory and simulations, and offers a precise calibration method for CPM simulations.

Thermal capillary waves have been the subject of intense research activity, both theoretically and experimentally,

With this aim in mind, we investigate in this paper the thermal fluctuations of a “two-dimensional” compressible droplet with low (but finite) compressibility; this may correspond to a droplet squeezed between two parallel plates. Analytical expressions for the average and variance of area, perimeter and energy are derived and compared with Cellular Potts Model (CPM) simulations, which is a standard numerical modelling of soft cellular systems. The present article is organized as follows. In Section 1, we evaluate the partition function of a 2D droplet with low but finite compressibility. The calculations are more complicated than in the incompressible case because the compressing energy introduces quartic terms of the contour position vector which cannot be neglected. In Sections 2 and 3 we derive analytical expressions for the average and standard deviations of the area, contour length, and energy, which fully characterize their normal probability distributions. In Section 4 we calculate the mean square amplitude of the contour fluctuation modes. Then, in Section 5 we confront the different expressions with numerical simulations performed with CPM simulations. The agreement between theory and numerical data is very good, and offers a robust method for calibrating surface tension in CPM. In particular, we explain in Section 6 that we are able to detect a slight thermal softening of CPM surface tension which is caused by shorter wavelength modes.

= _{b} + _{c}.
| (1) |

(2) |

To describe the vesicle contour, it is convenient to use the polar parametrization θ ↦ r_{c}(θ), as shown in Fig. 1, where the location of the center O is chosen arbitrarily for the moment. The perimeter and area of the droplet can be expressed as functionals of r_{c}(θ):

(3) |

The relative importance of boundary, compression, and thermal energies are characterized by two dimensionless quantities: the relative compressibility ε = γ/Br_{0} already introduced above, and the relative amplitude of fluctuations η = k_{B}T/γr_{0}, where . We consider here the low compressibility (ε ≪ 1) and low fluctuations (η ≪ 1) regime. In order to derive an analytical expression for the partition function, we need to expand the Hamiltonian around a reference configuration that we determine below.

(4) |

^{3} − + ε = 0.
| (5) |

(6) |

Expansion of the right-hand-side term in the low compressibility regime yields:

(7) |

(8) |

(9) |

(10) |

(11) |

(12) |

(13) |

(14) |

The Hamiltonian rewrites:

(15) |

(16) |

(17) |

(18) |

(19) |

(20) |

(21) |

(22) |

(23) |

(24) |

(25) |

Expressing the Hamiltonian (eqn (17)) in terms of the variables L,{ũ_{p≥1}} instead, one finds that the perimeter droplet also follows a Gaussian distribution: P(L) ∼ e^{−(L−L)2/(2ΔL2)}. The average perimeter can be easily expressed using the identity (considering A^{⋆}, L^{⋆}, r^{⋆} as constants). It comes:

(26) |

(27) |

(28) |

P(E) ∝ g(E)e^{−βE},
| (29) |

(30) |

g(E)dE ∝ Ω_{d}R^{d−1}dR,
| (31) |

P(E) = c(E − E^{⋆})^{N/2}e^{−βE},
| (32) |

(33) |

(34) |

(35) |

(36) |

Finally, the mean amplitude of mode p = 0, which corresponds to a uniform compression of the droplet, is obtained from eqn (16): as 〈A〉 = A^{⋆} (eqn (24)), and using the mathematical identity , it comes:

(37) |

The boundary energy term in CPM simulations is defined as:

(38) |

Definition of length on lattice-based models such as CPM must be done cautiously, as it is a common source of confusion and implementation errors.^{19} Here we define the length of the boundary from its energy. We consider that site k belongs to the boundary as soon as there is a site in its coupling neighborhood whose label is different from σ_{k}. Therefore, a boundary in CPM has a thickness which is directly related to the size of coupling neighborhood. Let us introduce z, the number of sites with ID ≠ σ_{k} in the coupling neighborhood of a boundary site k, averaged over all sites k that belong to that boundary. We then define the boundary length as , i.e. the boundary energy measured for a unit coupling constant, divided by the boundary thickness z. Accordingly, surface tension is given by γ = zJ (in units of energy per pixel^{2}). Roughly, z is equal to half of the number of sites in the coupling neighborhood. Actually, because of the anisotropy of the underlying lattice, z – and so γ – slightly depend on the contour shape: their values will be different for a square and a circular contour. Values of z reported in the literature range from 10.5 to 11.3^{19–21} for 4th order coupling neighborhood, which has 20 neighboring lattice sites (see Fig. 2). Here we propose to take advantage of the thermal fluctuations of the contour and perform statistical averaging to lessen this dependency on contour shape, and determine an accurate value for z.

Fig. 2 (a) Order I neighborhood, with 4 lattice sites (in orange), is used as target neighborhood in our simulations, in which the target value is chosen from in the updating rule. (b) Order IV neighborhood, with 20 lattice sites (in orange), is used as coupling neighborhood _{c} for computing the boundary energy (eqn (38)). |

The coupling neighborhood _{c}, used for calculation of the boundary energy is the order IV neighborhood, with 20 neighbors for each pixel. The target neighborhood, in which the target value is chosen from in the site updating rule of the Metropolis-like algorithm, is the order I neighborhood, with 4 neighbors per pixel (see Fig. 2).

We run simulations for different points in the phase space {J,B,T,A_{0}}, in a way that the two dimensionless parameters ε and η are varied independently. We ensure that values of variables are in ranges such that the hypotheses we made in our theoretical model stay valid in simulations. Namely, we ensured weak bulk compressibility (ε ≤ 0.1), small thermal fluctuations (η ≤ 0.1), and a sufficient number of degrees of freedom for the membrane (a/r_{0} ≪ 1), where here a and z are taken to be of order 1 and 10 pixels, respectively. Each simulation was run for 10^{6} Monte-Carlo steps, and statistics are performed over 2 × 10^{−4} measures. The number of uncorrelated samples is obtained using the pymbar.timeseries module,^{23} so that uncertainties of the measured quantities are based on fully independent samples.

We want to use the comparison with theory to determine the values of the two unknown parameters z and a. We then introduce the uncalibrated/raw length L′ = _{b}/J, the calibrated length being L = L′/z, and we rewrite eqn (24), (25), (27), (28), (33) and (34) in terms of the dimensionless variables ε′ = J/Br_{0} and η′ = k_{B}T/Jr_{0}, rather than ε = zε′, η = η′/z. Using expansion eqn (7) for r^{⋆}, we obtain:

(39) |

(40) |

(41) |

(42) |

(43) |

(44) |

Fig. 3 shows the typical histograms of A, L′ and E obtained for one simulation run. In agreement with our theoretical treatment, they all follow Gaussian distributions. We then report the variation of different quantities as function of the dimensionless variables ε′ and η′, and use expressions eqn (39)–(44) as fitting functions, where a and z are the two adjusting parameters. Fit values of z and a are summarized in Fig. 4.

Fig. 3 Histograms of perimeter L, area A and energy E. Orange curve are Gaussian functions whose mean and variance are calculated from the corresponding histogram. |

Fig. 4 Fit values for z and a obtained using expressions (39)–(44) with numerical data. Values reported with open symbols were obtained by setting the value of either z or a and using only the other one as a fitting parameter. |

Fig. 5 shows the variation of 〈A〉/A_{0} and (ΔA/A_{0})^{2} with ε′ and η′. In agreement with eqn (39), we observe in Fig. 5a clear linear dependence of 〈A〉/A_{0} on ε′ whose slope provides a first estimation for z. However, we observe a slight dependence of the slope on η′, which is not predicted by eqn (39). This suggests that the CPM calibration parameter z slightly decreases as η′ increases. This weak dependence can also be observed in Fig. 5b: a linear fit of the data for a given ε′ value yields a small positive slope, indicating that z – and hence, the surface tension in CPM simulations – is a slightly decreasing function of η′. We discuss in Section 6 of the origin of this thermal softening of surface tension. Expression of the area variance (eqn (40)) does not depend on z and a, and the comparison with numerics in Fig. 5c and d is done with no adjusting parameter. Nevertheless, the agreement is very good.

Fig. 5 (a) Variation of 〈A〉/A_{0} with dimensionless parameter ε′. Solid curves are fits using eqn (39) with adjusting parameter z. (b) Variation of 〈A〉/A_{0} with dimensionless parameter η′. Solid curves are averages. Dashed curves are linear fits. (c and d) Variation of (ΔA/A_{0})^{2} with dimensionless parameters ε′ and η′. Solid curves correspond to eqn (40) with no adjusting parameters. |

Fig. 6 shows the variation of 〈L′〉/L_{0} and (ΔL/L_{0})^{2} with ε′ and η′. Fits are in excellent agreement with eqn (41) and (42), giving credence to our theoretical treatment. For the variation of 〈L〉/L_{0} with η′, z derived from the fits are in the expected range. Values are slightly above those deduced from the mean area, because the decrease of the proportionality coefficient between γ and J at short wavelength has different impact on the coarse-grained curvature and perimeter of the droplet. Fits in Fig. 6a, c, and d do not allow to determine precisely both z and a, as a whole range of pairs of parameter values gives similar data fit. For these curves, the value of z was set to the one obtained from fitting the variation of 〈L〉/L_{0} with η′, and a was left as the only adjusting parameter. Values obtained for a from are in a narrow range (a ∼ 1.1), and in agreement with our expectation: a step of a = 1 pixel corresponds to a boundary aligned with one of the square lattice axes. The fact that values are slightly above 1 can be attributed to the finite curvature 1/r^{⋆} of the droplet, and to the contribution of other boundary orientations in the averaging. Values of a obtained from the variation of (ΔL/L_{0})^{2} with ε′ show a slight dependence with η′, and are more spread out than other values. This dependence can be explained from the variation of z with η′ revealed earlier, which has been ignored in the fits, as they were performed with the same value for z.

Fig. 6 (a and b) Variation of 〈L〉/L_{0} with dimensionless parameters ε′ and η′. Solid curves are fits using eqn (41) with adjusting parameters z and a. (c and d) Variation of (ΔL/L_{0})^{2} with dimensionless parameters ε′ and η′. Solid curves correspond to eqn (42) with adjusting parameter a. |

We now test our predictions on statistical distribution of energy (eqn (43) and (44)). Variations of 〈E〉/JL_{0} and (ΔE/JL_{0})^{2} are shown in Fig. 7. Again, we have perfect agreement between the numerical data and the fits from the analytical expressions. The values for z and a obtained from these fits are almost identical to the ones obtained from the statistical distribution of perimeter. The energy is a function of perimeter and surface area, and since the boundary energy is much larger than the bulk compressive energy, it comes as no surprise that we observe similar trends for the energy and the perimeter.

Fig. 7 (a and b) Variation of 〈E〉/JL_{0} with dimensionless parameters ε′ and η′. Solid curves are fits using eqn (43) with adjusting parameter z and a. (c and d) Variation of (ΔE/JL_{0})^{2} with dimensionless parameters ε′ and η′. Solid curves correspond to eqn (44) with adjusting parameter a. |

Finally, we test our predictions on the amplitudes of modes (eqn (35) and (37)). From the positions r_{c}(θ) of lattice sites localized at the droplet boundary for a given configuration, we deduce u(θ) = r_{c}(θ)/r^{⋆} − 1, where , and then evaluate ũ_{p} from eqn (10). Fig. 8 shows the variation of 〈|ũ_{p}^{2}|〉 with p. Again, data are well predicted by the theoretical expression eqn (35), and yields for the unique fit parameter z = 10.88, consistent with values obtained from other expressions (see Fig. 4). We also reported on the graph −4〈ũ_{0}〉/9 calculated with this value of z. In agreement with eqn (37), its value is very close to 〈|ũ_{2}|〉^{2}.

Fig. 8 Mean square amplitude of fluctuation modes ũ_{p}. Orange curve is a fit using eqn (35) with z as adjusting parameter. Fit value is z = 10.88. The simulation was done with parameters A_{0} = 5000, ε′ = 3 × 10^{−4}, and η′ = 0.03. |

In future extensions to this work, we will consider the case of compressible vesicles by including bending and stretching terms to the boundary energy. We will also extend our analysis to multicellular systems, for which contour fluctuations of every cellular unit are related in a sophisticated way, that involves in particular the topology of the cellular pattern.

(45) |

(46) |

The simplest expression for the boundary term reads _{b} = γ; _{b} does not depend on the specific shape of the droplet either, but only on its perimeter . Note that this expression implies that surface tension γ is uniform; if the droplet boundary contains surfactant molecules, it is then assumed that fluctuations of their surface concentration can be neglected.

(47) |

(48) |

(49) |

(50) |

(51) |

(52) |

(53) |

(54) |

(55) |

(56) |

Thus, ℑ_{ij} is triangular with one diagonal element equal to 1/2A^{⋆}, and the others equals to 1. Therefore, |ℑ| = 1/2A^{⋆}.

The list of unique labels amongst the neighbors is first recorded. If multiple neighboring sites have the same label, this label only appears once in the list. The target label is then randomly selected in that list, with an even probability. This setup ends up being different than the one in the standard CC3D code, as if a label appears multiple times amongst the neighboring sites, it will not be more likely to be chosen than if it only appears once.

- A. K. Doerr, M. Tolan, W. Prange, J.-P. Schlomka, T. Seydel, W. Press, D. Smilgies and B. Struth, Phys. Rev. Lett., 1999, 83, 3470–3473 CrossRef CAS.
- J. Rowlinson and B. Widom, Molecular Theory of Capillarity, Dover Publications, 2002 Search PubMed.
- D. G. A. L. Aarts, M. Schmidt and H. N. W. Lekkerkerker, Science, 2004, 304, 847–850 CrossRef CAS.
- D. Langevin, Colloids Surf., 1990, 43, 121–131 CrossRef CAS.
- S. T. Milner and S. A. Safran, Phys. Rev. A: At., Mol., Opt. Phys., 1987, 36, 4371–4379 CrossRef CAS.
- H. Gang, A. H. Krall and D. A. Weitz, Phys. Rev. Lett., 1994, 73, 3435–3438 CrossRef CAS.
- W. Cai, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1996, 54, 2780–2785 CrossRef CAS.
- E. Hannezo, J. Prost and J.-F. Joanny, Proc. Natl. Acad. Sci. U. S. A., 2014, 111, 27–32 CrossRef CAS.
- F. Graner and J. A. Glazier, Phys. Rev. Lett., 1992, 69, 2013–2016 CrossRef CAS.
- D. B. Staple, R. Farhadifar, J. C. Röper, B. Aigouy, S. Eaton and F. Jülicher, Eur. Phys. J. E: Soft Matter Biol. Phys., 2010, 33, 117–127 CrossRef CAS.
- D. Bi, J. H. Lopez, J. M. Schwarz and M. L. Manning, Nat. Phys., 2015, 11, 1074–1079 Search PubMed.
- M. Durand and J. Heu, Phys. Rev. Lett., 2019, 123, 188001 CrossRef CAS.
- S. Ljunggren and J. C. Eriksson, J. Chem. Soc., Faraday Trans. 2, 1984, 80, 489–497 RSC.
- M. Durand and E. Guesnet, Comput. Phys. Commun., 2016, 208, 54–63 CrossRef CAS.
- J. A. Glazier, M. P. Anderson and G. S. Grest, Philos. Mag. B, 1990, 62, 615–645 CAS.
- Y. Jiang, P. J. Swart, A. Saxena, M. Asipauskas and J. A. Glazier, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1999, 59, 5819 CrossRef CAS.
- A. J. Kabla, J. R. Soc., Interface, 2012, 9, 3268–3278 CrossRef.
- T. Hirashima, E. G. Rens and R. M. H. Merks, Dev., Growth Differ., 2017, 59, 329–339 CrossRef.
- R. Magno, V. A. Grieneisen and A. F. Marée, BMC Biophys., 2015, 8(8) DOI:10.1186/s13628-015-0022-x.
- J. Käfer, T. Hayashi, A. F. Marée, R. W. Carthew and F. Graner, Proc. Natl. Acad. Sci. U. S. A., 2007, 104, 18549–18554 CrossRef.
- P. Marmottant, A. Mgharbel, J. Käfer, B. Audren, J.-P. Rieu, J.-C. Vial, B. Van Der Sanden, A. F. Marée, F. Graner and H. Delanoë-Ayari, Proc. Natl. Acad. Sci. U. S. A., 2009, 106, 17271–17275 CrossRef CAS.
- M. H. Swat, G. L. Thomas, J. M. Belmonte, A. Shirinifard, D. Hmeljak and J. A. Glazier, Computational Methods in Cell Biology, Academic Press, 2012, vol. 110, pp. 325–366 Search PubMed.
- J. D. Chodera, J. Chem. Theory Comput., 2016, 12, 1799–1805 CrossRef CAS.
- O. Farago and P. Pincus, Eur. Phys. J. E: Soft Matter Biol. Phys., 2003, 11, 399–408 CrossRef CAS.

This journal is © The Royal Society of Chemistry 2020 |