Gavin
Melaugh
*,
Davide
Marenduzzo
,
Alexander
Morozov
and
Rosalind J.
Allen
SUPA, School of Physics and Astronomy, University of Edinburgh, Edinburgh, EH9 3FD, UK. E-mail: g.melaugh@ed.ac.uk; Tel: +44 (0)131 651 3456
First published on 23rd September 2019
Mechanical interactions between biological cells can be mediated by secreted products. Here, we investigate how such a scenario could affect the cells' collective behaviour. We show that if the concentration field of secreted products around a cell can be considered to be in steady state, this scenario can be mapped onto an effective attractive interaction that depends on the local cell density. Using a field-theory approach, this density-dependent attraction gives rise to a cubic term in the Landau–Ginzburg free energy density. In continuum field simulations this can lead to “nucleation-like” appearance of homogeneous clusters in the spinodal phase separation regime. Implementing the density-dependent cohesive attraction in Brownian dynamics simulations of a particle-based model gives rise to similar “spinodal nucleation” phase separation behaviour.
A system of living cells that actively secrete a cohesion-influencing chemical agent (be it a polymer or signal molecule) is clearly out of equilibrium. Yet, even for out-of-equilibrium living organisms, useful mappings can often be made to key concepts in soft matter and statistical physics.^{18–27} Importantly, the case we investigate here, in which the cohesion-inducing agent is actively produced, is quite different to “passive” scenarios in which a cohesion-inducing agent (e.g., a polymer) is added to a system of cells or colloidal particles.^{24,28–32} In the latter cases, the polymer concentration is conserved; in the case that we consider here, the cohesive agent is constantly produced and degraded, and hence is not conserved.
In this paper, we investigate the collective behaviour of a system of aggregating particles secreting a cohesion-inducing agent. We follow previous studies in active matter physics,^{19,21–24,26} by considering aggregation as a phase separation process, albeit here our focus is on cell–cell cohesion and not activity arising from cellular motility. First, we formulate a mathematical model that shows that a system of particles producing a cohesion-inducing agent under steady-state conditions can be coarse-grained via an effective attractive interaction that depends on the local particle density. Taking a field-theory approach, this density dependence gives rise to a cubic term in the Landau–Ginzburg free energy that alters the thermodynamic landscape in comparison to that of a standard, density-independent model. Computer simulations of the Cahn–Hilliard equation also show different dynamics during early-stage phase separation for the density-dependent versus density-independent systems. Specifically, for the density-dependent system we observe “nucleation-like” appearance of homogeneous particle clusters in the spinodal decomposition regime. Similar emergent behaviour arises when the density-dependent cohesion is implemented in particle-based Brownian Dynamics simulations.
(1) |
−D∂_{r′}c|_{r′=δ} = w, | (2) |
(3) |
Fig. 1 Schematic illustration of cohesion-inducing agents (red disks) emanating from a source particle of radius δ (green). At the particle's surface, the flux of the cohesion-inducing agent is fixed according to eqn (2); this boundary condition imposes a diffusive flux of cohesive agent away from the source. |
Consider now a collection of N such sources located at positions (r_{1}, r_{2},…r_{N}) in a volume V. The total concentration of the cohesion-inducing agent at a location R, sufficiently far away from the sources, can be approximated by
(4) |
(5) |
We now take a further step and assume that the cohesion-inducing agents can be treated in a coarse-grained manner, as an effective interaction between the particles. We suppose that the agents produce attractive forces between the particles, and that the degree of attraction is proportional to the local concentration of cohesive agents. Following our analysis above (in particular eqn (5)), this in turn implies that the degree of attraction depends on the local density of the particles themselves. Hence, our system can be modelled as a collection of particles interacting via a density-dependent attractive potential.
In soft matter physics, the concept of a density-dependent potential is not new: for example, density-dependent effective pair potentials arise when treating polymers as soft colloids, and when treating a two-component mixture as an effective one-component system, as in the Asakura-Oosawa model for polymer–colloid mixtures.^{34} It is well known that naive application of potentials that depend on the global particle density can lead to problems such as a lack of consistency between virial and compressibility routes to the equation of state.^{34} This problem appears to be alleviated in the case of interactions that depend on the local density, like those we consider here.^{34,35} Moreover, our model represents a non-equilibrium system, in which cohesive agents are continuously produced – for such a system we would not necessarily expect the rules of equilibrium thermodynamics to be obeyed.
We also note that in a real system of interacting biological cells, many other factors may be at play, including motility and cell growth and division. For the sake of simplicity, these factors are neglected here.
(6) |
In Landau–Ginzburg theory, the phenomenological coefficients a and c often depend on the temperature. The coefficient a, which is a function of the deviation of the system temperature and the critical temperature (T − T_{c}), is the critical parameter that governs the transition between the disordered and ordered states. In our model, however, we will not consider this temperature dependence explicitly. Rather, in our Landau-like free energy density, eqn (6), we consider the ξ^{2} term as an effective two-body interaction in which the critical parameter a governs the transition between the disordered state in which entropic/repulsive interactions dominate, and the ordered state in which enthalpic/cohesive interactions dominate (Fig. 2(a)). Thus in our model a encapsulates a tradeoff between repulsive (ψ) and cohesive (υ) contributions, such that a = υ − ψ (analogous to the quadratic coefficient in the Bragg–Williams approximation for the Ising model^{37}). We fix ψ arbitrarily to −0.05 so that the phase transition is governed by the value of υ, which represents the strength of the cohesive interactions. The quartic term in eqn (6) ensures that the equilibrium state has a bounded value of ξ, and thus must have a positive coefficient c > 0. The square gradient term imposes a free energy cost for any non-uniformity in ξ, and thus κ, which is related to the surface tension, must be positive. Note that we omit a positive cubic term in eqn (6), which is allowed by symmetry but would not affect our results, to facilitate numerical comparison to the density-dependent case described below.
Fig. 2 Homogeneous part of the Landau–Ginzburg free energy density and the associated phase diagrams. (a) Free energy density, f(ξ), for various values of the critical parameter, a, in the density-independent system (eqn (6)). Note that the point ξ = 0 undergoes a continuous transition from a minimum to a maximum as a decreases through a = 0. Thus, this free energy density would produce a second-order phase transition in a system where the order parameter was not globally conserved. (b) Free energy density, f(ξ), for various values of the critical parameter, b, in the density-dependent system (eqn (11)). Note the emergence of a secondary minimum at ξ ≥ 0 as b decreases. Upon decreasing through b ∼ −0.237, where two equal free energy minima appear, the system would undergo a first-order phase transition if the order parameter was not globally conserved. (c and d) Phase diagrams generated by applying the common-tangent construction to the density-independent (c) and density-dependent (d) free energy profiles (see the ESI† for details). The black solid and dashed lines denote the binodal and spinodal coexistence lines respectively. Red arrows represent constant density quenches to the points (ξ = 0.07, a = −1.0, and ξ = 0.07, b = −0.5) within the coexistence regions of (c) and (d) respectively, as discussed in the text. Green arrows represent constant density quenches to the points (ξ = 0.3, a = −1.0, and ξ = 0.3, b = −0.5) within the coexistence regions of (c) and (d) respectively. |
For positive values of the parameter a, the free energy density, eqn (6), has a single minimum at ξ = 0. In contrast, for negative values of a, eqn (6) has two minima at . For the model represented by eqn (6) to be physically realistic, the condition of positive density, ξ ≥ 0, has to be imposed throughout, ruling out any equilibrium phase with ξ < 0. Fig. 2(a) shows the homogeneous part of the free energy density (eqn (6)); the region shaded red is inaccessible because it does not meet the condition ξ ≥ 0. The continuous transition from a state whereby ξ = 0 is the global minimum of the free energy density to a state where the global minimum is a positive finite value, , is a signature of the fact that this free energy density would produce a second-order phase transition in a system where the order parameter is not globally conserved.
υ ∝ c(x) = Kξ, | (7) |
(8) |
(9) |
(10) |
(11) |
In our model, particles are not created or removed (since we neglect cell growth or death), and therefore the density must be globally conserved. To minimise the free-energy density subject to this constraint, we apply the common-tangent construction to the homogeneous part of the free energy densities, eqn (6) and (11), to generate the density-independent and density-dependent phase diagrams shown in Fig. 2(c) and (d) respectively (see the ESI† for details). Importantly, the cubic term in the density-dependent system gives rise to a region of metastability (black dashed line) at small values of ξ that is absent in the density-independent system.
(12) |
Fig. 4 The time evolution of the density-independent (a–c) and -dependent (d–f) systems following a quench highlights the effect of the free energy barrier at low density shown in Fig. 3. (a) Density-independent, ξ_{0} = 0.07, a = −1.0, at t = 0 SU. (b) Density-independent, ξ_{0} = 0.07, a = −1.0, at t = 250 SU. (c) Density-independent, ξ_{0} = 0.07, a = −1.0, at t = 47500 SU. (d) Density-independent, ξ_{0} = 0.07, b = −0.5, at t = 0 SU. (e) Density-independent, ξ_{0} = 0.07, b = −0.5, at t = 250 SU. (f) Density-independent, ξ_{0} = 0.07, b = −0.5, at t = 47500 SU. |
We also performed a similar quench for the density-dependent case. Here the parameter b was decreased as shown by the red arrow in Fig. 2(d) (b > 0 → b = −0.5), such that the resulting condensed phase density, ξ_{cp} ∼ 2, was similar to that of the density-independent system. In this case, the system is quenched into the low-density metastable region of the phase diagram (which does not exist for the density-independent case), and therefore we expect phase separation to occur via nucleation, requiring a large enough density fluctuation to overcome the free-energy barrier (red curve Fig. 3(a)). Indeed, the snapshots in Fig. 4(d)–(f) are consistent with nucleation: we observe the formation of a single aggregate that grows larger with time (see also Movie S2 in ESI†).
We would like to know whether phase separation in the spinodal region differs intrinsically between the density-dependent and density-independent scenarios. To test this, we performed a series of quenches in the density-independent case, for ξ = 0.3, a > 0 initially, and a range of final a values (−0.1,−0.2, −0.3,…,−1.0; dark green points along the green line in Fig. 2(c)). The rationale for performing a range of simulations was that the density-dependent simulation effectively samples a range of υ values (since υ = 2/3υ_{0}Kξ); thus many values of a should be sampled in the density-independent system in order to compare with the density-dependent case. Fig. 5(b) and (c) show snapshots of these simulations, for a = −0.2 and a = −1.0 respectively, taken at time t = 2.5 × 10^{5} SU following the quench. Clearly, the spinodal decomposition process is affected by the value of a: a more negative value of a produces a more rapid contraction of the early-stage elongated structures into the later-stage rounded clusters (see also Movie S4 (a = −0.2) and S5 (a = −1.0) in ESI†). The coloured lines in Fig. 5(a) show the distribution of density values at this same time point, for each of the density-independent simulations. These distributions are also bimodal, with a gas-phase peak at ξ ∼ 0.07 and a condensed-phase peak at a value ξ_{cp} that increases with decreasing a.
Comparing the density distributions in Fig. 5(a) between the density-dependent and density-independent cases (black vs. coloured lines), we see a notable difference: aggregate formation appears to be hindered in the density-dependent system. Specifically for the density-dependent case, the gas phase peak is larger, and the condensed phase peak is smaller. This points to a higher fraction of the gas phase coexisting with the condensed-phase aggregates. This affects the aggregate size distribution in the condensed phase: comparing for example Fig. 5(c) and (d), the aggregates in the density-dependent simulation (Fig. 5(d)) appear smaller and more narrowly distributed in size. This is confirmed by the aggregate size distributions shown in Fig. 6. We also observe that density-dependence gives rise to aggregates with more diffuse interfaces. These results suggest that the existence of the low-density metastable region (Fig. 2(d)) can affect the phase separation dynamics even upon quenching into the spinodal region. Specifically, local density fluctuations can bring the system locally into the low-density metastable region, where there is a free energy barrier to phase separation.
Fig. 6 Field simulations. Aggregate size distributions for the density-dependent (red) and -independent (blue) systems computed at t_{run} = 2.5 × 10^{5} SU. The distributions are normalised such that they represent the probability that a given pixel belongs to an aggregate of size A. All distributions were generated from 5 final configurations at t = 2.5 × 10^{5} SU, produced in independent replicate simulations. The construction of these distributions is described in more detail in the ESI.† |
(13) |
Fig. 3(b) shows that this curvature is significantly less negative in the density-dependent case than the density-independent case (comparing (ξ = 0.3, b = −0.5) and (ξ = 0.3, a = −1.0)). We therefore expect spinodal decomposition to happen more slowly for the density-dependent case. Indeed, Fig. 7(a) shows that this is the case; here we plot the maximum density ξ_{max} in our simulations as a function of time following the quench to (ξ = 0.3, b = −0.5) or (ξ = 0.3, a = −1.0) for the density-dependent and -independent cases respectively. We can also measure the time taken to reach the condensed phase density ξ_{cp}, which is ∼1.75 in both cases: this time is indeed much shorter for the density-independent simulation (t_{ps} ∼ 550 SU versus t_{ps} ∼ 3750 SU in the density-dependent case); see also Movies S6 (density-independent) and S7 (density-dependent) in ESI.† From Fig. 7(b), which shows representative density distributions at t_{ps} of the density-dependent (∼3750 SU) and -independent systems (∼550 SU), we see that the density-dependent system contains a much larger proportion of the noncondensed phase at this early stage.
Fig. 7 Field simulations. Phase separation in the density-dependent system (red curves) is slower than in the density-independent system (blue curves). (a) Time evolution of the maximum density, ξ_{max}. The dashed horizontal line corresponds to the phase separation density ξ_{cp} = 1.75, the location of the condensed phase peak in Fig. 5(a). Black vertical lines represent the time, t_{ps}, that it takes for the systems to reach ξ_{cp}. (b) The corresponding distribution of densities, p(ξ), at t_{ps} ∼ 3750 SU (density-dependence) and t_{ps} ∼ 550 SU (density-independence). The black vertical line corresponds to ξ_{cp} = 1.75. |
(14) |
To test this, we performed density-dependent simulations initialised with the configurations from density-independent simulations taken at t_{ps} (i.e., once the condensed phase density had been reached), and vice versa. Fig. 9 shows representative snapshots during the time evolution of these “switched” simulations (2nd and 4th rows), alongside, for comparison, representative snapshots from simulations (initiated from the same configurations) that were not switched (1st and 3rd rows). Rows 1 and 2 illustrate the effect of turning on density-dependent interactions in systems initialised from density-independent simulations. Shortly after switching (compare (b) and (e)), the density becomes redistributed in a manner that increases the concentration of the gas phase whilst lowering that of the dense phase (see also Fig. S13, ESI†). This redistribution persists at later times (see Fig. S13, ESI†), leading to smaller aggregates (compare (c) and (f)) that are more narrowly distributed in size (see Fig. S14, ESI†). In other words, the system takes on the characteristics of a density-dependent free-energy function.
Switching from density-dependent interactions to density-independent interactions (rows 3 and 4) has the opposite effect, i.e., the aggregates become less diffuse and the density is redistributed to decrease the gas-phase concentration whilst increasing that of the condensed phase aggregates; this behaviour persists at later times, leading to large aggregates that are more broadly distributed in size – i.e., the system takes on the characteristics of a density-independent free-energy function. Therefore we conclude that it is the underlying thermodynamics which govern the phase separation process; the trajectory history has little influence on the long term dynamics of phase-separation.
U(r,ρ) = 4ε(ρ)[(σ/r)^{12} − (σ/r)^{6} − U_{c}], | (15) |
Following eqn (5), ε(ρ) is assumed to take the linear form
ε = ρ(x)ε′, | (16) |
In this particle-based model, as in our field simulations, the total density is conserved, i.e., particles do not grow, reproduce or die. The particles are also assumed to be non-motile.
The position, x_{i} of an individual particle, i, evolves in our simulations via numerical integration of the over-damped Langevin equation
(17) |
To characterise the system with density-dependent interactions, we systematically varied the interaction strength parameter, ε′, for three values of the total number of particles: N = 3600, corresponding to area fraction θ = 0.21, N = 4900 (θ = 0.29), and N = 6400 (θ = 0.38).^{47} We observed similar phase separation behaviour in all three systems – we therefore focus on the results for the intermediate area fraction θ = 0.29 (N = 4900). Results for the other area fractions can be found in ESI.†
After an equilibration period (1000 τ) without inter-particle attractions, our simulations were run for times t_{run} = 1250 τ (in total 2250 τ), which is the approximate time required for the fraction of particles in the condensed phase to reach a steady state. After time t_{run}, coarsening and coalescence change the distribution of aggregate sizes (ultimately leading to the formation of one large phase-separated domain), but the aggregated fraction remains constant (see the ESI† for details). The time-scales used here, although not long enough to observe fully phase separated states, are long enough to capture the rich early aggregation behaviour resulting from density-dependent cohesion.
Fig. 10 shows configurations of the system for θ = 0.29 and for various values of ε′, after simulation time t_{run} = 1250τ. For ε′ = 1 (weak density-dependent attraction), the system is in a “gas-like” phase (Fig. 10(a)) whose ordering, characterised by the structure factor S(q), is consistent with that of a simple colloidal dispersion or hard sphere fluid (see the ESI† for details). For ε′ = 70 (strong density-dependent attraction), elongated and interconnected aggregates give rise to a “gel-like” state (Fig. 10(c)).
For ε′ = 40 (intermediate density-dependent attraction), we obtain interesting behaviour arising from the density-dependent attraction. Our simulation snapshots at t_{run} = 1250τ (Fig. 10(b)) show clearly that the system is undergoing phase separation into condensed and non-condensed phases, with the emergence of the former being confirmed by the presence of well-defined peaks in S(q) (see the ESI† for details). For this value of ε′, inspection of the simulation trajectory shows that phase separation proceeds via the formation of aggregates that appear to grow only once they reach a threshold size (see Movie S8 in ESI†). This is reminiscent of the “spinodal nucleation” phenomenon observed in our field simulations. Fig. 11(a) shows the aggregate size distribution for ε′ = 40, computed at time t_{run} (i.e., the probability that a particle belongs to an aggregate of a given size). Consistent with the snapshot of Fig. 10(b), the distribution is bimodal, with two peaks corresponding to aggregates of >300 particles and individual particles (clusters of <2 particles).
Fig. 11 Particle-based simulations. Aggregate size distributions, computed at t_{run} = 1250τ, for density-dependent and -independent systems with area fraction θ = 0.29. The distributions are weighted by the aggregate size, N_{a}, such that they represent the probability p(N_{a}) that a particle belongs to an aggregate of size N_{a} particles. In plots (a) to (c), the distributions are plotted vs. the logarithm of N_{a}. All distributions were generated from 10 final configurations at t_{run} = 1250τ generated in independent replicate simulations. ε′ values are in units of k_{B}Tσ^{2} and ε values are in units of k_{B}T. Aggregate sizes were computed using a linked-list method.^{48} (a) p(N_{a}) for the density-dependent system with ε′ =40. (b) p(N_{a}) for density-independent systems with ε values in the range ε = 2.84 → 42.68 that is explored in a density-dependent simulation with . Note that the apparent peaks located in the range 1 < log(N_{a}) < 2.0 (green and cyan curves) correspond to transient aggregates, i.e., aggregate “seeds” which dissolve soon after formation. For comparison, the red curve shows p(N_{a}) for the density-dependent system with ε′ = 40. (c) A weighted linear superposition (blue curve) of the p(N_{a}) curves computed in the density-independent simulations of (b) (see the ESI† for details). Again, the red curve shows p(N_{a}) for the density-dependent system with ε′ = 40. (d) Aggregate size distribution, including only aggregates of size N_{a} > 10 particles, for ε = 28.44 (density-independent) and ε′ = 40 (density-dependent). |
To determine whether this behaviour is intrinsic to the density-dependent interaction potential, we also performed a series of simulations with a standard, density-independent, Lennard-Jones interaction potential. We used values of the interaction strength ε in the range 2.84 → 42.68 k_{B}T: this corresponds to the ε-range that is sampled (viaeqn (16)) in the final configuration of our density-dependent simulation with ε′ = 40 (see the ESI† for details). A snapshot from the density-independent simulation with ε = 28.4 k_{B}T is shown in Fig. 10(d). Consistent with our observations from the field simulations, the aggregates are more varied in size and shape than those of the density-dependent system, and the non-condensed phase contains small clusters rather than single particles (compare to Fig. 10(b)). Phase separation in this system proceeds in a manner typical of spinodal decomposition, in that clusters form rapidly throughout the entire system (see Movie S9 in ESI†). Fig. 11(b) shows the aggregate size distribution (at time t_{run}) for density-independent simulations with ε = 2.84 → 42.68 (see also the ESI†). In contrast to the density-dependent case, none of these distributions are bimodal. Rather, the density-independent system forms either large aggregates (for ε > 25.6), or predominantly small aggregates (for ε ≤ 25.6). We also checked that a weighted linear superposition of the aggregate size distributions for different ε values cannot reproduce the aggregate size distribution produced by the density-dependent simulations (Fig. 11(c)) (see the ESI† for details). Thus, as we observed in our field simulations, the bimodal aggregate size distribution appears to be an intrinsic feature of the density-dependent system. For the particle-based simulations, we also observe that the peaks in the aggregate size distribution are broader for the density-independent simulations than in the density-dependent case (Fig. 11(b); see also Fig. 11(d) which confirms that the large-aggregate peak is narrower in the density-dependent system). Thus, density-dependent aggregation appears to confer increased “control” of aggregate size.
Given that our particle-based simulations exhibit qualitatively similar behaviour to that of our field simulations, we can conclude that “spinodal-nucleation”-like cluster formation, consisting of rather homogeneous aggregates that appear to show a barrier to growth, bimodal aggregate size distribution and a narrow cluster size distribution, are characteristic of systems with density-dependent attractive interactions, whether manifested as a cubic term in the Landau–Ginzburg free energy, or as a density-dependent potential in a particle-based model.
Combining a field theory approach with particle-based simulations, we observe behaviour that appears to be characteristic of systems with density-dependent attractive interactions. Specifically, we see an apparent barrier to cluster growth even in the spinodal regime, that we term “spinodal nucleation”, and bimodal aggregate size distributions with narrow peaks – corresponding to aggregates of rather uniform size in a “sea” of single particles. These characteristic behaviours can be understood as originating from a local free-energy barrier to cluster formation at low particle density, even in the spinodal region of the phase diagram. This barrier arises because local density fluctuations drive the system into the low-density coexistence region of the phase diagram, which in turn exists as a consequence of the cubic term that emerges from the Landau–Ginzburg free energy when the secretion of cohesion-inducing agents are considered, like that observed in nematic liquid crystals.^{39}
From the point of view of real biological cells, our model is highly simplified: we neglect, for example, cell motility, proliferation and death. Moreover, our assumption that the concentration field of cohesive agent is at steady state around a given particle is idealistic: in reality cohesive agents such as extracellular polymers are likely to degrade only slowly and may well accumulate in the system. Even for rapidly degrading agents the concentration field might be perturbed, for example by local fluid motion. Taking account of these factors would of course require a different modelling approach. Nevertheless, our work suggests that production of cohesive molecules might provide a route to control of aggregate size, and reveals interesting physics associated with the kind of density-dependent interactions that could be generated by the non-equilibrium process of cohesive agent production.
Footnote |
† Electronic supplementary information (ESI) available: For movies, and details on simulation methods, systems at different area fraction, structure factors, and phase diagrams for continuum model. See DOI: 10.1039/c9sm01462d |
This journal is © The Royal Society of Chemistry 2019 |