Proteins at curved fluid–fluid interfaces in a coarse-grained model

We employ an empirical coarse-grained model with a proposed Gaussian-like interfacial potential to describe proteins at curved fluid–fluid interfaces such as occurring in bubbles and droplets. We consider the air–water and oil–water interfaces. We study the mass distributions and the geometry of the aqueous proteins as a function of the radius of curvature for protein G and two lipid transfer proteins. At curved interfaces the distortion of the proteins is different than at flat interfaces. We find that the proteins come closer to the surface of a bubble than to the surface of similarly curved droplet. In addition, the bubbles adsorb more proteins. We identify the pinning residues. We demonstrate the existence of the second layer in the density profile for sufficiently dense solutions.


Introduction
Understanding of the behavior of proteins at fluid-fluid interfaces has many ramifications. It is important for the food industry, for instance, in the context of stabilization of the oil-water emulsions [1][2][3][4] and of the beer foam [5][6][7]. It is relevant in physiology as it impacts stabilization of lung alveoli against collapse [10], defence mechanisms of the lungs against inhaled pathogens [11], and stabilization of saliva [14]. It is important for explaining the viscous properties of biofilms [12,13]. The problem also connects to the liquid-liquid phase transitions leading to the formation of proteinaceous liquid droplets that may act as the membraneless organelles [15][16][17][18].
Such interfaces may become curved when they form boundaries of droplets or bubbles. Their radii of curvature may span a big range of values, starting from the nanoscale. Even the nominally flat interfaces undergo shape fluctuations that introduce nanoscale curvature. This is illustrated in figure 1 which shows snapshots of 'air'-water and oil-water interfaces derived by our all-atom simulations outlined in reference [19] (the figure was not shown there).We used the NAMD simulational package [20] at room temperature with the CHARM22 force field [21]. The water molecules were described by the TIP3P model 3 Author to whom any correspondence should be addressed. [22] and the oil molecules were represented by the chains of oleic acid of 18 carbon atoms each [23]. The panels on the right show that vertical shape fluctuations range between -6 and +6 Å for the air-water case and between -10 and +10 Å for the oil-water situation. The interface have been obtained by setting the bin size as 1 Å in both the x-and y-axis directions and 1.2 Å in the z-direction. They are seen to produce an instantaneous landscape of hills and valleys, both curved, in the interface. The oil droplets have properties similar to those of lipid vesicles and membranes [24,25] which are also curved. Thus it is interesting to ask what is the impact of the curvature on the proteins at interfaces.
In our previous studies [7,26,27], we have proposed an empirical structure-based coarse-grained molecular model to study globular proteins at flat fluid-fluid interfaces. The crucial feature of the model is that the interactions between beads representing residues are defined through effective native contacts. These contacts are determined from atomic-level considerations pertaining to the native conformation [8,9]. (See also a lattice version of the model [28] and another realization [29]). It should be noted that there are many types of coarsegrained models of proteins [30][31][32][33]. In particular, the physicsbased approach of Liwo et al [32] introduces anisotropic sub-units similar to those employed in the coarse-grained modelling of liquid crystals [34]. What distinguishes our model Left: snapshots of the air-water (top) and oil-water (bottom) interfaces at T = 300 K. The simulation box extends between −L x and L x in the x and y directions, and −L z and L z along the z-direction. We take L x = 50 Å and L z is equal to 100 Å for the air-water interface and 50 Å for the oil-water one. In both cases, the water molecules are placed initially in the space corresponding to z 0. The air/oil molecules are placed in the upper half of the box. Right: the distribution of the z coordinates of the surface atoms as obtained from 10 frames that are 5 ns apart. The number density profiles of water molecules at the air-water oil-water interfaces along z-axis and the radial direction in the x-y plane can be found in figure 1 of reference [19]. is its simplicity, resulting from the structure-based approach, the usage of spherical interaction centers, and, perhaps more importantly, the possibility to include the interfaces in an empirical way.
The interfaces are introduced through hydropathy-related forces that are maximal at the center of the interface and extend away from the center in a Gaussian fashion with a width, W, that depends on the system. We have also provided all-atom justification of the model [19]. All-atom density profiles indicate W to be of order 5 Å and 10 Å for the air-water and oil-water cases respectively. We have found that the hydropathy forces result in surface trapping and deformation of the proteins since the polar and charged groups seek bulk water whereas the hydrocarbon groups avoid water. The proteins remain close to their native states in bulk water and their interfacial behavior can be different at various interfaces as has been discussed in details in our previous work [19]. In addition, the proteinic films exhibit glassy behavior [26].
The deformation of the proteins may cause denaturation and, as pointed out by D'Imprima et al [35] may severely affect structural measurements in the cryo-EM studies that involve vitrified aqueous solutions. Their studies on yeast fatty acid synthase indicate that about 90% of complexes are partly or fully denatured because proteins tend to adsorb to the air-water menisci in unsupported films. D'Imprima et al demonstrate that this difficulty can be alleviated by plunge freezing on a substrate of hydrophilized graphene.
Menisci arising in the cryo-EM studies, like the droplets and bubbles, are also endowed with a curvature. Often, the size of the globular proteins is small compared to the corresponding radii of curvature and the approximation of the flat interface is adequate. However, for small droplets and/or large protein complexes the sizes may become comparable and the effects of the curvature need to be assessed. Our simplified empirical model offers an opportunity to perform such a pilot assessment.
Here, we consider air-water and oil-water spherical interfaces of radius R where water is either inside or outside the sphere. We generalize our model to deal with the curved situations, as described in the methods section. We generate ensembles of proteins and determine the protein number density profiles in the flat case. The density profiles display large peaks just at the interface. However, for a sufficiently large density, secondary peaks develop at distances that are further away from the interface.
We also consider single proteins and characterize the dependence of their deformation on R. For films of proteins, additional deformation may arise due to the interprotein interactions but the dominant effect is due to the hydropathy forces [26].
As an illustration, we focus on three proteins: protein G, with 56 residues and the PDB structure code 1GB1, and two variants of the barley lipid transfer protein LTP1 comprising 91 residues. One variant corresponds to the PDB structure 1LIP has all 4 possible disulfide bonds present and will be denoted by 1LIP (4). The second has its all disulfide bonds reduced and will be denoted by 1LIP(0). Both variants are present in beer and contribute to the longevity of the foam making bubbles in beer [7]. The native radius of gyration for 1GB1 and 1LIP(4) is 8.2 and 11.2 Å respectively.

Methods
In the coarse-grained model we use, the proteins are represented by the α-C atom of residues and the interactions are derived from the native structure. The version we implement is described in references [8,9,36]. The interactions of bonded residues (the peptide and disulfide bonds) are described by the harmonic potentials and the backbone stiffness is accounted for by the chirality potential. The non-bonded native interactions of proteins are represented by contacts. Their list is determined based on the existence of overlaps between effective spheres associated with the heavy atoms in the native state. The contact potentials are taken in the Lennard-Jones form with the depth of that is approximately equal to 1.6 kcal/mol as determined by benchmarking the stretching simulations against the experimental data for proteins in bulk water [8]. The pairs of residues that do not form native contacts interact through the repulsive part of the Lennard-Jones potential. The solvent in our model is implicit, and it is introduced through the random kicking forces. The temperature, T, is controlled by the overdamped Langevin thermostat. The resulting room temperature is close to 0.3-0.35 /k B , where k B is the Boltzmann constant. All of our simulations in this work are performed at The characteristic time scale, τ , is of order 1 ns.
The inter-protein interactions are governed by the parameter λ. When λ is zero, they are purely repulsive. Otherwise, we introduce the attractive Lennard-Jones interactions between pairs of the hydrophobic residues with the energy parameter λ . We consider λ = 0 and 1.
The flat fluid-fluid interfaces are introduced through the hydropathy force that acts normal to the interface-this is the z direction. Here, i is the sequential label of the residue and z i is the Cartesian coordinate of the residue along the z-axis. It is assumed that the center of the planar interface is at z = 0. The region with the bulk water corresponds to the negative z-axis. The region corresponding to the positive z-axis is either air or oil. q i is the hydropathy index of the residue: it is positive for hydrophobic residues and negative for hydrophilic ones. We take the values of q i from the work of Kyte and Doolittle [37], in which q i scales from −4.5 for argine to +4.5 for isoleucine (a discussion of other scales can be found in reference [38]). The overall degree of hydrophobicity for a single protein is given by where N is the number of residues. It is equal to −0.63 and −0.38 for 1GB1 and 1LIP respectively. Thus most of the residues, for both proteins, are more likely to stay below the center of the interface.
The usage of equation (1) when the interfacial 'air' (or vapor) is replaced by oil is justified by the following argument. Suppose a protein approaches the interface from the water side. Its hydropathic residues prefer to get out of water and join the oil phase. Thus hydropathy is conceptually identical to oil-liking. We do not know the actual values of q i when the protein is surrounded by oil, but expect that their absolute values do not change much and the equation applies approximately. On the other hand, when the protein approaches the interface from the oil side, it will be driven to it by oil-hating, i.e. hydrophilic, residues. Thus the signs of q i (and possibly the magnitudes) should be reversed. Here, however, we consider only the situations, in which the protein is initially fully in the water phase.
The parameter A controls the strength of the hydropathic forces. We chose A to be equal to 10 so that when a protein that comes to the interface, it usually stays at it instead of getting depinned. The width parameter W distinguishes between the two kinds of the interfaces that we study, and, in principle, A may also be distinct. However, Danov and Kralchevsky [39] have argued that the free energy of surfactant adsorption for oil-water is merely 10% higher than that for air-water. The free energy was determined by considering surfactants (nhexadecanes) moving from the bulk to the interface. This result can be captured by taking A for the oil-water case to be the same as for air-water. Figure 2 shows the single-residue force for both interfaces and the corresponding potential energy which is given by where erf denotes the error function. For the identical values of A, the energies at the center of the interfaces are the same. Figure 3 illustrates the effect of the hydropathy forces on 1GB1 and 1LIP [identical to 1LIP(4)] at the two kinds of the interfaces.Throughout the paper, the hydrophobic residues are indicated in red. It is seen that the oil-water conformations are more extended in the vertical direction compared to the air-water case.
For curved interfaces, formula (1) still applies but z needs to be replaced by the normal coordinate z along the direction that connects the center of curvature with the interfacial point that is closest to the residue under consideration. The direction of the force is normal. Notice that now the forces act on the residues in multitude of directions. In our studies, we take the sphere of radius R so that its center, placed at (0, 0, 0), serves as a center of curvature for all points on the surface of the sphere.
We consider four spherical situations: water droplet in air, denoted as w/a, air bubble in water (a/w), water droplet in oil (w/o), and oil droplet in water (o/w). In each case, a protein or a number of proteins (up to 100) are initially placed in the water phase. In particular, when considering the bubbles, the proteins are placed outside the sphere. The outside region is then bounded by lateral walls at the distance of L x in the x-direction and L y in the y-direction. We take L x = L y = R + d = l. In most cases, d = 50 Å, but d is an adjustable parameter, when studying systems of varying densities. There is a similar bouncing wall in the z-direction. We consider R ranging between 10 and 150 Å-larger values of R yield the similar adsorbed conformations as the flat case, at least for the proteins studied. When considering the droplets, the proteins are initially placed inside the sphere so no restraining boxes are needed. Examples of conformations of 100 1GB1 molecules for the w/a and a/w cases and R = 60 Å are shown in figure 4. Figure 5 shows conformations obtained for a single 1GB1 molecule at w/a and w/o interfaces with R = 15 Å.
In all cases, the initial conformation of each protein is native [the native structure of 1LIP(0) is taken to be the same as for 1LIP(4)] and the protein/proteins are released away from the interface. Whenever possible, and especially in the flat situation, the distance of the topmost atom in the initial state from the closest point on the interface is not smaller than 20 Å. In this location a protein is not yet affected by the interface. The diffusion towards the flat interface and equilibration takes place within a time up to of order 100 τ as in the flat situation [26]. This time depends on H and on the nature of the interface. To prevent the situation that proteins diffuse too deep into bulk water, we set a repulsive wall at z = −100 Å.
In the process of arriving at the interface, independent of whether the interface is flat or curved, the proteins get distorted. Here, we focus only on the stationary state obtained when the interface is occupied maximally for a given number In order to characterize the distortion of a single protein, we determine which residue penetrates the interface the most (see figure 5), calculate the moment of inertia around the axis that connects this residue with the center, and determine the corresponding radius of gyration around this axis. Its time averaged value is denoted as R z . The averaging interval is 1000 τ after achieving the stationary state. Figure 6 shows the density profiles for the flat and curved cases when the protein is 1GB1. The top panels correspond to the dense situation that arises when the lateral length l is 120 Å. These profiles are plotted against the z coordinate. It should be noted that the inter-proteinic interaction, as regulated by the parameter λ, does not affect the profiles in any significant manner except for the first layer. Most of the proteins gather close to z = 0 which shows as a major peak. However, in the air-water situation this peak is closer to the center of the interface at z = 0 Å than in the oil-water case. This is because W for the oil-water interface is twice as big than for the air-water one. There is also a secondary peak at either 30 Å into the bulk or 40 Å for the air-water and oil-water cases respectively. The corresponding snapshots are shown in figures 7 and 8 for l = 120 and 180. We have determined that the proteins in the first layer are fully alligned whereas they point in random directions when away from the first layer. As a result, there is no clear boundary between the second and the third layers, see figure 7. At lower densities (l > 180), we observe that the second layer disappears.

Results
In figure 6, the three lowest panels on the left correspond to w/a, and on the right-to a/w. The values of R are indicated. The profiles are radially symmetric. The histograms are plotted against the radial distance r from which R is subtracted. Thus the negative values of r − R correspond to the inside of the droplet and the positive ones to the outside of the bubble. It is seen that the bubble brings the proteins closer to the center of the interface than the corresponding droplet. This is because the crowdedness level of residues near the surface of a bubble is smaller than for a droplet, especially for a small  R. This is because the hydrophobic residues pointing from the outside of the sphere (proteins making a 'fan') generate fewer steric constraints than the residues pointing from the inside. We observe that the profiles shift with R, approaching the flat limit gradually (large R). Figure 9 characterizes the protein overall shape, as measured by R z (the effective adsorption radius of the protein in the plane perpendicular to the normal coordinate z ) as a function of R for the w/a, a/w, w/o, and o/w situations and for the three proteins studied. The R-dependence is fairly strong, and  The effective radius R z as a function of R for the three proteins studied, as indicated. The black lines, solid and broken, correspond to the w/a and a/w cases respectively. The blue lines, solid and broken,-to the o/w and w/o cases respectively. not necessarily monotonic, for R's not exceeding about 100 Å. We observe that the dependencies are not symmetric with respect to the replacement of w/a by a/w and of w/o by o/w. These observations agree with the results present in figure 6. For example, the density profiles get closer to that of flat interface (shown in the left bottom panel) as R increases in both w/a and a/w cases. There are several reasons for the lack of the symmetry. The first is that the pinning centers are usually distinct, as shown in figure 10 for the three proteins studied. For 1GB1 there is a good level of symmetry between w/a and a/w but this is not so in the other cases. The pinning centers are defined as residues that cross the center of the interface on going from the water side. Figure 10 shows the probabilities of doing so in the stationary state. The second is related to the way we calculate the effective radius R z -it is done by projecting all residues to the spherical interface with respect to the water phase. Thus, R z is also affected by the location of the residues. In cases of w/a and w/o, the hydrophobic residues are fanning out of the sphere at small R (such as R < 35 Å for 1GB1) as shown in figure 5. This results in a larger R z . As R increases further, this effect is reduced, which leads to a decrease in R z . For a/w and o/w, no fanning out of the proteins is observed, because then the hydrophobic residues are located inside of the sphere.
It is expected that for large values of R, one should obtain results corresponding to the flat interface (the data points placed at the right hand edges of the panels), i.e., R z ∼ = R z . In most cases this is indeed so when R = 150 Å. However, this does not seem to be valid for w/o in the case of 1GB1. As one can see from figure 9 for 1GB1, the value of R z for w/o at R = 150 Å is larger than that for the flat case. We find that, in this case, the protein can be pinned at the interface by different well separated segments (see figure 10). These segments remain still well separated even when R = 150 Å. This phenomenon is not observed for w/a. We propose that the larger width W for w/o makes the effective curvature of the sphere smaller and drives the separation of the pinned segments of 1GB1, resulting in larger R z .

Conclusions
In our previous works, we have studied proteins at flat air-water interfaces within a (native) structure coarse-grained model. Here, we have generalized this empirical approach to be applicable to the oil-water systems and to the curved interfaces, such as occuring in bubbles. Our simple model is primarily meant to elucidate the essential physics of the behavior of proteins at fluid-fluid interfaces. In the flat case, the protein density profiles show the existence of the second layer if the concentration of the proteins is sufficiently high. Unlike the first layer located right at the interface, the second layer is not ordered orientationally. The density profiles move with the radius of curvature and are not sensitive to the nature of the interactions between the proteins. However, due to steric constraints, the aqueous proteins come closer to the surface of a bubble than to the surface of a similarly curved droplet, i.e. with the nominally same surface area. The bubbles also adsorb more proteins and could thus acquire a stronger stabilization.
The curved interfaces generate different deformations of the proteins than the flat ones and the differences show for R smaller than about 150 Å. Our results seem to echo the findings of studies of protein adsorption on silica nanoparticles [40]-the analogs of our bubbles-that conformational changes depend on the particles' curvature. Our model may perhaps be applied to situations that involve interactions with solids, provided one derives rules for finding which segments of the proteins get anchored at the solid surface. This task most likely requires all-atom simulations and is not easy, especially in the case of gold [41]. The case of the curved membranes is also complicated because the cellular membranes are typically densely strewn with built-in proteins. This is not so with synthetic membranes.
One possible further generalization of our model is to consider the radius to be time dependent. In particular this feature could mimicks droplet evaporation, allowing for studies of protein flows under such conditions.