On the critical Casimir interaction between anisotropic inclusions on a membrane

Abstract

Using a lattice model and a versatile thermodynamic integration scheme, we study the critical Casimir interactions between inclusions embedded in a two-dimensional critical binary mixtures. For single-domain inclusions we demonstrate that the interactions are very long range, and their magnitudes strongly depend on the affinity of the inclusions with the species in the binary mixtures, ranging from repulsive when two inclusions have opposing affinities to attractive when they have the same affinities. When one of the inclusions has no preference for either of the species, we find negligible critical Casimir interactions. For multiple-domain inclusions, mimicking the observations that membrane proteins often have several domains with varying affinities to the surrounding lipid species, the presence of domains with opposing affinities does not cancel the interactions altogether. Instead we can observe both attractive and repulsive interactions depending on their relative orientations. With increasing number of domains per inclusion, the range and magnitude of the effective interactions decrease in a similar fashion to those of electrostatic multipoles. Finally, clusters formed by multiple-domain inclusions can result in an effective affinity patterning due to the anisotropic character of the Casimir interactions between the building blocks.

---

1 Introduction

Originally, the cell membrane was considered just as a physical barrier that kept cell components together. However, in the last decades, advances in experimental techniques, including atomic force and fluorescence microscopies, have enabled the probing of membranes' inner structure and composition to such extents that it has been necessary to rethink our understanding of its physics, chemistry and its role in biology.¹ Far from being a simple barrier, cell membranes are in fact very complex mixtures in which lipids and proteins meet, interact and self-assemble. A number of studies have shown that lipids are organized in the lateral dimension² and that most membrane proteins organize in clusters.^3–6 Further, there is now a growing consensus that such lateral organization and compartmentalisation play important roles in biological processes, such as in cell signalling and membrane trafficking;^7–9 and it has also been suggested to be involved in a number of diseases from HIV¹⁰ to liver¹¹ and prion diseases like Alzheimer.¹²

Biologically-motivated scenarios, such as the fence and pickets model¹³ or the lipid raft hypothesis,¹⁴ have been put forward to explain lateral membrane organizations. However, to date it remains unclear whether such scenarios are feasible from a chemical physics standpoint and, if possible, to what extent they are harnessed by biological processes in vivo. In this context, our aim here is to contribute towards understanding how inclusions (e.g. membrane proteins) may assemble into clusters within a cell membrane. Assuming the aggregation process occurs in thermodynamic equilibrium conditions, two possibilities come to mind. Firstly, it may occur by the free diffusion, collision and irreversible chemical bonding of the inclusions in a process reminiscent of a Diffusion Limited Aggregation.^15–17 Secondly, the inclusions may instead self-organise into metastable clusters, with a given equilibrium size distribution, in order to minimise the overall free energy. In this work we will concentrate on the latter.

There are various candidate mechanisms for interactions between somewhat large protein-like inclusions surrounded by smaller lipid species. They range from direct, specific interactions between the inclusions, such as via van der Waals or electrostatic forces,^18,19 to indirect interactions mediated by the lipid membranes. For the latter, one possible mechanism arises from minimising membrane deformations. Both theory and experiments have demonstrated that, by inducing curvature, proteins and colloidal particles can self-assemble into complex structures such as lines, rings and lattices.^20–23 More recently, ideas for a new class of indirect mechanisms have emerged, namely Casimir-like forces,^24–32 which could arise both from shape and composition fluctuations of the membrane. Here we will focus on composition fluctuations. Initially envisioned for bulk critical mixtures in three dimensions,³³ Casimir-like forces are expected to arise between any inclusions embedded in a fluctuating fluid when separated by a distance shorter than the correlation length of the fluid; and for fluid mixtures near a critical point, this correlation length diverges as we approach the critical temperature. Interestingly, close to room temperature, lipid membranes with critical composition pass through a miscibility critical point whose nature closely follows that of the two-dimensional Ising universality class.³⁴ This behaviour has been observed in giant plasma membrane vesicles isolated directly from living cells,³⁵ as well as in synthetic lipid mixtures.³⁶ It thus becomes plausible that Critical Casimir (CC) interactions between membrane proteins embedded in critical lipid mixtures could pave ways for their aggregation.

Machta et al.²⁴ recently demonstrated that two “like” inclusions immersed in a 2-dimensional critical binary lipid mixture have an attractive interaction, while “unlike” ones have a repulsive interaction. In their model each protein inclusion has a set uniform affinity with one of the lipid components. These findings are in line with those previously reported both theoretically and experimentally for CC interactions in 3 dimensions.^37–41 Recent extensive simulations of colloids immersed in a critical solvent further suggest strong similarities between the 2 and 3-dimensional cases whereby the observed variations in the phase behaviour appear to be accountable to differences in critical exponent and dimension of the solvent-inclusion interface.^42,43

Besides uniform inclusions, inhomogeneities and anisotropies at the fluid-inclusion interface have been theoretically shown to affect substantially CC interactions in 3 dimensions^40,44–46 to the extent that patterning can now be tuned and exploited for targeted colloidal self-assembly.^46,47 This observation is particularly interesting and relevant for CC forces between membrane proteins as proteins often have multiple domains with varying affinities with respect to the lipid species.^48–50 Thus, if CC interactions are to play a role in membrane protein aggregations, it is key to understand how their sign and magnitude can be tailored by heterogeneities in the protein domains and different affinities between the inclusions and the lipid species.

To address this issue, in this work we employ Monte Carlo (MC) simulations in a square lattice to simulate a binary mixture of lipids with two inclusions embedded in it. The free energy corresponding to the CC interactions between these inclusions is explicitly computed by using a thermodynamic integration scheme. The flexibility of our scheme allows us to (i) explore the effect of gradually modifying the inclusion boundary conditions, and (ii) study how the presence of different domains (anisotropy) can affect the interactions.

2 Methods

The Critical Casimir (CC) interaction is characterised by the free energy F_CC(d) of a critical fluid mixture within which lie two inclusions separated by a fixed distance d.³³ There are two common approaches to compute this CC free energy. Firstly, the free energy can be defined as minus the work done by the ensemble average of the mechanical force exerted by the fluid on the inclusions when brought from infinity to a finite separation d.^24,51 The equilibrium fluctuations of this mechanical force also inform experiments relying on mechanical probes to evidence CC interactions; however, some care must be taken when interpreting these fluctuations, depending on how the inclusions interact with the fluid mixtures.^26,52

Secondly, from an equilibrium statistical thermodynamics viewpoint, the statistics of configurations will be governed by a Boltzmann weight that depends on F_CC(d). Thus, a fruitful strategy is to compute the probability weights for inclusions at varying distances, which are readily available from simulations, and invert them to recover the CC free energy as a function of distance.⁴²

However, both the force and relative probability approaches usually need to be supplemented by a theoretical estimate at a reference distance because 2D CC interactions are very long range and it can be impractical to extend the calculations to get an effectively zero interaction at large distances.^24,42 In this work, we devise a non-distance-based thermodynamic integration scheme precisely to avoid the need for a reference point, which can be difficult to obtain for complex inclusion geometries. We will detail our approach in the following subsections.

2.1 The model

Our system consists of a critical lipid binary mixture in which two proteins are embedded. We use a lattice gas model of this system with a square lattice of N = 150 × 150 cells. Each cell is occupied either by a species of the binary mixture or part of a protein. This lattice size is chosen as it gives a good trade-off between reducing computational costs and finite size effects. By this we mean that up to N = 500 × 500, larger lattice sizes yield results whose difference with our N = 150 × 150 lattice size are not statistically significant.

We simulate the lipid mixture by making use of a two-dimensional Ising model, notorious for its critical behaviour, in which the spin variables can take values s = ±1 corresponding to lipid species A and B. The protein inclusions are modelled as two blocky patches occupying n_p = 12 cells each, whose size is determined by their effective radius, r, as shown in Fig. 1. For the inclusions, the average value of their spin variable is set to a target value s_t which can be different for the two proteins. To avoid a large standard deviation from this mean, spins belonging to a protein region are modelled with Potts-like spin variables s = k/2, with k an integer such that k ∈ [−2,2].

Fig. 1 Representation of the various stages of the model. Top-left figure: state α is a binary mixture at the critical point T_c. Top-right figure: state β with protein 1, alone, embedded within a critical mixture. Bottom-left figure: state γ with protein 2, alone, embedded within a critical mixture. Bottom-right: state δ with proteins 1 and 2 embedded in a critical mixture. The proteins' effective radius, r, is determined by the distance from the center of the patch to the furthest vertex.

To get the free energy of the system with two protein inclusions we make use of the 4 states represented in Fig. 1. The distance between proteins, d, is measured as the distance between their centers. These states are the following:

• State α comprising only the binary mixture at criticality (top left of Fig. 1) with hamiltonian:

• State β comprising the critical mixture and protein 1 (top right of Fig. 1) with hamiltonian:

• State γ comprising the critical mixture and protein 2 (bottom left of Fig. 1) with hamiltonian:

• State δ comprising the two proteins 1 and 2 separated by a distance d/r within the critical mixture (bottom right of Fig. 1) and with hamiltonian:

The first term in eqn (1)–(4) is equivalent to a 2D Ising model, while the second and third terms are quadratic terms whose role is to impose the average spin values inside the regions where the proteins are. Here J > 0 is the coupling parameter which characterises the energy cost for having two different species next to each other, stands for a sum over the 4 closest neighbours of cell i, h is a parameter setting the strength of the external potential imposing the average spin values inside the proteins, and s_t1 and s_t2 are the spin target values for proteins 1 and 2, respectively.

2.2 Thermodynamic Integration

Thermodynamic integration is a versatile method to compute free energy differences between two thermodynamic states. Let us consider an initial state with hamiltonian H_ini, a final state with hamiltonian H_fin, and a parameter-dependent hamiltonian H(λ) such that H(λ = 0) = H_ini and H(λ = 1) = H_fin. It can be shown that the free energy difference between the final and initial states can be computed by:⁵¹The equilibrium averages can readily be obtained from standard MC methods. It is important to note that, as a thermodynamic quantity, the free energy difference is only a function of the final and initial states, but not the path we have taken between them. Thus, for simplicity, we have taken a linear interpolation with crossover Hamiltonian

H(λ) = H_ini + λ(H_fin − H_ini). (6)

This then yields

Specific to the problem at hand, here we have four thermodynamic states α, β, γ and δ (see Fig. 1), and we denote their free energies as F_α, F_β, F_γ and F_δ(d/r). We finally denote F_CC(d/r) as the CC free energy and define it as the work done to bring the two proteins to within a distance d/r from infinity. It follows from this definition that:

F_CC(d/r) = F_δ(d/r) − (F_β + F_γ) + F_α, (8)

where (F_β + F_γ) − F_α represents the free energy of the system with the two proteins being infinitely far apart in the critical mixture. From eqn (7), F_CC(d/r) can then be expressed as a function of thermodynamic integrals only:which can be readily estimated from our MC simulations. We note that, for a given set of protein inclusions, the second thermodynamic integral only needs to be computed once, while the first integral has to be repeated for various values of d/r. The advantage of eqn (9), which resemble to some extent the parameter variation method used in ref. 44, is that it avoids the problem of having to supply a theoretical estimate of the reference free energy, which tends to be very much system-dependent and difficult to compute analytically.

2.3 Simulation details

We perform standard Metropolis Monte Carlo simulations at the critical temperature of the Ising model with non-conserved order parameter. In our simulations, it corresponds to J/k_BT = 2.27, in agreement with analytical predictions by Onsager and previous simulation results.^53,54 For each system we employ three random different initial configurations, which are equilibrated for 5 × 10⁵ cycles, and sampled for 5 × 10⁶ cycles. Each cycle consists of N random single spin flips, which ensures that all particles in the system can be selected in each cycle, and configurations are saved every 1000 cycles. In order to determine the error in our calculations we use block analysis⁵⁵ to determine the number of independent measurements of the free energy. The error is obtained from the standard deviation with a confidence interval of 95%.

When performing the thermodynamic integration, there are two technical aspects worth commenting. The first aspect is the choice for the parameter h in eqn (1)–(4). On the one hand, its value must be large enough to effectively pin the spin values inside an inclusion such that the associated energy variation 〈H_fin − H_ini〉_λ → 0 as λ → 1. On the other hand, the accuracy of our estimate of the CC free energy, F_CC, is better when we have more points substantially contributing to the numerical integration in eqn (7). As a consequence, h must be in a window of values guaranteeing both fidelity to the model we want to implement and reliability of the integration method.

We further illustrate this issue in Fig. 2, where we show a typical result of our thermodynamic integration method for three values of h = 50k_BT, 250k_BT and 500k_BT. In this example, we compute the CC free energy difference between states γ and δ. For the largest value of h, the contribution to the integral in eqn (7) primarily comes from a small window of λ close to λ → 0 and only involves a couple of integration points, leading to a poor estimate of the integral. From this perspective, smaller values of h are better. However, for small values of h, we notice that the integrand in eqn (7) does not converge to zero as λ → 1 (see the inset of Fig. 2). Physically this means the spin value of protein inclusions is not correctly set to the target value s_t. In this work, we find that we obtain equivalent results for the Critical Casimir free energy if we use h = 100–250k_BT. For the rest of this paper, we use h = 250k_BT.

Fig. 2 Average values of the integrand in eqn (7) as a function of λ for different values of the external field: h = 50k_BT (red circles), h = 250k_BT (green squares) and h = 500k_BT (blue diamonds). The initial state corresponds to state γ with s_t = −1, and the final state corresponds to state δ with s_t = −1 for the first inclusion and s_t = 0 for the second inclusion. The distance between the two inclusions is d/r = 2.24. Inset: Magnification in the region of high values of λ. The lines for h = 250k_BT (green squares) and h = 500k_BT (blue diamonds) are indistinguishable.

The second aspect concerns the difference in degrees of freedom when a given lattice point represents a lipid species or part of a protein inclusion. For instance, when computing ΔF_α/β, the protein sites in state β have five possible spin values (Potts model with s = k/2, k ∈ [−2,2]), whereas the equivalent sites in state α only have two possible spin values (Ising model with s = ±1), because in state α they represent lipid species. Since we are computing the CC free energy at varying distances between the protein inclusions, in principle we must account for corrections due to variations in degrees of freedom explicitly in the computations for F_CC(d/r). To do this, we carry out two sets of thermodynamic integration calculations. In the first set, see sketch in Fig. 3(b), the initial state corresponds to the case where there is no interaction neither between the cells where the inclusions are located nor with their surrounding cells. The final state is where all cells are interacting and they all have two possible spin values (Ising model). We denote this free energy difference as ΔF_Ising(d/r), where d once again is the separation between the two inclusions and r is the radius of the inclusions. For the second set of calculations, see sketch in Fig. 3(b), the initial state is the same as before. However, for the final state, the cells where the inclusions are located can now have five possible spin values (Potts model). The surrounding sites still have s = ±1. The free energy difference for this case is denoted as ΔF_Potts(d/r). The difference ΔF_Ising(d/r) − ΔF_Potts(d/r) is therefore the correction in the CC free energy due to changes in the spin degrees of freedom as a function of the normalised separation d/r, and the typical results are shown in Fig. 3(a) for a 150 × 150 square lattice. In practise, these corrections are small and, in fact, within the uncertainty of our typical thermodynamic integration results.

Fig. 3 (a) ΔF_Ising(d/r) − ΔF_Potts(d/r) is the correction in the CC free energy due to variations in the spin degrees of freedom between lattice points representing lipid species and the protein inclusions, as a function of the normalised separation d/r. The corrections are negligible. (b) Sketch of the systems simulated. In the middle, the inclusion cells (inside the circles) do not interact with their surrounding cells or with each other. On the top right, all lattice points, including the inclusion cells, interact via standard Ising model. On the bottom right, the inclusion cells are modelled with Potts-like spin variables, while the other lattice sites have Ising spin variables.

3 Results and discussion

We study the effective interaction between two protein inclusions embedded in a lipid membrane depending on the binding affinity that the proteins have for the lipid species A and B present in the membrane. For instance, a protein (or domain in a protein) labelled as −1 strongly prefers being surrounded by lipid species A, while a protein (or domain in a protein) labelled as 1 strongly prefers lipid species B. A value of 0 means that it has no preference for either of the species.

3.1 Single-domain proteins

We first illustrate the approach to the critical temperature T_c by plotting two-dimensional maps of the ensemble average spin value denoted φ for various temperatures in the neighbourhood of T_c. The presence of uniform inclusions with set spin value s_t = −1 introduces inhomogeneities into the spin system which persist over a distance of the order of the correlation length. Plotting the φ-maps enables us to have a direct visualisation of the correlation length. Indeed, the divergence of the correlation length at T_c gives rise to long-range CC forces to exist.

As shown in Fig. 4, approaching T_c from above, the distance over which the inclusions influence the rest of the spins in the system is initially rather short (panel (a)), but increases more and more (panels (b) and (c)) until the correlation length covers the whole box at T_c.

Fig. 4 Visualising the approach to T_c. 2D maps of the average spin value φ in the vicinity of 2 inclusions (disks) with target spin value set at −1 (corresponding to a dark colour in the map). (a) T = 1.33T_c, (b) T = 1.06T_c and (c) T = T_c.

We next test our thermodynamic integration scheme by computing the CC free energy for proteins that consist of one single domain. We explore the variations of the protein–protein interactions as a function of their separation for different protein binding affinities to the lipid species. In our approach, we gradually change the binding affinity of one of the proteins to the lipid species, while keeping the other one constant. More specifically, we consider systems in which protein 1 has a spin target value of −1, while protein 2 can get spin target values of 1, 0.5, 0, −0.5 or −1 (see Fig. 5 for a description). The results for these systems are summarised in Fig. 5.

Fig. 5 CC interactions for single-domain inclusions. (a) From top to bottom: 2D maps of the average spin value φ for cases (−1,1), (−1,0.5), (−1,0), (−1,−0.5) and (−1,−1). (b) Effective interactions as a function of distance for varying affinities of the 2nd inclusion: system (−1,1) (orange diamonds), system (−1,0.5) (blue triangles), system (−1,0) (green stars), system (−1,−0.5) (red squares) and system (−1,−1) (black circles). (c) Fitting parameters to F_CC(d/r)/k_BT = ζ(d/r)^−ν for the systems shown in (b).

We start by looking at the limit in which both proteins have opposite binding affinities for the two lipid species, system (−1,1). Here, protein 1 has strong preference for lipid A and protein 2 has strong preference for lipid B, which we call the “unlike” limit. In this limit the effective interaction is clearly repulsive (orange diamonds in Fig. 5) and can be quite strong with a free energy barrier rising up to 3.5k_BT. Then, as the binding affinity of protein 2 is gradually varied, and its preference becomes weaker for lipid species B and stronger for A, the strength of the repulsive interaction decreases (blue triangles in Fig. 5). At the point in which protein 2 has no preference for either of the lipid species, system (−1,0), we find a crossover between the two regimes, and CC interactions become negligible with respect to k_BT (green stars in Fig. 5). Finally, as protein 2 keeps on increasing its affinity for lipid A we get into the attractive regime (red squares in Fig. 5). In the “like” limit in which both proteins have the same strong affinity to one of the lipid species, system (−1,−1), we find that the magnitude of attraction is the strongest with a well depth of about 1k_BT (black circles in Fig. 5). However, this attraction is much weaker than the magnitude of the repulsion in the “unlike” limit, by more than a 3 fold difference. That like proteins attract supports the idea that proteins of the same kind tend to coalesce in critical lipid mixtures in order to minimise the total free energy of the system.

We note that the free energy curves in Fig. 5 appear to have an offset with respect to the zero free energy baseline. However, this apparent offset is just a consequence of the very slow decay of these curves to zero as d/r → ∞ and thus of the very long range nature of CC interactions in two dimensions. We further find that the free energy dependence with distance can be well fitted to a power law F_CC(d/r)/k_BT = ζ(d/r)^−ν, shown as plain lines in Fig. 5. The values for ζ and ν are tabulated in panel (c) in Fig. 5. Smaller values of ν signify longer-range of interactions.

Our results for the like (black circles) and unlike (orange diamonds) limits in Fig. 5 corroborate the findings of Machta et al.²⁴ where we find excellent agreement without any fitting parameter. In their simulation work using Bennett algorithm,^56,57 they take advantage of asymptotic results from conformal field theory for the reference free energy at large distances. Our thermodynamic integration approach here circumvents the need for supplying the reference free energy at large distances, and thus it allows more complex scenarios to be computed readily, as we shall show in the next section for inclusions with multiple domains. Our results are also in good qualitative agreement with the experimental findings of ref. 41, where the CC forces between colloids and a wall were tuned by varying the preferential absorption properties of the wall for the critical mixture's components.

Focusing on the attractive regime which is the strongest for like inclusions, we next look at the magnitude of the CC interactions between like proteins, as we vary their binding affinity with the surrounding matrix from strong preference for lipid A (system (−1,−1)) to no preference for either lipid (system (0,0)). Note that we do not investigate explicitly the case where the two proteins prefer lipid B since it is equivalent to the system (−1,−1). Our results are summarised in Fig. 6. It is found that the greater the affinity of the two proteins for one of the lipid species, the greater the magnitude of the resulting CC interactions. In fact, for system (0,0) where the inclusions have no preference for either of the lipid species, we cannot measure any net interaction within the limits of our computational accuracy. For system (−0.5,−0.5) we clearly observe an attractive interaction (negative free energy). However, the CC force (the slope of the free energy with distance) is very small. It is akin to the case where the system sits in an effective negative square well potential.

Fig. 6 CC interactions for single-domain inclusions. (a) From top to bottom systems: (0,0), (−0.5,−0.5) and (−1,−1). (b) Effective interactions as a function of distance for the different systems in (a): system (0,0) (green diamonds), system (−0.5,−0.5) (red squares) and system (−1,−1) (black circles).

3.2 Multiple-domain proteins

In the previous section we noted the absence of CC interaction at any distance when one or both of the proteins do not have any preference for either of the lipid species. In this section we investigate whether this still remains the case for identical proteins with multiple domains each of which being given a target affinity of either −1 or 1—i.e. with strong binding affinity for either lipids A or B—but such that the overall net affinity of each protein is zero. To some extent, this approach with multiple-domain proteins represents a more realistic model of proteins which usually comprise several domains with different binding affinities to the surrounding lipids,^48–50 thus making the present study more relevant to biological systems.

We start by studying interactions between two proteins with two different domains each. These proteins can have different relative orientations depending on which domains of the proteins are facing each other (see Fig. 7(a) for a description). Our aim is to study how orientation can affect these interactions and if there are preferential orientations that could lead to patterning in protein clusters. Given that these proteins do not have a net preference for any lipid type, one could naively expect no interaction between them as in the cases of the (0,0) and the (−1,0) pairs of single-domain proteins. However, our results as shown in Fig. 7(b) show again the presence of two different regimes and a crossover between them. We find that systems 1 and 2 appear in the repulsive regime while systems 4 and 5 fall in the attractive one. System 3, corresponding to a 90° relative orientation, is slightly attractive, though its magnitude is very small compared to the other systems.

Fig. 7 CC interactions for two-domain inclusions. (a) 2D maps of the average spin field φ for different relative orientations of the affinity dipoles. From top to bottom: system 1, system 2, system 3, system 4 and system 5. The black domains are set with s_t = −1, while the white domains have s_t = 1. (b) Effective interactions as a function of distance for the different systems in (a): system 1 (black circles), system 2 (orange diamonds), system 3 (green stars), system 4 (red squares), and system 5 (blue triangles). (c) Fitting parameters to F_CC(d/r)/k_BT = ζ(d/r)^−ν for the systems shown in (b).

A closer look at systems 1 and 2 shows that, in both cases, proteins are facing each other with domains that have opposite preferences for the lipid species, giving rise to predominantly “unlike” interactions. Therefore, this situation is qualitatively analogous to the system (−1,1) for single-domain proteins. On the other hand, systems 4 and 5 are confronting domains of the same kind, leading to predominantly “like” interactions as in system (−1,−1) for single-domain proteins. Lastly, in system 3, proteins are orientated in such a way that half of protein 1 has “like” interactions with protein 2, while the other half has “unlike” interactions. This case is a good example of the non-additive nature of CC interactions. If they were simply additive, taking advantage of the results in Fig. 5, we would have expected the net interaction to be repulsive. The magnitude of repulsive interactions between unlike uniform inclusions is stronger than attractions between like ones. However, as shown in Fig. 7, for system 3, we observe a very weak attractive interaction instead.

When compared with single-domain proteins, the range of the emerging CC interaction is clearly shorter in all cases, getting close to zero at distances of about ten times the radius of the protein; we interpret this as being a partial screening effect analogous to electrostatic screening in dielectric systems and owing to proteins comprising domains with opposing affinities. Panel (c) in Fig. 7 tabulates the power law fit of the CC free energy with distance, F_CC(d/r)/k_BT = ζ(d/r)^−ν. Here, the exponents ν are clearly larger for the values for single-domain inclusions in Fig. 5.

In order to further investigate the effect of the domains we next turn to proteins with four distinct domains. In this case only two possible orientations are of interest, see Fig. 8(a). We present our results in Fig. 8(b). Again, in spite of having no net preference for either of the species, proteins show either attractive or repulsive interactions. A similar analysis to that done with two-domain proteins show that system 6 confronts domains with opposite preferences for the lipid species, falling, therefore, in the repulsive regime. On the other hand system 7 has domains with the same preference facing each other, leading to attractive interactions. It is worth pointing out that, in these cases, both the range and the magnitude of the interactions are smaller than those observed for single- or two-domain inclusions. Fig. 8(c) tabulates the power law fit for the CC free energy with distance between two four-domain inclusions. Eventually in the limit of the domain size going to zero, we expect to recover the uniform case of s_t = 0. Comparing the results in Fig. 5, 7 and 8, we also observe that the magnitude of the repulsive interactions falls faster than attraction interactions with increasing number of domains.

Fig. 8 CC interactions between four-domain inclusions. (a) 2D maps of the average spin value φ for two different orientations of the affinity quadrupoles. From top to bottom: system 6, system 7. The black domains are set with s_t = −1, while the white domains have s_t = 1. (b) Effective interactions as a function of distance for the different systems in (a): system 6 (black circles), system 7 (red squares). (c) Fitting parameters to F_CC(d/r)/k_BT = ζ(d/r)^−ν for the systems shown in (b).

The observed behaviours for multiple-domain inclusions are in good qualitative agreement with the findings of Toldin et al. for CC forces between chemically striped surfaces immersed in critical liquid mixtures.⁴⁴ In their study, they confront a surface with a strong preference for one component of the binary liquid mixtures with a striped surface in which the preference of the stripes for the liquid components is alternating. They observe a crossover between a regime with noticeable CC interaction and one with no noticeable CC interaction (similar to our (0,0) and (−1,0) cases for single-domain proteins) as the strip width is decreased from large values to zero.

Finally, we note that the range and strength of the CC interactions depend only on the lattice points which are at the boundary of the inclusions, and not the inner ones. This is as expected because we only have nearest neighbour interactions, and thus the inner sites do not interact directly with the surrounding cells representing the lipid mixtures.

4 Conclusions and future work

In this work, we have numerically computed the CC free energy between inclusions in two-dimensional critical binary mixtures. To this end, we have introduced a versatile thermodynamic integration scheme whose calculations for single-domain inclusions compare very favourably with the existing literature on the subject.²⁴ This has enabled us to identify repulsive and attractive CC regimes respectively for unlike inclusions—with opposing affinities for the surrounding species comprising the binary mixture—or like inclusions—with the same affinity for the surrounding species. We noted that the CC interactions for single-domain inclusions are very long range and their magnitude is increasing with the strength of the affinity of the domain with either of the surrounding species.

A key advantage of our thermodynamic integration approach is that we do not need to supply a reference free energy at large distances, which can be notoriously difficult to compute analytically for complex scenarios. In turn this allows us to study inclusions comprising multiple domains such that, overall, they do not have any net preference for either of the surrounding species in the mixture.

The presence of domains with opposing affinities does not cancel the CC interactions altogether but instead leads to a decrease in the range and magnitude of the effective interactions in a manner reminiscent of electrostatic multipoles. For such multiple-domain inclusions, both repulsive and attractive CC interactions were observed depending on their relative orientations.

If biologically relevant lipid membranes are close to criticality as suggested by some studies,^34,35 then these findings may prove useful to better understand lateral spatial organisation within cellular membranes. As a matter of fact, attractive CC interactions between the single-domain inclusions of the same kind provide both sufficient magnitude and range for protein clustering to occur in equilibrium conditions.^42,47 Utilising multiple-domain inclusions as more realistic models of membrane proteins lead us to the conclusion that the emerging CC interactions between them also suffice for them to aggregate, even where there is no overall net affinity to the surrounding species in the mixture.

There are several avenues for future work. Firstly, it would be interesting to study the clustering, and more generally the phase diagram, of the multiple-domain inclusions. For instance, in Fig. 9(a) we show the relative stability of the different orientations for two-domain proteins. Since some orientations are energetically favourable, we would expect proteins to rotate in order to minimise the overall free energy of the system. Our preliminary results suggest that, as these inclusions start to cluster, they will form an alternating patterning in affinities with respect to the surrounding species (Fig. 9(b)). In this context, it would be exciting to explore how CC interactions could be tailored to allow components in biological membranes achieve specificity as suggested in ref. 58. Such strategies could also be exploited for the self-assembly of anisotropic (e.g. Janus) colloidal particles.^46,59 Secondly, Casimir interaction due to shape fluctuations of the membrane has also been proposed as a non-specific mechanism for membrane proteins to interact.^28–32 Here it will be instructive to compare the strengths of the two possible Casimir interactions, and map regions in parameter space (temperature, composition and bending rigidity) where one mechanism dominates over the other, and where they interfere either constructively or destructively. For instance, Dean et al.⁶⁰ suggests that taking into account membrane shape fluctuations can result in a shift in the critical temperature at which phase separation occurs for the lipid mixtures. Stress tensor-based methods following ref. 61 could help inform on these competing effects.

Fig. 9 Stability, aggregation and patterning for two-domain proteins.

Conflicts of interest

There are no conflicts of interest to declare.

Acknowledgements

We thank Pietro Cicuta, Mark Miller and Lorenzo di Michele for helpful discussions. We acknowledge funding from EPSRC, grant EP/J017566/1. Research data supporting this article can be accessed at https://doi.org/10.15128/r26d56zw62r.

References

1. Y. Shan and H. Wang, Chem. Soc. Rev., 2015, 44, 3617–3638.
2. A. Kusumi and Y. Sako, Curr. Opin. Cell Biol., 1996, 8, 566–574.
3. S. Uhles, T. Moede, B. Leibiger, P.-O. Berggren and I. B. Leibiger, J. Cell Biol., 2003, 163, 1327–1337.
4. M. Kai, F. Sakane, Y.-J. Jia, S.-I. Imai, S. Yasuda and H. Kanoh, J. Biochem., 2006, 140, 677.
5. S. H. Low, A. Vasanji, J. Nanduri, M. He, N. Sharma, M. Koo, J. Drazba and T. Weimbs, Mol. Biol. Cell, 2006, 17, 977–989.
6. J. J. Sieber, K. I. Willig, R. Heintzmann, S. W. Hell and T. Lang, Biophys. J., 2006, 90, 2843–2851.
7. M. A. Alonso and J. Millán, J. Cell Sci., 2001, 114, 3957–3965.
8. E. Ikonen, Nat. Rev. Mol. Cell Biol., 2008, 9, 125–138.
9. U. E. Martinez-Outschoorn, F. Sotgia and M. P. Lisanti, Nat. Rev. Cancer, 2015, 15, 225–237.
10. K. Simons and R. Ehehalt, J. Clin. Invest., 2002, 110, 597–603.
11. J. C. Holthuis and A. K. Menon, Nature, 2014, 510, 48–57.
12. P. Critchley, J. Kazlauskaite, R. Eason and T. J. Pinheiro, Biochem. Biophys. Res. Commun., 2004, 313, 559–567.
13. A. Kusumi, C. Nakada, K. Ritchie, K. Murase, K. Suzuki, H. Murakoshi, R. S. Kasai, J. Kondo and T. Fujiwara, Annu. Rev. Biophys. Biomol. Struct., 2005, 34, 351–378.
14. K. Simons and E. Ikonen, Nature, 1997, 569–572.
15. M. J. Saxton, Biophys. J., 1993, 61, 119–128.
16. M. J. Saxton, Biophys. J., 1993, 64, 1053–1062.
17. X. Tian, H. Zheng, P. T. Matsudaira and U. Mirsaidov, Microsc. Microanal., 2016, 22, 808.
18. B. H. Honig, W. L. Hubbell and R. F. Flewelling, Annu. Rev. Biophys. Biophys. Chem., 1986, 15, 163–193.
19. A. Lee, Biochim. Biophys. Acta, 2003, 1612, 1–40.
20. B. J. Reynwar, G. Illya, V. A. Harmandaris, M. M. Muller, K. Kremer and M. Deserno, Nature, 2007, 447, 461–464.
21. J. A. G. Briggs, J. D. Riches, B. Glass, V. Bartonova, G. Zanetti and H.-G. Kräusslich, Proc. Natl. Acad. Sci. U. S. A., 2009, 106, 11090–11095.
22. A. Saric and A. Cacciuto, Soft Matter, 2013, 9, 6677–6695.
23. A. Vahid and T. Idema, Phys. Rev. Lett., 2016, 117, 138102.
24. B. B. Machta, S. L. Veatch and J. P. Sethna, Phys. Rev. Lett., 2012, 109, 138101.
25. G. Bimonte, T. Emig and M. Kardar, Europhys. Lett., 2013, 104, 21001.
26. A.-F. Bitbol and J.-B. Fournier, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2011, 83, 061107.
27. L. M. Woods, D. Dalvit, A. Tkatchenko, P. Rodriguez-Lopez, A. W. Rodriguez and R. Podgornik, Rev. Mod. Phys., 2016, 88, 045003.
28. A.-F. Bitbol, P. G. Dommersnes and J.-B. Fournier, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2010, 81, 050903.
29. R. Golestanian, M. Goulian and M. Kardar, Europhys. Lett., 1996, 33, 241.
30. C. Yolcu, I. Z. Rothstein and M. Deserno, Europhys. Lett., 2011, 96, 20003.
31. H.-K. Lin, R. Zandi, U. Mohideen and L. P. Pryadko, Phys. Rev. Lett., 2011, 107, 228104.
32. T. R. Weikl, Europhys. Lett., 2001, 54, 547.
33. M. E. Fisher and P. G. de Gennes, C. R. Acad. Sci. Paris B, 1987, 287, 207–209.
34. A. R. Honerkamp-Smith, P. Cicuta, M. D. Collins, S. L. Veatch, M. den Nijs, M. Schick and S. L. Keller, Biophys. J., 2008, 95, 236–246.
35. S. L. Veatch, P. Cicuta, P. Sengupta, A. Honerkamp-Smith, D. Holowka and B. Baird, ACS Chem. Biol., 2008, 3, 287–293.
36. S. D. Connell, G. Heath, P. D. Olmsted and A. Kisil, Faraday Discuss., 2013, 161, 91–111.
37. M. Krech, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1997, 56, 1642.
38. C. Hertlein, L. Helden, A. Gambassi, S. Dietrich and C. Bechinger, Nature, 2008, 451, 172–175.
39. A. Gambassi, A. Maciołek, C. Hertlein, U. Nellen, L. Helden, C. Bechinger and S. Dietrich, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2009, 80, 061143.
40. M. Trondle, S. Kondrat, A. Gambassi, L. Harnau and S. Dietrich, Europhys. Lett., 2009, 88, 40004.
41. U. Nellen, L. Helden and C. Bechinger, Europhys. Lett., 2009, 88, 26001.
42. J. R. Edison, N. Tasios, S. Belli, R. Evans, R. van Roij and M. Dijkstra, Phys. Rev. Lett., 2015, 114, 038301.
43. N. Tasios and M. Dijkstra, J. Chem. Phys., 2017, 146, 134903.
44. F. Parisen Toldin, M. Tröndle and S. Dietrich, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2013, 88, 052110.
45. P. Nowakowski and M. Napiórkowski, J. Chem. Phys., 2014, 141, 064704.
46. M. Labbe-Laurent and S. Dietrich, Soft Matter, 2016, 12, 6621–6648.
47. V. D. Nguyen, M. T. Dang, T. A. Nguyen and P. Schall, J. Phys.: Condens. Matter, 2016, 28, 043001.
48. G. B. Cohen, R. Ren and D. Baltimore, Cell, 1995, 80, 237–248.
49. T. Pawson and P. Nash, Science, 2003, 300, 445–452.
50. B. J. Mayer, Nat. Rev. Mol. Cell Biol., 2015, 16, 691–698.
51. D. Frenkel and B. Smit, Understanding Molecular Simulation, Academic Press, 2002.
52. D. Bartolo, A. Ajdari and J.-B. Fournier, Phys. Rev. Lett., 2002, 89, 230601.
53. L. Onsager, Phys. Rev., 1944, 65, 117–149.
54. K. K. Mon and K. Binder, J. Chem. Phys., 1992, 96, 6989–6995.
55. H. Flyvbjerg and H. G. Petersen, J. Chem. Phys., 1989, 91, 461–466.
56. C. H. Bennett, J. Comput. Phys., 1976, 22, 245–268.
57. C. Jarzynski, Phys. Rev. Lett., 1997, 78, 2690–2693.
58. J. J. Sieber, K. I. Willig, R. Heintzmann, S. W. Hell and T. Lang, Biophys. J., 2006, 90, 2843–2851.
59. E. Noruzifar and M. Oettel, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2009, 79, 051401.
60. D. S. Dean, V. A. Parsegian and R. Podgornik, J. Phys.: Condens. Matter, 2015, 27, 214004.
61. A.-F. Bitbol, L. Peliti and J.-B. Fournier, Eur. Phys. J. E: Soft Matter Biol. Phys., 2011, 34, 53.

---